bp_bin/denoise_seq

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2012 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  22
  23 # This program is part of the Biopieces framework (www.biopieces.org).
  24
  25 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  26
  27 # Denoises sequences with quality scores in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31 require 'pp'
  32 require 'maasha/biopieces'
  33 require 'maasha/seq'
  34 require 'maasha/fastq'
  35 require 'maasha/fasta'
  36 require 'maasha/align'
  37 require 'maasha/usearch'
  38
  39 class Seq
  40   ALPH_DNA  = %w{A T C G}
  41   ALPH_AMBI = %w{A T C G M R W S Y K V H D B N}
  42
  43   BIT_INDEL = 0
  44   BIT_A = 1 << 0
  45   BIT_T = 1 << 1
  46   BIT_C = 1 << 2
  47   BIT_G = 1 << 3
  48
  49   BIT_M = BIT_A | BIT_C
  50   BIT_R = BIT_A | BIT_G
  51   BIT_W = BIT_A | BIT_T
  52   BIT_S = BIT_C | BIT_G
  53   BIT_Y = BIT_C | BIT_T
  54   BIT_K = BIT_G | BIT_T
  55   BIT_V = BIT_A | BIT_C | BIT_G
  56   BIT_H = BIT_A | BIT_C | BIT_T
  57   BIT_D = BIT_A | BIT_G | BIT_T
  58   BIT_B = BIT_C | BIT_G | BIT_T
  59   BIT_N = BIT_G | BIT_A | BIT_T | BIT_C
  60
  61   BITMAP = [
  62     BIT_A,
  63     BIT_T,
  64     BIT_C,
  65     BIT_G,
  66     BIT_M,
  67     BIT_R,
  68     BIT_W,
  69     BIT_S,
  70     BIT_Y,
  71     BIT_K,
  72     BIT_V,
  73     BIT_H,
  74     BIT_D,
  75     BIT_B,
  76     BIT_N
  77   ]
  78
  79   TR_NUC = "-" + ALPH_AMBI.join("").downcase
  80   TR_HEX = [BIT_INDEL].pack("C") + BITMAP.pack("C*")
  81
  82   def to_na
  83     entry      = Seq.new
  84     entry.seq  = NArray.to_na(self.seq.downcase.tr(TR_NUC, TR_HEX), "byte")
  85     entry.qual = NArray.to_na(self.qual, "byte") - Seq::SCORE_BASE if self.qual
  86     entry
  87   end
  88 end
  89
  90 class Align
  91   def to_na
  92     cols = self.length
  93     rows = self.members
  94
  95     na_seq  = NArray.byte(cols, rows)
  96     na_qual = NArray.byte(cols, rows)
  97
  98     self.entries.each_with_index do |entry, i|
  99       na_entry = entry.to_na
 100       na_seq[true, i]  = na_entry.seq
 101       na_qual[true, i] = na_entry.qual
 102     end
 103
 104     return na_seq, na_qual
 105   end
 106 end
 107
 108 class Denoise
 109   attr_reader :align
 110
 111   ROW_A = 0
 112   ROW_T = 1
 113   ROW_C = 2
 114   ROW_G = 3
 115
 116   def initialize(align, options)
 117     @align   = align
 118     @options = options
 119     @cols    = align.length
 120     @rows    = align.members
 121     @na_seq, @na_qual = align.to_na
 122     @na_rescored = nil
 123   end
 124
 125   def denoise_sequences
 126     freq = NArray.int(@cols, 4)
 127
 128     freq[true, ROW_A] = (@na_seq & Seq::BIT_A > 0).to_type("int").sum(1)
 129     freq[true, ROW_T] = (@na_seq & Seq::BIT_T > 0).to_type("int").sum(1)
 130     freq[true, ROW_C] = (@na_seq & Seq::BIT_C > 0).to_type("int").sum(1)
 131     freq[true, ROW_G] = (@na_seq & Seq::BIT_G > 0).to_type("int").sum(1)
 132
 133     mask_freq = freq.eq freq.max(1)
 134
 135     mask_freq[true, ROW_A] *= Seq::BIT_A
 136     mask_freq[true, ROW_T] *= Seq::BIT_T
 137     mask_freq[true, ROW_C] *= Seq::BIT_C
 138     mask_freq[true, ROW_G] *= Seq::BIT_G
 139
 140     mask_replace = mask_freq.max(1)
 141
 142     mask_bad = @na_qual <= @options[:quality_min]
 143
 144     @na_rescored = mask_bad.dup
 145
 146     new_values = mask_replace * mask_bad
 147     old_values = mask_bad.eq(0) * @na_seq
 148
 149     old_scores = @na_qual * mask_bad.eq(0)
 150
 151     sum   = old_scores.to_type("float").sum(1)
 152     count = (old_scores > 0).to_type("int").sum(1)
 153     mean  = (sum / count).to_type("byte")
 154
 155     new_scores = mask_bad * mean
 156
 157     @na_seq  = new_values + old_values
 158     @na_qual = new_scores + old_scores
 159
 160     self
 161   end
 162
 163   # Method that lowercase residues that have been removed during
 164   # the determination of the consensus sequence.
 165   def mask_sequences
 166     @align.each_with_index do |entry, i|
 167       entry.qual = (@na_qual[true, i] + Seq::SCORE_BASE).to_s
 168
 169       j = 0
 170
 171       while entry.seq[j] do
 172         if @na_rescored[j, i] == 0
 173           entry.seq[j] = entry.seq[j].upcase
 174         else
 175           entry.seq[j] = entry.seq[j].downcase
 176         end
 177
 178         j += 1
 179       end
 180     end
 181   end
 182 end
 183
 184 casts = []
 185 casts << {:long=>'cluster_ident', :short=>'i', :type=>'float', :mandatory=>true,  :default=>0.97, :allowed=>nil, :disallowed=>nil}
 186 casts << {:long=>'sequence_min',  :short=>'s', :type=>'uint',  :mandatory=>true,  :default=>1,    :allowed=>nil, :disallowed=>"0"}
 187 casts << {:long=>'residue_min',   :short=>'r', :type=>'float', :mandatory=>true,  :default=>0.3,  :allowed=>nil, :disallowed=>nil}
 188 casts << {:long=>'gap_max',       :short=>'g', :type=>'float', :mandatory=>true,  :default=>0.4,  :allowed=>nil, :disallowed=>nil}
 189 casts << {:long=>'quality_min',   :short=>'q', :type=>'uint',  :mandatory=>true,  :default=>10,   :allowed=>nil, :disallowed=>nil}
 190 casts << {:long=>'quality_mean',  :short=>'Q', :type=>'uint',  :mandatory=>true,  :default=>15,   :allowed=>nil, :disallowed=>nil}
 191 casts << {:long=>'cpus',          :short=>'C', :type=>'uint',  :mandatory=>true,  :default=>1,    :allowed=>nil, :disallowed=>"0"}
 192
 193 options          = Biopieces.options_parse(ARGV, casts)
 194 tmpdir           = Biopieces.mktmpdir
 195 fastq_file       = File.join(tmpdir, "test.fq")
 196 fasta_file       = File.join(tmpdir, "test.fna")
 197 fasta_file_align = File.join(tmpdir, "test.aln.fna")
 198
 199 options[:identity] = options[:cluster_ident]
 200 options[:msa]      = true
 201
 202 def alignment_to_fastq(entries, index)
 203   entries.each do |entry|
 204     seq_name = entry.seq_name.sub(/^\*/, '')
 205     elem     = index.get(seq_name)   # disk based lookup
 206
 207     entry.seq_name = elem.seq_name
 208     entry.qual     = elem.qual
 209
 210     entry.seq.scan(/-+/) do |m|
 211       entry.qual = entry.qual[0 ... $`.length] + ('!' * m.length) + entry.qual[$`.length .. -1]
 212     end
 213   end
 214 end
 215
 216 index     = FastqIndex.new
 217 seq_count = 0
 218
 219 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 220   Fasta.open(fasta_file, "w") do |fasta_io|
 221     Fastq.open(fastq_file, "w") do |fastq_io|
 222       input.each_record do |record|
 223         if record[:SEQ] and record[:SCORES]
 224           entry = Seq.new_bp(record)
 225           orig_name = entry.seq_name.dup
 226           entry.seq_name = seq_count.to_s
 227
 228           fasta_io.puts entry.to_fasta
 229           fastq_io.puts entry.to_fastq
 230
 231           index.add(entry, orig_name)
 232
 233           seq_count += 1
 234         else
 235           output.puts record
 236         end
 237       end
 238     end
 239   end
 240
 241   fastq_io  = File.open(fastq_file, "r")
 242   index.ios = fastq_io
 243
 244   uc = Usearch.new(fasta_file, fasta_file_align, options)
 245   uc.sortbylength
 246   uc.cluster_smallmem
 247
 248   uc.each_alignment do |align|
 249     align.options = options
 250     alignment_to_fastq(align.entries, index)
 251
 252     dn = Denoise.new(align, options)
 253     dn.denoise_sequences
 254     dn.mask_sequences
 255
 256     if options[:verbose]
 257       puts dn.align
 258     else
 259       dn.align.each do |entry|
 260         output.puts entry.to_bp
 261       end
 262     end
 263   end
 264 end
 265
 266 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 267