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   = 20
  33 INDEL_RATIO = 0.5
  34 QUAL_MIN    = 15
  35 QUAL_BASE   = 64
  36 ALPH_DNA    = %w{A T C G}
  37 ALPH_AMBI   = %w{A T C G M R W S Y K V H D B N}
  38
  39 BIT_INDEL = 0
  40 BIT_A = 1 << 0
  41 BIT_T = 1 << 1
  42 BIT_C = 1 << 2
  43 BIT_G = 1 << 3
  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 # Adding an initialize method to the existing Fasta class.
  83 class Fasta
  84   def initialize(io)
  85     @io = io
  86   end
  87 end
  88
  89 # Class for aligning sequences.
  90 class Align
  91   # Class method to align sequences in a list of Seq objects and
  92   # return these as a new list of Seq objects.
  93   def self.muscle(entries)
  94     cmd     = "muscle"
  95     aligned = []
  96
  97     Open3.popen3(cmd) do |stdin, stdout, stderr|
  98       entries.each { |entry| stdin.puts entry.to_fasta }
  99       stdin.close
 100
 101       fa = Fasta.new(stdout)
 102
 103       fa.each do |align|
 104         if block_given?
 105           yield align
 106         else
 107           aligned << align
 108         end
 109       end
 110     end
 111
 112     return aligned unless block_given?
 113   end
 114
 115   # Class method to create a consensus sequence from a list of Seq objects and
 116   # return a new Seq object with the consensus.
 117   def self.consensus(entries)
 118     cons = Consensus.new(entries.first.length)
 119
 120     entries.each { |entry| cons.add(entry) }
 121
 122     cons.to_seq
 123   end
 124
 125   class Consensus
 126     # Method to initialize a Consensus object given a Seq object, which is added
 127     # to the Consensus.
 128     def initialize(length)
 129       @length = length
 130
 131       @count    = NArray.int(@length)
 132       @freq     = NArray.int(@length, ALPH_DNA.size)
 133       @qual     = NArray.int(@length, ALPH_DNA.size)
 134       @mask_A   = NArray.byte(@length).fill!(BIT_A)
 135       @mask_T   = NArray.byte(@length).fill!(BIT_T)
 136       @mask_C   = NArray.byte(@length).fill!(BIT_C)
 137       @mask_G   = NArray.byte(@length).fill!(BIT_G)
 138       @has_qual = false
 139     end
 140
 141     # Method to add a Seq entry to a Consensus object.
 142     def add(entry)
 143       seq  = NArray.to_na(entry.seq.downcase.tr(TR_NUC, TR_HEX), "byte")
 144
 145       @count += seq > 0
 146
 147       mask_A = (seq & @mask_A) > 0
 148       mask_T = (seq & @mask_T) > 0
 149       mask_C = (seq & @mask_C) > 0
 150       mask_G = (seq & @mask_G) > 0
 151
 152       @freq[true, ROW_A] += mask_A
 153       @freq[true, ROW_T] += mask_T
 154       @freq[true, ROW_C] += mask_C
 155       @freq[true, ROW_G] += mask_G
 156
 157       unless entry.qual.nil?
 158         @has_qual = true
 159
 160         qual = NArray.to_na(entry.qual, "byte") - QUAL_BASE
 161
 162         @qual[true, ROW_A] += mask_A * qual
 163         @qual[true, ROW_T] += mask_T * qual
 164         @qual[true, ROW_C] += mask_C * qual
 165         @qual[true, ROW_G] += mask_G * qual
 166       end
 167     end
 168
 169     def to_seq
 170       #qual_filter
 171       new_seq      = Seq.new
 172       new_seq.seq  = consensus_seq
 173       new_seq.qual = consensus_qual if @has_qual
 174       new_seq.type = "dna"
 175
 176       new_seq
 177     end
 178
 179     private
 180
 181     # Method that calculates the sum for each column except the indel row and
 182     # returns the sum in an NArray.
 183     def freq_total
 184       @freq.sum(1)
 185     end
 186
 187     # Method that locates the maximum value for each column and
 188     # returns this in an NArray.
 189     def freq_max
 190       @freq.max(1)
 191     end
 192
 193     def qual_mean
 194       freq = @freq[true, [0 ... ALPH_DNA.size]]
 195       @qual / (freq + freq.eq(0))
 196     end
 197
 198     def freq_percent
 199       (@freq / freq_total.to_type("float") * 100).round
 200     end
 201
 202     # Method to calculate a consensus sequence from a Consensus object.
 203     def consensus_seq
 204       fp    = freq_percent > CONSENSUS
 205       cons  = NArray.byte(@length)
 206       cons |= fp[true, ROW_A] *= BIT_A
 207       cons |= fp[true, ROW_T] *= BIT_T
 208       cons |= fp[true, ROW_C] *= BIT_C
 209       cons |= fp[true, ROW_G] *= BIT_G
 210
 211       cons *= @count > @count.max * INDEL_RATIO
 212
 213       cons.to_s.tr!(TR_HEX, TR_NUC).upcase
 214     end
 215
 216     def consensus_qual
 217       q = (@qual / (@count + @count.eq(0))).max(1)
 218
 219       (q.to_type("byte") + QUAL_BASE).to_s
 220     end
 221
 222     def qual_filter
 223       filter = qual_mean > QUAL_MIN
 224
 225   #    pp @freq
 226   #    pp @qual
 227
 228   #    @freq *= filter
 229   #    @qual *= filter
 230
 231   #    pp @freq
 232   #    pp @qual
 233     end
 234   end
 235 end
 236
 237
 238
 239 __END__