code_ruby/lib/maasha/align.rb

   1 # Copyright (C) 2007-2011 Martin A. Hansen.
   2
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software
  15 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16
  17 # http://www.gnu.org/copyleft/gpl.html
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # This software is part of the Biopieces framework (www.biopieces.org).
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25 require 'pp'
  26 require 'open3'
  27 require 'narray'
  28 require 'maasha/fasta'
  29
  30 class AlignError < StandardError; end;
  31
  32 CONSENSUS   = 0.20
  33 FREQ_MIN    = 2
  34 QUAL_MIN    = 20
  35 ALPH_DNA    = %w{A T C G}
  36 ALPH_AMBI   = %w{A T C G M R W S Y K V H D B N}
  37
  38 BIT_INDEL = 0
  39 # BIT_A = 1 << 0
  40 # BIT_T = 1 << 1
  41 # BIT_C = 1 << 2
  42 # BIT_G = 1 << 3
  43
  44 BIT_M = BIT_A | BIT_C
  45 BIT_R = BIT_A | BIT_G
  46 BIT_W = BIT_A | BIT_T
  47 BIT_S = BIT_C | BIT_G
  48 BIT_Y = BIT_C | BIT_T
  49 BIT_K = BIT_G | BIT_T
  50 BIT_V = BIT_A | BIT_C | BIT_G
  51 BIT_H = BIT_A | BIT_C | BIT_T
  52 BIT_D = BIT_A | BIT_G | BIT_T
  53 BIT_B = BIT_C | BIT_G | BIT_T
  54 BIT_N = BIT_G | BIT_A | BIT_T | BIT_C
  55
  56 BITMAP = [
  57   BIT_A,
  58   BIT_T,
  59   BIT_C,
  60   BIT_G,
  61   BIT_M,
  62   BIT_R,
  63   BIT_W,
  64   BIT_S,
  65   BIT_Y,
  66   BIT_K,
  67   BIT_V,
  68   BIT_H,
  69   BIT_D,
  70   BIT_B,
  71   BIT_N
  72 ]
  73
  74 TR_NUC = "-"                   + ALPH_AMBI.join("").downcase
  75 TR_HEX = [BIT_INDEL].pack("C") + BITMAP.pack("C*")
  76
  77 ROW_A = 0
  78 ROW_T = 1
  79 ROW_C = 2
  80 ROW_G = 3
  81
  82 class Fasta
  83   def initialize(io)
  84     @io = io
  85   end
  86
  87   def close
  88     @io.close
  89   end
  90 end
  91
  92 # Class for aligning sequences.
  93 class Align
  94   attr_reader :entries
  95
  96   # Class method to align sequences in a list of Seq objects and
  97   # return these as a new list of Seq objects.
  98   def self.muscle(entries, verbose = false)
  99     result   = []
 100     index    = {}
 101
 102     Open3.popen3("muscle", "-quiet") do |stdin, stdout, stderr|
 103       entries.each do |entry|
 104         raise AlignError, "Duplicate sequence name: #{entry.seq_name}" if index.has_key? entry.seq_name
 105
 106         index[entry.seq_name] = entry
 107
 108         stdin.puts entry.to_fasta
 109       end
 110
 111       stdin.close
 112
 113       stderr.each_line { |line| $stderr.puts line } if verbose
 114
 115       aligned_entries = Fasta.new(stdout)
 116
 117       aligned_entries.each do |fa_entry|
 118         fq_entry = index[fa_entry.seq_name]
 119
 120         fa_entry.seq.scan(/-+/) do |m|
 121           fq_entry.seq  = fq_entry.seq[0 ... $`.length]  + ('-' * m.length) + fq_entry.seq[$`.length .. -1]
 122           fq_entry.qual = fq_entry.qual[0 ... $`.length] + ('@' * m.length) + fq_entry.qual[$`.length .. -1] unless fq_entry.qual.nil?
 123         end
 124
 125         result << fq_entry
 126       end
 127     end
 128
 129     self.new(result)
 130   end
 131
 132   # Method to initialize an Align object with a list of aligned Seq objects.
 133   def initialize(entries)
 134     @entries = entries
 135   end
 136
 137   # Method that returns the length of the alignment.
 138   def length
 139     @entries.first.length
 140   end
 141
 142   # Method that returns the number of members or sequences in the alignment.
 143   def members
 144     @entries.size
 145   end
 146
 147   # Method to create a consensus sequence from an Align object and
 148   # return a new Seq object with the consensus.
 149   def consensus
 150     cons = Consensus.new(@entries)
 151     cons.to_seq
 152   end
 153
 154   # Method to pretty print an alignment from an Align object.
 155   def to_s
 156     cons = self.consensus
 157     cons.mask_seq_soft!(QUAL_MIN) unless cons.qual.nil?
 158
 159     max_name = @entries.group_by { |entry| entry.seq_name.length }.max.first
 160
 161     @entries.each do |entry|
 162       entry.mask_seq_soft!(QUAL_MIN) unless entry.qual.nil?
 163       puts entry.seq_name + (" " * (max_name + 3 - entry.seq_name.length )) + entry.seq
 164     end
 165
 166     output = ""
 167     output << " " * (max_name + 3) + cons.seq + $/
 168     output << " " * (max_name + 3) + cons.qual.tr("[@-h]", "           ..........ooooooooooOOOOOOOOOO") unless cons.qual.nil?
 169     output
 170   end
 171
 172   private
 173
 174   class Consensus
 175     # Method to initialize a Consensus object given a list of aligned Seq object.
 176     def initialize(entries)
 177       @entries = entries
 178       @cols    = entries.first.length
 179       @rows    = entries.size
 180
 181       @has_qual = entries.first.qual.nil? ? false : true
 182
 183       @na_seq  = NArray.byte(@cols, @rows)
 184       @na_qual = NArray.byte(@cols, @rows) if @has_qual
 185
 186       na_add_entries
 187
 188       if @has_qual
 189         m1 = mask_high_qual
 190         m2 = mask_conserved_columns
 191         m3 = mask_conserved_residues
 192
 193         @na_seq  *= m3
 194         @na_qual *= m3
 195
 196         @na_seq  *= (m1 | m2)
 197         @na_qual *= (m1 | m2)
 198       end
 199     end
 200
 201     def to_seq
 202       new_seq      = Seq.new
 203       new_seq.seq  = consensus_seq
 204       new_seq.qual = consensus_qual if @has_qual
 205       new_seq.type = "dna"
 206
 207       new_seq
 208     end
 209
 210     private
 211
 212     def na_add_entries
 213       @entries.each_with_index do |entry, i|
 214         @na_seq[true, i]  = NArray.to_na(entry.seq.downcase.tr(TR_NUC, TR_HEX), "byte")
 215         @na_qual[true, i] = NArray.to_na(entry.qual, "byte") - SCORE_ILLUMINA if @has_qual
 216       end
 217     end
 218
 219     def mask_high_qual
 220       @na_qual > QUAL_MIN
 221     end
 222
 223     def mask_conserved_columns
 224       @na_seq * (@na_seq - @na_seq[true, 0]).sum(1).eq(0) > 0
 225     end
 226
 227     def mask_conserved_residues
 228       factor = (1 / @rows.to_f)
 229       mask_A = @na_seq & BIT_A > 0
 230       mask_T = @na_seq & BIT_T > 0
 231       mask_C = @na_seq & BIT_C > 0
 232       mask_G = @na_seq & BIT_G > 0
 233
 234       mask_A = (mask_A * mask_A.sum(1)).to_f * factor > CONSENSUS
 235       mask_T = (mask_T * mask_T.sum(1)).to_f * factor > CONSENSUS
 236       mask_C = (mask_C * mask_C.sum(1)).to_f * factor > CONSENSUS
 237       mask_G = (mask_G * mask_G.sum(1)).to_f * factor > CONSENSUS
 238
 239       mask_A | mask_T | mask_C | mask_G
 240     end
 241
 242     # Method to calculate a consensus sequence from a Consensus object.
 243     def consensus_seq
 244       cons  = NArray.byte(@cols)
 245       cons |= (@na_seq & BIT_A).max(1)
 246       cons |= (@na_seq & BIT_T).max(1)
 247       cons |= (@na_seq & BIT_C).max(1)
 248       cons |= (@na_seq & BIT_G).max(1)
 249
 250       cons.to_s.tr!(TR_HEX, TR_NUC).upcase
 251     end
 252
 253     # Method to calculate a consensus quality from a Consensus object.
 254     def consensus_qual
 255       (@na_qual.mean(1).round.to_type("byte") + SCORE_ILLUMINA).to_s
 256     end
 257   end
 258 end
 259
 260
 261
 262 __END__
 263
 264 cons |= ((@na_seq & BIT_A > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_A
 265 cons |= ((@na_seq & BIT_T > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_T
 266 cons |= ((@na_seq & BIT_C > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_C
 267 cons |= ((@na_seq & BIT_G > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_G