forester/ruby/scripts/bioruby_examples/msa_1.rb

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