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 CONS_MIN    = 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_accessor :options
  95   attr_reader   :entries
  96
  97   # Class method to align sequences in a list of Seq objects and
  98   # return these as a new list of Seq objects.
  99   def self.muscle(entries, verbose = false)
 100     result   = []
 101     index    = {}
 102
 103     Open3.popen3("muscle", "-quiet") do |stdin, stdout, stderr|
 104       entries.each do |entry|
 105         raise AlignError, "Duplicate sequence name: #{entry.seq_name}" if index.has_key? entry.seq_name
 106
 107         index[entry.seq_name] = entry
 108
 109         stdin.puts entry.to_fasta
 110       end
 111
 112       stdin.close
 113
 114       stderr.each_line { |line| $stderr.puts line } if verbose
 115
 116       aligned_entries = Fasta.new(stdout)
 117
 118       aligned_entries.each do |fa_entry|
 119         fq_entry = index[fa_entry.seq_name]
 120
 121         fa_entry.seq.scan(/-+/) do |m|
 122           fq_entry.seq  = fq_entry.seq[0 ... $`.length]  + ('-' * m.length) + fq_entry.seq[$`.length .. -1]
 123           fq_entry.qual = fq_entry.qual[0 ... $`.length] + ('@' * m.length) + fq_entry.qual[$`.length .. -1] unless fq_entry.qual.nil?
 124         end
 125
 126         result << fq_entry
 127       end
 128     end
 129
 130     self.new(result)
 131   end
 132
 133   # Method to initialize an Align object with a list of aligned Seq objects.
 134   def initialize(entries, options = {})
 135     @entries = entries
 136     @options = options
 137   end
 138
 139   # Method that returns the length of the alignment.
 140   def length
 141     @entries.first.length
 142   end
 143
 144   # Method that returns the number of members or sequences in the alignment.
 145   def members
 146     @entries.size
 147   end
 148
 149   # Method to create a consensus sequence from an Align object and
 150   # return a new Seq object with the consensus.
 151   def consensus
 152     cons = Consensus.new(@entries, @options)
 153     cons.to_seq
 154   end
 155
 156   # Method to pretty print an alignment from an Align object.
 157   def to_s
 158     qual_min = @options[:quality_min] || QUAL_MIN
 159
 160     cons = self.consensus
 161     cons.mask_seq_soft!(qual_min) unless cons.qual.nil?
 162
 163     max_name = @entries.group_by { |entry| entry.seq_name.length }.max.first
 164
 165     @entries.each do |entry|
 166       entry.mask_seq_soft!(qual_min) unless entry.qual.nil?
 167       puts entry.seq_name + (" " * (max_name + 3 - entry.seq_name.length )) + entry.seq
 168     end
 169
 170     output = ""
 171     output << " " * (max_name + 3) + cons.seq + $/
 172     output << " " * (max_name + 3) + cons.qual.tr("[@-h]", "           ..........ooooooooooOOOOOOOOOO") unless cons.qual.nil?
 173     output
 174   end
 175
 176   private
 177
 178   class Consensus
 179     # Method to initialize a Consensus object given a list of aligned Seq object.
 180     def initialize(entries, options)
 181       @entries = entries
 182       @options = options
 183
 184       @cons_min = options[:consensus_min] || CONS_MIN
 185       @freq_min = options[:frequency_min] || FREQ_MIN
 186       @qual_min = options[:quality_min]   || QUAL_MIN
 187
 188       @cols = entries.first.length
 189       @rows = entries.size
 190
 191       @has_qual = entries.first.qual.nil? ? false : true
 192
 193       @na_seq  = NArray.byte(@cols, @rows)
 194       @na_qual = NArray.byte(@cols, @rows) if @has_qual
 195
 196       na_add_entries
 197
 198       if @has_qual
 199         m1 = mask_high_qual
 200         m2 = mask_conserved_columns
 201         m3 = mask_conserved_residues
 202
 203         @na_seq  *= m3
 204         @na_qual *= m3
 205
 206         @na_seq  *= (m1 | m2)
 207         @na_qual *= (m1 | m2)
 208       end
 209     end
 210
 211     def to_seq
 212       new_seq      = Seq.new
 213       new_seq.seq  = consensus_seq
 214       new_seq.qual = consensus_qual if @has_qual
 215       new_seq.type = "dna"
 216
 217       new_seq
 218     end
 219
 220     private
 221
 222     def na_add_entries
 223       @entries.each_with_index do |entry, i|
 224         @na_seq[true, i]  = NArray.to_na(entry.seq.downcase.tr(TR_NUC, TR_HEX), "byte")
 225         @na_qual[true, i] = NArray.to_na(entry.qual, "byte") - SCORE_ILLUMINA if @has_qual
 226       end
 227     end
 228
 229     def mask_high_qual
 230       @na_qual > @qual_min
 231     end
 232
 233     def mask_conserved_columns
 234       @na_seq * (@na_seq - @na_seq[true, 0]).sum(1).eq(0) > 0
 235     end
 236
 237     def mask_conserved_residues
 238       factor = (1 / @rows.to_f)
 239       mask_A = @na_seq & BIT_A > 0
 240       mask_T = @na_seq & BIT_T > 0
 241       mask_C = @na_seq & BIT_C > 0
 242       mask_G = @na_seq & BIT_G > 0
 243
 244       mask_A = (mask_A * mask_A.sum(1)).to_f * factor > @cons_min
 245       mask_T = (mask_T * mask_T.sum(1)).to_f * factor > @cons_min
 246       mask_C = (mask_C * mask_C.sum(1)).to_f * factor > @cons_min
 247       mask_G = (mask_G * mask_G.sum(1)).to_f * factor > @cons_min
 248
 249       mask_A | mask_T | mask_C | mask_G
 250     end
 251
 252     # Method to calculate a consensus sequence from a Consensus object.
 253     def consensus_seq
 254       cons  = NArray.byte(@cols)
 255       cons |= (@na_seq & BIT_A).max(1)
 256       cons |= (@na_seq & BIT_T).max(1)
 257       cons |= (@na_seq & BIT_C).max(1)
 258       cons |= (@na_seq & BIT_G).max(1)
 259
 260       cons.to_s.tr!(TR_HEX, TR_NUC).upcase
 261     end
 262
 263     # Method to calculate a consensus quality from a Consensus object.
 264     def consensus_qual
 265       (@na_qual.mean(1).round.to_type("byte") + SCORE_ILLUMINA).to_s
 266     end
 267   end
 268 end
 269
 270
 271
 272 __END__
 273
 274 cons |= ((@na_seq & BIT_A > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_A
 275 cons |= ((@na_seq & BIT_T > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_T
 276 cons |= ((@na_seq & BIT_C > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_C
 277 cons |= ((@na_seq & BIT_G > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_G