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