bp_bin/pcr_seq

   1 #!/usr/bin/env ruby
   2
   3 # Copyright (C) 2007-2010 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 # Runs virtual PCR on sequences in the stream.
  28
  29 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  30
  31
  32 require 'maasha/biopieces'
  33 require 'maasha/fasta'
  34 require 'maasha/seq'
  35
  36 class Pcr
  37   def initialize(tmpdir, infile, options)
  38     @infile = infile
  39
  40     pattern    = Pattern.new(options[:forward], options[:reverse], options[:max_dist])
  41     @pat_files = pattern.save(tmpdir)
  42   end
  43
  44   # Run scan_for_matches using different pattern files on the subject
  45   # FASTA file and return a list of the result files.
  46   def run
  47     outfiles = []
  48
  49     @pat_files.each do |pat_file|
  50       outfile = pat_file.sub("pat", "fna")
  51
  52       command =  "scan_for_matches"
  53       # command << " -c"
  54       command << " -o 1"
  55       command << " #{pat_file}"
  56       command << " < #{@infile}"
  57       command << " > #{outfile}"
  58
  59       system(command)
  60       raise "Command failed: #{command}" unless $?.success?
  61
  62       outfiles << outfile
  63     end
  64
  65     outfiles
  66   end
  67 end
  68
  69 class Pattern
  70   attr_accessor :forward, :reverse
  71
  72   def initialize(forward, reverse, max_dist)
  73     @forward  = forward
  74     @reverse  = reverse
  75     @max_dist = max_dist
  76   end
  77
  78   # For each primer pair we need to check 4 possible
  79   # combinations that can give rise to PCR products:
  80   # - forward and reverse (the intended product)
  81   # - forward and revcomp forward
  82   # - revcomp reverse and reverse
  83   # - revcomp reverse and revcomp forward
  84   # Thus we create 4 pattern files and return
  85   # the file names.
  86   def save(tmpdir)
  87     forward = @forward
  88     reverse = @reverse
  89     revcomp_forward = revcomp(forward)
  90     revcomp_reverse = revcomp(reverse)
  91
  92     files = []
  93
  94     file = File.join(tmpdir, "forward_reverse.pat")
  95     self.forward = forward
  96     self.reverse = reverse
  97     save_pattern(file)
  98     files << file
  99
 100     file = File.join(tmpdir, "forward_forward.pat")
 101     self.forward = forward
 102     self.reverse = revcomp_forward
 103     save_pattern(file)
 104     files << file
 105
 106     file = File.join(tmpdir, "reverse_reverse.pat")
 107     self.forward = revcomp_reverse
 108     self.reverse = reverse
 109     save_pattern(file)
 110     files << file
 111
 112     file = File.join(tmpdir, "reverse_forward.pat")
 113     self.forward = revcomp_reverse
 114     self.reverse = revcomp_forward
 115     save_pattern(file)
 116     files << file
 117
 118     files
 119   end
 120
 121   private
 122
 123   # Method to output a pattern.
 124   def to_s
 125     "#{@forward} 1 ... #{@max_dist} #{@reverse}"
 126   end
 127
 128   # Save a pattern to file
 129   def save_pattern(file)
 130     File.open(file, mode="w") do |ios|
 131       ios.puts self
 132     end
 133   end
 134
 135   # Split a primer pattern in the form of ATCG[3,2,1] into
 136   # sequence and match descriptor, reverse complement the
 137   # primer and append the match descriptor: CGAT[3,2,1].
 138   def revcomp(pattern)
 139     if pattern.match(/^(\w+)(\[.+\])?/)
 140       primer     = $1
 141       descriptor = $2
 142     else
 143       raise "Failed splitting pattern: #{pattern}"
 144     end
 145
 146     seq      = Seq.new
 147     seq.seq  = primer
 148     seq.type = 'dna'
 149     seq.reverse!.complement!
 150
 151     descriptor ? seq.seq + descriptor : seq.seq
 152   end
 153 end
 154
 155 casts = []
 156 casts << {:long=>'forward',    :short=>'f', :type=>'string', :mandatory=>false, :default=>nil,  :allowed=>nil, :disallowed=>nil}
 157 casts << {:long=>'forward_rc', :short=>'F', :type=>'string', :mandatory=>false, :default=>nil,  :allowed=>nil, :disallowed=>nil}
 158 casts << {:long=>'reverse',    :short=>'r', :type=>'string', :mandatory=>false, :default=>nil,  :allowed=>nil, :disallowed=>nil}
 159 casts << {:long=>'reverse_rc', :short=>'R', :type=>'string', :mandatory=>false, :default=>nil,  :allowed=>nil, :disallowed=>nil}
 160 casts << {:long=>'max_dist',   :short=>'m', :type=>'uint',   :mandatory=>true, :default=>5000, :allowed=>nil, :disallowed=>"0"}
 161
 162 options = Biopieces.options_parse(ARGV, casts)
 163 tmpdir  = Biopieces.mktmpdir
 164 infile  = File.join(tmpdir, "in.fna")
 165
 166 if options[:forward_rc]
 167   options[:forward] = Seq.new("test", options[:forward_rc], 'dna').reverse.complement.seq
 168 end
 169
 170 if options[:reverse_rc]
 171   options[:reverse] = Seq.new("test", options[:reverse_rc], 'dna').reverse.complement.seq
 172 end
 173
 174 raise ArgumentError, "no adaptor specified" unless options[:forward] or options[:reverse]
 175 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 176   Fasta.open(infile, mode="w") do |ios|
 177     input.each_record do |record|
 178       output.puts record
 179
 180       if record[:SEQ]
 181         entry = Seq.new_bp(record)
 182         ios.puts entry.to_fasta
 183       end
 184     end
 185   end
 186
 187   outfiles = Pcr.new(tmpdir, infile, options).run
 188
 189   outfiles.each do |outfile|
 190     Fasta.open(outfile, mode="r") do |ios|
 191       ios.each do |entry|
 192         record = entry.to_bp
 193         record[:REC_TYPE] = "PCR"
 194         record[:STRAND]   = "+"
 195         record[:TYPE]     = File.basename(outfile).sub(".fna", "").upcase
 196         record[:SEQ_NAME].match(/(.+):\[(\d+),(\d+)\]$/)
 197         record[:SEQ_NAME] = $1
 198         record[:PCR_BEG]  = $2.to_i
 199         record[:PCR_END]  = $3.to_i
 200
 201         if record[:PCR_BEG] > record[:PCR_END]
 202           record[:PCR_BEG], record[:PCR_END] = record[:PCR_END], record[:PCR_BEG]
 203           record[:STRAND] = "-"
 204         end
 205
 206         output.puts record
 207       end
 208     end
 209   end
 210 end
 211
 212
 213 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 214
 215
 216 __END__