forester/ruby/scripts/bioruby_examples/msa_1.rb

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