forester/ruby/scripts/bioruby_examples/msa_1.rb

   1 require 'bio'
   2
   3 seq_ary = Array.new
   4 ff = Bio::FlatFile.auto('bcl2.fasta')
   5 ff.each_entry do |entry|
   6   seq_ary.push(entry)
   7   puts entry.entry_id          # prints the identifier of the entry
   8   puts entry.definition        # prints the definition of the entry
   9   puts entry.seq               # prints the sequence data of the entry
  10 end
  11
  12 # Creates a multiple sequence alignment (possibly unaligned) named
  13 # 'seqs' from array 'seq_ary'.
  14 seqs = Bio::Alignment.new( seq_ary )
  15
  16
  17 seqs.each { |seq| puts seq.to_s }
  18
  19
  20 puts seqs.consensus
  21
  22 # Writes multiple sequence alignment (possibly unaligned) 'seqs'
  23 # to a file in PHYLIP format.
  24 File.open('out0.phylip', 'w') do |f|
  25   f.write(seqs.output(:phylip))
  26 end
  27
  28 File.open('out0.fasta', 'w') do |f|
  29   f.write(seqs.output(:fasta))
  30 end
  31
  32 exit
  33 #############
  34
  35 # Reads in a FASTA-formatted multiple sequence alignment (which does
  36 # not have to be aligned, though) and stores its sequences in
  37 # array 'seq_ary'.
  38 seq_ary = Array.new
  39 fasta_seqs = Bio::Alignment::MultiFastaFormat.new(File.open('bcl2.fasta').read)
  40 fasta_seqs.entries.each do |seq|
  41   seq_ary.push( seq )
  42 end
  43
  44 # Creates a multiple sequence alignment (possibly unaligned) named
  45 # 'seqs' from array 'seq_ary'.
  46 seqs = Bio::Alignment.new( seq_ary )
  47 seqs.each { |seq| puts seq.to_s }
  48
  49
  50 puts seqs.consensus
  51
  52 # Writes multiple sequence alignment (possibly unaligned) 'seqs'
  53 # to a file in PHYLIP format.
  54 File.open('out1.phylip', 'w') do |f|
  55   f.write(seqs.output(:phylip))
  56 end
  57
  58 File.open('out1.fasta', 'w') do |f|
  59   f.write(seqs.output(:fasta))
  60 end
  61
  62 exit
  63 #################
  64
  65 #ff = Bio::FlatFile.new(Bio::FastaFormat, 'bcl2.fasta')
  66 #ff.each_entry do |f|
  67 #  puts "definition : " + f.definition
  68 #  puts "nalen      : " + f.nalen.to_s
  69 #  puts "naseq      : " + f.naseq
  70 #end
  71 #exit
  72
  73 seq_ary = Array.new
  74 Bio::FastaFormat.open('bcl2.fasta') do | file |
  75   file.each do |entry|
  76     puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
  77     # puts entry.definition         # Gets the complete fasta description line.
  78     #puts entry.seq                # Gets the actual sequence.
  79     seq_ary.push( entry )
  80   end
  81 end
  82 seqs =Bio::Alignment.new( seq_ary )
  83 seqs.each { |x| puts x }
  84 puts seqs.consensus
  85 exit
  86
  87
  88
  89 #will be obsolete!
  90 #seqs =Bio::Alignment.readfiles(File.open('bcl2.fasta'))
  91 #seqs.entries.each do |seq|
  92 #  puts seq.to_biosequence
  93 #end
  94
  95 #Bio::Alignment.
  96
  97
  98
  99 #exit
 100 #############
 101
 102 seqs = Bio::Alignment::MultiFastaFormat.new(File.open('bcl2.fasta').read)
 103 seqs.entries.each do |seq|
 104   puts seq.to_seq.output(:genbank)
 105 end
 106 puts
 107 puts
 108 puts :genbank
 109 puts seqs.entries[0].to_seq.output(:genbank)
 110 puts
 111 puts :fasta
 112 puts seqs.entries[0].to_seq.output(:fasta)
 113 puts
 114 puts :embl
 115 puts seqs.entries[0].to_seq.output(:embl)
 116 puts
 117 puts :raw
 118 puts seqs.entries[0].to_seq.output(:raw)
 119 puts
 120 puts :fasta_ncbi
 121 puts seqs.entries[0].to_seq.output(:fasta_ncbi)
 122 puts
 123 puts :fastq
 124 puts seqs.entries[0].to_seq.output(:fastq)
 125 puts
 126 puts :fastq_sanger
 127 puts seqs.entries[0].to_seq.output(:fastq_sanger)
 128 puts
 129 puts :fastq_solexa
 130 puts seqs.entries[0].to_seq.output(:fastq_solexa)
 131 puts
 132 puts :fastq_illumina
 133 puts seqs.entries[0].to_seq.output(:fastq_illumina)
 134 puts
 135 puts :fasta_numeric
 136 puts seqs.entries[0].to_seq.output(:fasta_numeric)
 137 puts
 138 puts :qual
 139 puts seqs.entries[0].to_seq.output(:qual)
 140 #exit
 141 ##############
 142
 143
 144 # Reads in a ClustalW formatted multiple sequence alignment
 145 # from a file named "infile_clustalw.aln" and stores it in 'report'.
 146 report = Bio::ClustalW::Report.new(File.read('infile_clustalw.aln'))
 147
 148 # Accesses the actual alignment.
 149 msa = report.alignment
 150
 151 # Goes through all sequences in 'msa' and prints the
 152 # actual molecular sequence.
 153 msa.each do |entry|
 154  # puts entry.seq
 155 end
 156
 157 ##############
 158
 159 DEFAULT_PARSER = Bio::Alignment::MultiFastaFormat
 160 puts DEFAULT_PARSER.to_s
 161
 162 #file = Bio::Alignment.readfiles('bcl2.fasta', Bio::Alignment::MultiFastaFormat)
 163 #file.each do |entry|
 164 #  puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
 165 #  puts entry.definition         # Gets the complete fasta description line.
 166 #  puts entry.seq                # Gets the actual sequence.
 167   #puts entry.aaseq.composition  # Gets the amino acid composition.
 168 #end
 169 #puts 'OK'
 170 #puts
 171
 172 file = Bio::FastaFormat.open('bcl2.fasta')
 173 file.each do |entry|
 174    puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
 175   puts entry.definition         # Gets the complete fasta description line.
 176   puts entry.seq                # Gets the actual sequence.
 177   # do something on each fasta sequence entry
 178 end
 179
 180
 181 ##############
 182
 183 # Creates a new file named "outfile.fasta" and writes
 184 # multiple sequence alignment 'msa' to it in fasta format.
 185 File.open('outfile.fasta', 'w') do |f|
 186   f.write(msa.output(:fasta))
 187 end
 188
 189 # Other formats
 190 File.open('outfile.clustal', 'w') do |f|
 191   f.write(msa.output(:clustal))
 192 end
 193 File.open('outfile.phylip', 'w') do |f|
 194   f.write(msa.output(:phylip))
 195 end
 196 File.open('outfile.phylipnon', 'w') do |f|
 197   f.write(msa.output(:phylipnon))
 198 end
 199 File.open('outfile.msf', 'w') do |f|
 200   f.write(msa.output(:msf))
 201 end
 202 File.open('outfile.molphy', 'w') do |f|
 203   f.write(msa.output(:molphy))
 204 end
 205
 206
 207
 208 #############
 209
 210 seq1 = Bio::Sequence.auto("gggggg")
 211
 212
 213 puts seq1.output(:fasta)
 214 #seq2 = Bio::Sequence::AA.new("ggggt")
 215 #seq3 = Bio::Sequence::AA.new("ggt")
 216
 217
 218
 219 seqs = ['MFQIPEFEPSEQEDSSSAER',
 220         'MGTPKQPSLAPAHALGLRKS',
 221         'PKQPSLAPAHALGLRKS',
 222         'MCSTSGCDLE']
 223
 224
 225 # MAFFT
 226 options = [ '--maxiterate', '1000', '--localpair' ]
 227 mafft = Bio::MAFFT.new('mafft', options )
 228 report = mafft.query_align( seqs)
 229
 230 # Accesses the actual alignment
 231 align = report.alignment
 232
 233 # Prints each sequence to the console.
 234 align.each { |s| puts s.to_s }
 235
 236
 237 puts 'MAFFT OK'
 238 puts
 239
 240 #clustalw
 241 clustalw = Bio::ClustalW.new('/home/zma/SOFTWARE/clustalw-2.1/src/clustalw2' )
 242 report = clustalw.query_align( seqs)
 243 #puts report.alignment.output_fasta.to_s
 244 report.alignment.each { |x| puts x.to_s }
 245 puts 'OK'
 246 puts
 247
 248 #muscle
 249 options = [ '-quiet', '-maxiters', '64' ]
 250 muscle = Bio::Muscle.new('/home/zma/SOFTWARE/muscle3.8.31/src/muscle', options )
 251 report = muscle.query_align( seqs)
 252 #puts report.alignment.output_fasta.to_s
 253 report.alignment.each { |x| puts x.to_s }
 254 puts 'OK'
 255 puts
 256
 257 file = Bio::FastaFormat.open('bcl2.fasta')
 258 file.each do |entry|
 259   puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
 260   puts entry.definition         # Gets the complete fasta description line.
 261   puts entry.seq                # Gets the actual sequence.
 262   puts entry.aaseq.composition  # Gets the amino acid composition.
 263 end
 264 puts 'OK'
 265 puts
 266
 267 Bio::FlatFile.auto('bcl2.fasta') do |ff|
 268   ff.each do |entry|
 269     puts entry.entry_id           # Gets the identifier, e.g. 'sp|O35147|BAD_RAT'.
 270     puts entry.definition         # Gets the complete fasta description line.
 271     puts entry.seq                # Gets the actual sequence.
 272     puts entry.aaseq.composition  # Gets the amino acid composition.
 273   end
 274 end
 275 puts 'OK'
 276 puts
 277
 278
 279
 280
 281
 282