forester/ruby/scripts/bioruby_examples/msa_1.rb

   1 require 'rubygems'
   2 require 'bio'
   3
   4 #############
   5
   6 # Reads in a ClustalW formatted multiple sequence alignment
   7 # from a file named "infile_clustalw.aln" and stores it in 'report'.
   8 report = Bio::ClustalW::Report.new(File.read('infile_clustalw.aln'))
   9
  10 # Accesses the actual alignment.
  11 msa = report.alignment
  12
  13 # Goes through all sequences in 'msa' and prints the
  14 # actual molecular sequence.
  15 msa.each do |entry|
  16   puts entry.seq
  17 end
  18
  19 ##############
  20
  21 DEFAULT_PARSER = Bio::Alignment::MultiFastaFormat
  22 puts DEFAULT_PARSER.to_s
  23
  24 #file = Bio::Alignment.readfiles('bcl2.fasta', Bio::Alignment::MultiFastaFormat)
  25 #file.each do |entry|
  26 #  puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
  27 #  puts entry.definition         # Gets the complete fasta description line.
  28 #  puts entry.seq                # Gets the actual sequence.
  29   #puts entry.aaseq.composition  # Gets the amino acid composition.
  30 #end
  31 #puts 'OK'
  32 #puts
  33
  34 file = Bio::FastaFormat.open('bcl2.fasta')
  35 file.each do |entry|
  36    puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
  37   puts entry.definition         # Gets the complete fasta description line.
  38   puts entry.seq                # Gets the actual sequence.
  39   # do something on each fasta sequence entry
  40 end
  41
  42 ##############
  43
  44 # Creates a new file named "outfile.fasta" and writes
  45 # multiple sequence alignment 'msa' to it in fasta format.
  46 File.open('outfile.fasta', 'w') do |f|
  47   f.write(msa.output(:fasta))
  48 end
  49
  50 # Other formats
  51 File.open('outfile.clustal', 'w') do |f|
  52   f.write(msa.output(:clustal))
  53 end
  54 File.open('outfile.phylip', 'w') do |f|
  55   f.write(msa.output(:phylip))
  56 end
  57 File.open('outfile.phylipnon', 'w') do |f|
  58   f.write(msa.output(:phylipnon))
  59 end
  60 File.open('outfile.msf', 'w') do |f|
  61   f.write(msa.output(:msf))
  62 end
  63 File.open('outfile.molphy', 'w') do |f|
  64   f.write(msa.output(:molphy))
  65 end
  66
  67
  68
  69 #############
  70
  71 seq1 = Bio::Sequence.auto("gggggg")
  72
  73
  74 puts seq1.output(:fasta)
  75 #seq2 = Bio::Sequence::AA.new("ggggt")
  76 #seq3 = Bio::Sequence::AA.new("ggt")
  77
  78
  79
  80 seqs = ['MFQIPEFEPSEQEDSSSAER',
  81         'MGTPKQPSLAPAHALGLRKS',
  82         'PKQPSLAPAHALGLRKS',
  83         'MCSTSGCDLE']
  84
  85
  86 # MAFFT
  87 options = [ '--maxiterate', '1000', '--localpair' ]
  88 mafft = Bio::MAFFT.new('mafft', options )
  89 report = mafft.query_align( seqs)
  90
  91 # Accesses the actual alignment
  92 align = report.alignment
  93
  94 # Prints each sequence to the console.
  95 align.each { |s| puts s.to_s }
  96
  97
  98 puts 'MAFFT OK'
  99 puts
 100
 101 #clustalw
 102 clustalw = Bio::ClustalW.new('/home/zma/SOFTWARE/clustalw-2.1/src/clustalw2' )
 103 report = clustalw.query_align( seqs)
 104 #puts report.alignment.output_fasta.to_s
 105 report.alignment.each { |x| puts x.to_s }
 106 puts 'OK'
 107 puts
 108
 109 #muscle
 110 options = [ '-quiet', '-maxiters', '64' ]
 111 muscle = Bio::Muscle.new('/home/zma/SOFTWARE/muscle3.8.31/src/muscle', options )
 112 report = muscle.query_align( seqs)
 113 #puts report.alignment.output_fasta.to_s
 114 report.alignment.each { |x| puts x.to_s }
 115 puts 'OK'
 116 puts
 117
 118 file = Bio::FastaFormat.open('bcl2.fasta')
 119 file.each do |entry|
 120   puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
 121   puts entry.definition         # Gets the complete fasta description line.
 122   puts entry.seq                # Gets the actual sequence.
 123   puts entry.aaseq.composition  # Gets the amino acid composition.
 124 end
 125 puts 'OK'
 126 puts
 127
 128 Bio::FlatFile.auto('bcl2.fasta') do |ff|
 129   ff.each do |entry|
 130     puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
 131     puts entry.definition         # Gets the complete fasta description line.
 132     puts entry.seq                # Gets the actual sequence.
 133     puts entry.aaseq.composition  # Gets the amino acid composition.
 134   end
 135 end
 136 puts 'OK'
 137 puts
 138
 139
 140
 141
 142
 143