# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-require 'biopieces'
-require 'fasta'
-require 'seq'
-require 'pp'
+require 'maasha/biopieces'
+require 'maasha/fasta'
+require 'maasha/seq'
class Pcr
def initialize(tmpdir, infile, options)
outfile = pat_file.sub("pat", "fna")
command = "scan_for_matches"
- command << " -c"
+ # command << " -c"
+ command << " -o 1"
command << " #{pat_file}"
command << " < #{@infile}"
command << " > #{outfile}"
class Pattern
attr_accessor :forward, :reverse
+ # Split a primer pattern in the form of ATCG[3,2,1] into
+ # sequence and match descriptor, reverse complement the
+ # primer and append the match descriptor: CGAT[3,2,1].
+ def self.revcomp(pattern)
+ if pattern.match(/^(\w+)(\[.+\])?/)
+ primer = $1
+ descriptor = $2
+ else
+ raise "Failed splitting pattern: #{pattern}"
+ end
+
+ seq = Seq.new
+ seq.seq = primer
+ seq.type = :dna
+ seq.reverse!.complement!
+
+ descriptor ? seq.seq + descriptor : seq.seq
+ end
+
def initialize(forward, reverse, max_dist)
@forward = forward
@reverse = reverse
def save(tmpdir)
forward = @forward
reverse = @reverse
- revcomp_forward = revcomp(forward)
- revcomp_reverse = revcomp(reverse)
+ revcomp_forward = Pattern.revcomp(forward)
+ revcomp_reverse = Pattern.revcomp(reverse)
files = []
ios.puts self
end
end
-
- # Split a primer pattern in the form of ATCG[3,2,1] into
- # sequence and match descriptor, reverse complement the
- # primer and append the match descriptor: CGAT[3,2,1].
- def revcomp(pattern)
- if pattern.match(/^(\w+)(\[.+\])?/)
- primer = $1
- descriptor = $2
- else
- raise "Failed splitting pattern: #{pattern}"
- end
-
- seq = Seq.new
- seq.seq = primer
- seq.type = 'dna'
- seq.revcomp
-
- descriptor ? seq.seq + descriptor : seq.seq
- end
end
casts = []
-casts << {:long=>'forward', :short=>'f', :type=>'string', :mandatory=>true, :default=>nil, :allowed=>nil, :disallowed=>nil}
-casts << {:long=>'reverse', :short=>'r', :type=>'string', :mandatory=>true, :default=>nil, :allowed=>nil, :disallowed=>nil}
-casts << {:long=>'max_dist', :short=>'m', :type=>'uint', :mandatory=>true, :default=>5000, :allowed=>nil, :disallowed=>"0"}
-
-bp = Biopieces.new
-
-options = bp.parse(ARGV, casts)
-tmpdir = bp.mktmpdir
+casts << {:long=>'forward', :short=>'f', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'forward_rc', :short=>'F', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'reverse', :short=>'r', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'reverse_rc', :short=>'R', :type=>'string', :mandatory=>false, :default=>nil, :allowed=>nil, :disallowed=>nil}
+casts << {:long=>'max_dist', :short=>'m', :type=>'uint', :mandatory=>true, :default=>5000, :allowed=>nil, :disallowed=>"0"}
+
+options = Biopieces.options_parse(ARGV, casts)
+tmpdir = Biopieces.mktmpdir
infile = File.join(tmpdir, "in.fna")
-Fasta.open(infile, mode="w") do |ios|
- bp.each_record do |record|
- bp.puts record
- ios.puts record
- end
+if options[:forward_rc]
+ options[:forward] = Pattern.revcomp(options[:forward_rc])
+end
+
+if options[:reverse_rc]
+ options[:reverse] = Pattern.revcomp(options[:reverse_rc])
end
-outfiles = Pcr.new(tmpdir, infile, options).run
-
-outfiles.each do |outfile|
- Fasta.open(outfile, mode="r") do |ios|
- ios.each do |entry|
- record = entry.to_bp
- record[:REC_TYPE] = "PCR"
- record[:STRAND] = "+"
- record[:TYPE] = File.basename(outfile).sub(".fna", "").upcase
- record[:SEQ_NAME].match(/(.+):\[(\d+),(\d+)\]$/)
- record[:SEQ_NAME] = $1
- record[:PCR_BEG] = $2
- record[:PCR_END] = $3
-
- if record[:PCR_BEG] > record[:PCR_END]
- record[:PCR_BEG], record[:PCR_END] = record[:PCR_END], record[:PCR_BEG]
- record[:STRAND] = "-"
+raise ArgumentError, "no adaptor specified" unless options[:forward] or options[:reverse]
+Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
+ Fasta.open(infile, mode="w") do |ios|
+ input.each_record do |record|
+ output.puts record
+
+ if record[:SEQ]
+ entry = Seq.new_bp(record)
+ ios.puts entry.to_fasta
end
+ end
+ end
- bp.puts record
+ outfiles = Pcr.new(tmpdir, infile, options).run
+
+ outfiles.each do |outfile|
+ Fasta.open(outfile, mode="r") do |ios|
+ ios.each do |entry|
+ record = entry.to_bp
+ record[:REC_TYPE] = "PCR"
+ record[:STRAND] = "+"
+ record[:TYPE] = File.basename(outfile).sub(".fna", "").upcase
+ record[:SEQ_NAME].match(/(.+):\[(\d+),(\d+)\]$/)
+ record[:SEQ_NAME] = $1
+ record[:PCR_BEG] = $2.to_i
+ record[:PCR_END] = $3.to_i
+
+ if record[:PCR_BEG] > record[:PCR_END]
+ record[:PCR_BEG], record[:PCR_END] = record[:PCR_END], record[:PCR_BEG]
+ record[:STRAND] = "-"
+ end
+
+ output.puts record
+ end
end
end
end