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     has_qual = false
  95     result   = []
  96     index    = {}
  97
  98     Open3.popen3("muscle") do |stdin, stdout, stderr|
  99       entries.each do |entry|
 100         index[entry.seq_name] = entry
 101
 102         stdin.puts entry.to_fasta
 103       end
 104
 105       stdin.close
 106
 107       aligned_entries = Fasta.new(stdout)
 108
 109       aligned_entries.each do |fa_entry|
 110         fq_entry = index[fa_entry.seq_name]
 111
 112         fa_entry.seq.scan(/-+/) do |m|
 113           fq_entry.seq  = fq_entry.seq[0 ... $`.length]  + ('-' * m.length) + fq_entry.seq[$`.length .. -1]
 114           fq_entry.qual = fq_entry.qual[0 ... $`.length] + ('@' * m.length) + fq_entry.qual[$`.length .. -1]
 115         end
 116
 117         result << fq_entry
 118       end
 119     end
 120
 121     result
 122   end
 123
 124   # Class method to create a consensus sequence from a list of Seq objects and
 125   # return a new Seq object with the consensus.
 126   def self.consensus(entries)
 127     cons = Consensus.new(entries.first.length)
 128
 129     entries.each { |entry| cons.add(entry) }
 130
 131     cons.to_seq
 132   end
 133
 134   class Consensus
 135     # Method to initialize a Consensus object given a Seq object, which is added
 136     # to the Consensus.
 137     def initialize(length)
 138       @length = length
 139
 140       @count    = NArray.int(@length)
 141       @freq     = NArray.int(@length, ALPH_DNA.size)
 142       @qual     = NArray.int(@length, ALPH_DNA.size)
 143       @mask_A   = NArray.byte(@length).fill!(BIT_A)
 144       @mask_T   = NArray.byte(@length).fill!(BIT_T)
 145       @mask_C   = NArray.byte(@length).fill!(BIT_C)
 146       @mask_G   = NArray.byte(@length).fill!(BIT_G)
 147       @has_qual = false
 148     end
 149
 150     # Method to add a Seq entry to a Consensus object.
 151     def add(entry)
 152       seq  = NArray.to_na(entry.seq.downcase.tr(TR_NUC, TR_HEX), "byte")
 153
 154       @count += seq > 0
 155
 156       mask_A = (seq & @mask_A) > 0
 157       mask_T = (seq & @mask_T) > 0
 158       mask_C = (seq & @mask_C) > 0
 159       mask_G = (seq & @mask_G) > 0
 160
 161       @freq[true, ROW_A] += mask_A
 162       @freq[true, ROW_T] += mask_T
 163       @freq[true, ROW_C] += mask_C
 164       @freq[true, ROW_G] += mask_G
 165
 166       unless entry.qual.nil?
 167         @has_qual = true
 168
 169         qual = NArray.to_na(entry.qual, "byte") - QUAL_BASE
 170
 171         @qual[true, ROW_A] += mask_A * qual
 172         @qual[true, ROW_T] += mask_T * qual
 173         @qual[true, ROW_C] += mask_C * qual
 174         @qual[true, ROW_G] += mask_G * qual
 175       end
 176     end
 177
 178     def to_seq
 179       #qual_filter
 180       new_seq      = Seq.new
 181       new_seq.seq  = consensus_seq
 182       new_seq.qual = consensus_qual if @has_qual
 183       new_seq.type = "dna"
 184
 185       new_seq
 186     end
 187
 188     private
 189
 190     # Method that calculates the sum for each column except the indel row and
 191     # returns the sum in an NArray.
 192     def freq_total
 193       @freq.sum(1)
 194     end
 195
 196     # Method that locates the maximum value for each column and
 197     # returns this in an NArray.
 198     def freq_max
 199       @freq.max(1)
 200     end
 201
 202     def qual_mean
 203       freq = @freq[true, [0 ... ALPH_DNA.size]]
 204       @qual / (freq + freq.eq(0))
 205     end
 206
 207     def freq_percent
 208       (@freq / freq_total.to_type("float") * 100).round
 209     end
 210
 211     # Method to calculate a consensus sequence from a Consensus object.
 212     def consensus_seq
 213       fp    = freq_percent > CONSENSUS
 214       cons  = NArray.byte(@length)
 215       cons |= fp[true, ROW_A] *= BIT_A
 216       cons |= fp[true, ROW_T] *= BIT_T
 217       cons |= fp[true, ROW_C] *= BIT_C
 218       cons |= fp[true, ROW_G] *= BIT_G
 219
 220       cons *= @count > @count.max * INDEL_RATIO
 221
 222       cons.to_s.tr!(TR_HEX, TR_NUC).upcase
 223     end
 224
 225     def consensus_qual
 226       q = (@qual / (@count + @count.eq(0))).max(1)
 227
 228       (q.to_type("byte") + QUAL_BASE).to_s
 229     end
 230
 231     def qual_filter
 232       filter = qual_mean > QUAL_MIN
 233
 234   #    pp @freq
 235   #    pp @qual
 236
 237   #    @freq *= filter
 238   #    @qual *= filter
 239
 240   #    pp @freq
 241   #    pp @qual
 242     end
 243   end
 244 end
 245
 246
 247
 248 __END__