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 << " #{pat_file}"
  55       command << " < #{@infile}"
  56       command << " > #{outfile}"
  57
  58       system(command)
  59       raise "Command failed: #{command}" unless $?.success?
  60
  61       outfiles << outfile
  62     end
  63
  64     outfiles
  65   end
  66 end
  67
  68 class Pattern
  69   attr_accessor :forward, :reverse
  70
  71   def initialize(forward, reverse, max_dist)
  72     @forward  = forward
  73     @reverse  = reverse
  74     @max_dist = max_dist
  75   end
  76
  77   # For each primer pair we need to check 4 possible
  78   # combinations that can give rise to PCR products:
  79   # - forward and reverse (the intended product)
  80   # - forward and revcomp forward
  81   # - revcomp reverse and reverse
  82   # - revcomp reverse and revcomp forward
  83   # Thus we create 4 pattern files and return
  84   # the file names.
  85   def save(tmpdir)
  86     forward = @forward
  87     reverse = @reverse
  88     revcomp_forward = revcomp(forward)
  89     revcomp_reverse = revcomp(reverse)
  90
  91     files = []
  92
  93     file = File.join(tmpdir, "forward_reverse.pat")
  94     self.forward = forward
  95     self.reverse = reverse
  96     save_pattern(file)
  97     files << file
  98
  99     file = File.join(tmpdir, "forward_forward.pat")
 100     self.forward = forward
 101     self.reverse = revcomp_forward
 102     save_pattern(file)
 103     files << file
 104
 105     file = File.join(tmpdir, "reverse_reverse.pat")
 106     self.forward = revcomp_reverse
 107     self.reverse = reverse
 108     save_pattern(file)
 109     files << file
 110
 111     file = File.join(tmpdir, "reverse_forward.pat")
 112     self.forward = revcomp_reverse
 113     self.reverse = revcomp_forward
 114     save_pattern(file)
 115     files << file
 116
 117     files
 118   end
 119
 120   private
 121
 122   # Method to output a pattern.
 123   def to_s
 124     "#{@forward} 1 ... #{@max_dist} #{@reverse}"
 125   end
 126
 127   # Save a pattern to file
 128   def save_pattern(file)
 129     File.open(file, mode="w") do |ios|
 130       ios.puts self
 131     end
 132   end
 133
 134   # Split a primer pattern in the form of ATCG[3,2,1] into
 135   # sequence and match descriptor, reverse complement the
 136   # primer and append the match descriptor: CGAT[3,2,1].
 137   def revcomp(pattern)
 138     if pattern.match(/^(\w+)(\[.+\])?/)
 139       primer     = $1
 140       descriptor = $2
 141     else
 142       raise "Failed splitting pattern: #{pattern}"
 143     end
 144
 145     seq      = Seq.new
 146     seq.seq  = primer
 147     seq.type = 'dna'
 148     seq.revcomp
 149
 150     descriptor ? seq.seq + descriptor : seq.seq
 151   end
 152 end
 153
 154 casts = []
 155 casts << {:long=>'forward',  :short=>'f', :type=>'string', :mandatory=>true, :default=>nil,  :allowed=>nil, :disallowed=>nil}
 156 casts << {:long=>'reverse',  :short=>'r', :type=>'string', :mandatory=>true, :default=>nil,  :allowed=>nil, :disallowed=>nil}
 157 casts << {:long=>'max_dist', :short=>'m', :type=>'uint',   :mandatory=>true, :default=>5000, :allowed=>nil, :disallowed=>"0"}
 158
 159 options = Biopieces.options_parse(ARGV, casts)
 160 tmpdir  = Biopieces.mktmpdir
 161 infile  = File.join(tmpdir, "in.fna")
 162
 163 Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
 164   Fasta.open(infile, mode="w") do |ios|
 165     input.each_record do |record|
 166       output.puts record
 167
 168       if record.has_key? :SEQ
 169         entry = Seq.new_bp(record)
 170         ios.puts entry.to_fasta
 171       end
 172     end
 173   end
 174
 175   outfiles = Pcr.new(tmpdir, infile, options).run
 176
 177   outfiles.each do |outfile|
 178     Fasta.open(outfile, mode="r") do |ios|
 179       ios.each do |entry|
 180         record = entry.to_bp
 181         record[:REC_TYPE] = "PCR"
 182         record[:STRAND]   = "+"
 183         record[:TYPE]     = File.basename(outfile).sub(".fna", "").upcase
 184         record[:SEQ_NAME].match(/(.+):\[(\d+),(\d+)\]$/)
 185         record[:SEQ_NAME] = $1
 186         record[:PCR_BEG]  = $2
 187         record[:PCR_END]  = $3
 188
 189         if record[:PCR_BEG] > record[:PCR_END]
 190           record[:PCR_BEG], record[:PCR_END] = record[:PCR_END], record[:PCR_BEG]
 191           record[:STRAND] = "-"
 192         end
 193
 194         output.puts record
 195       end
 196     end
 197   end
 198 end
 199
 200
 201 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 202
 203
 204 __END__