Permalink
Browse files

amplify: return exact sequence

  • Loading branch information...
1 parent dc6e3a1 commit df548639b42d54943f3093133e4d8ae7af31a063 greg committed Jun 7, 2011
Showing with 32 additions and 9 deletions.
  1. +1 −0 TODO
  2. +30 −8 amplify
  3. +1 −1 fastastack.c
View
1 TODO
@@ -1 +1,2 @@
++ Bug: fastastack seems to choke on blank lines in input
+ consider abstracting over fastahead, fastatail
View
38 amplify
@@ -2,24 +2,46 @@
#
# Usage: amplify <template.fasta> <left-primer.fasta> <right-primer.fasta>
#
-# Cuts out the portion of 'template.fasta' _roughly_ matching the region
-# amplified by left-primer and right-primer.
-#
-# To get the exact amplicon -- use your editor! Sorry!
-#
+# Performs a 'virtual PCR reaction' on template.fasta using left and right
+# primers; returns amplified sequence
+#
PREFIX=`dirname $0`;
TEMPLATE=$1;
LEFT_PRIMER=$2;
RIGHT_PRIMER=$3;
+ANNEAL="kalign -e100 -gpo 100"
+
RIGHT_PRIMER_LENGTH=`${PREFIX}/fastacount $RIGHT_PRIMER`;
LEFT_PRIMER_LENGTH=`${PREFIX}/fastacount $LEFT_PRIMER`;
-start=`cat $TEMPLATE $LEFT_PRIMER | kalign -e100 -gpo 100 2>/dev/null | ${PREFIX}/fastagap`;
+# do an initial cut with sloppy primer ends
+primer_align_start=`cat $LEFT_PRIMER $TEMPLATE | $ANNEAL 2>/dev/null \
+ | ${PREFIX}/fastagap`;
+
+end=$((`(cat $RIGHT_PRIMER | ${PREFIX}/fastacomplement; cat $TEMPLATE) \
+ | $ANNEAL 2>/dev/null | ${PREFIX}/fastagap` + RIGHT_PRIMER_LENGTH));
+
+tmpfile=`mktemp`;
+
+cat $TEMPLATE | ${PREFIX}/fastahead $end \
+ | ${PREFIX}/fastatail $primer_align_start > $tmpfile;
+
+tmpfile2=`mktemp`;
-end=$((`(cat $RIGHT_PRIMER | ${PREFIX}/fastacomplement; cat $TEMPLATE) | kalign -e100 2>/dev/null | ${PREFIX}/fastagap` + RIGHT_PRIMER_LENGTH));
+# snip left leader sequence
+(cat $LEFT_PRIMER;
+ cat $tmpfile $LEFT_PRIMER | $ANNEAL 2>/dev/null \
+ | ${PREFIX}/fastatail $LEFT_PRIMER_LENGTH ) \
+ | ${PREFIX}/fastacat > $tmpfile2;
-cat $TEMPLATE | ${PREFIX}/fastahead $end | ${PREFIX}/fastatail $start | ${PREFIX}/fastalint;
+# where does the right primer start?
+right_junction=`(${PREFIX}/fastacomplement $RIGHT_PRIMER; cat $tmpfile2) \
+ | $ANNEAL 2>/dev/null | ${PREFIX}/fastagap`;
+(cat $tmpfile2 | ${PREFIX}/fastahead $((right_junction - 1)); \
+ ${PREFIX}/fastacomplement $RIGHT_PRIMER) | ${PREFIX}/fastacat;
+rm $tmpfile;
+rm $tmpfile2;
View
@@ -37,7 +37,7 @@ void post_process()
comment_ptr = comment_ptr->next, sequence_ptr = sequence_ptr->next)
{
int i;
- /* FIXME report comment */
+ /* report comment */
char *comment = (char*)(comment_ptr->data);
for (i = 0; i < comment_max_width; ++i)
if (*comment) putchar(*comment++); else putchar('.');

0 comments on commit df54863

Please sign in to comment.