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 align = report.alignment
  12
  13 # Goes through all sequences in 'align' and prints the
  14 # actual molecular sequence.
  15 align.each do |entry|
  16   puts entry.seq
  17 end
  18
  19
  20 # Creates a new file named "outfile.fasta" and writes
  21 # multiple sequence alignment 'align' to it in fasta format.
  22 File.open('outfile.fasta', 'w') do |f|
  23   f.write(align.output(:fasta))
  24 end
  25
  26 # Creates a new file named "outfile.aln" and writes
  27 # multiple sequence alignment 'align' to it in clustal format.
  28 File.open('outfile.aln', 'w') do |f|
  29   f.write(align.output(:clustal))
  30 end
  31
  32 #############
  33
  34 seq1 = Bio::Sequence.auto("gggggg")
  35
  36
  37 puts seq1.output(:fasta)
  38 #seq2 = Bio::Sequence::AA.new("ggggt")
  39 #seq3 = Bio::Sequence::AA.new("ggt")
  40
  41
  42
  43 seqs = [ "MFQIPEFEPSEQEDSSSAER",
  44          "MGTPKQPSLAPAHALGLRKS",
  45          "PKQPSLAPAHALGLRKS",
  46          "MCSTSGCDLE" ]
  47
  48
  49 # MAFFT
  50 options = [ '--maxiterate', '1000', '--localpair' ]
  51 mafft = Bio::MAFFT.new('/home/zma/SOFTWARE/mafft-6.847-without-extensions/scripts/mafft', options )
  52 report = mafft.query_align( seqs)
  53
  54 # Accesses the actual alignment
  55 align = report.alignment
  56
  57 # Prints each sequence to the console.
  58 align.each { |s| puts s.to_s }
  59
  60
  61 puts 'MAFFT OK'
  62 puts
  63
  64 #clustalw
  65 clustalw = Bio::ClustalW.new('/home/zma/SOFTWARE/clustalw-2.1/src/clustalw2' )
  66 report = clustalw.query_align( seqs)
  67 #puts report.alignment.output_fasta.to_s
  68 report.alignment.each { |x| puts x.to_s }
  69 puts 'OK'
  70 puts
  71
  72 #muscle
  73 options = [ '-quiet', '-maxiters', '64' ]
  74 muscle = Bio::Muscle.new('/home/zma/SOFTWARE/muscle3.8.31/src/muscle', options )
  75 report = muscle.query_align( seqs)
  76 #puts report.alignment.output_fasta.to_s
  77 report.alignment.each { |x| puts x.to_s }
  78 puts 'OK'
  79 puts
  80
  81 file = Bio::FastaFormat.open('bcl2.fasta')
  82 file.each do |entry|
  83   puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
  84   puts entry.definition         # Gets the complete fasta description line.
  85   puts entry.seq                # Gets the actual sequence.
  86   puts entry.aaseq.composition  # Gets the amino acid composition.
  87 end
  88 puts 'OK'
  89 puts
  90
  91 Bio::FlatFile.auto('bcl2.fasta') do |ff|
  92   ff.each do |entry|
  93     puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
  94     puts entry.definition         # Gets the complete fasta description line.
  95     puts entry.seq                # Gets the actual sequence.
  96     puts entry.aaseq.composition  # Gets the amino acid composition.
  97   end
  98 end
  99 puts 'OK'
 100 puts
 101
 102
 103
 104
 105
 106