r718: run-fermi.pl supports BCR

traymond · Jul 11, 2012 · 7e61e5a · 7e61e5a
1 parent 3b37fcf
commit 7e61e5a
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 9 deletions.
diff --git a/fermi.h b/fermi.h
@@ -4,7 +4,7 @@
 #include <stdint.h>
 #include <stdlib.h>
 
-#define FERMI_VERSION "1.0-r716"
+#define FERMI_VERSION "1.0-r718"
 
 #define FM_MASK30 0x3fffffff
 

diff --git a/ropebwt.c b/ropebwt.c
@@ -72,7 +72,7 @@ int main_ropebwt(int argc, char *argv[])
 		bcr = bcr_init(flag&FLAG_THR, tmpfn);
 		if (!(flag&FLAG_NON)) fprintf(stderr, "Warning: With bcr, an ambiguous base will be converted to a random base\n");
 	} else if (algo == BPR) bpr = bpr_init(max_nodes, max_runs);
-	fp = gzopen(argv[optind], "rb");
+	fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "rb") : gzdopen(fileno(stdin), "rb");
 	ks = kseq_init(fp);
 	while (kseq_read(ks) >= 0) {
 		int l = ks->seq.l;

diff --git a/run-fermi.pl b/run-fermi.pl
@@ -8,7 +8,7 @@
 
 sub main {
 	my %opts = (e=>'fermi', t=>2, p=>'fmdef', k=>50);
-	getopts('e:t:p:Pck:Dl:', \%opts);
+	getopts('e:t:p:Pck:Dbl:', \%opts);
 	$opts{P} = 1 if defined($opts{c});
 	$opts{l} = defined($opts{l})? "-l$opts{l}" : "";
 	my $min_clean_o = int($opts{k} * 1.2 + .499);
@@ -18,6 +18,7 @@ sub main {
 Options: -P        the input files are paired (ends in separate files)
          -c        the input is collated\/initerleaved FASTQ (two ends in the same file)
          -D        halve the number of jobs for building the split index
+         -b        use bcr for constructing FM-index
          -e FILE   fermi executable [$opts{e}]
          -t INT    number of threads [$opts{t}]
          -p STR    prefix of output files [$opts{p}]
@@ -50,19 +51,29 @@ sub main {
 
 	push(@lines, "# Construct the FM-index for raw sequences");
 	my $pre = "$opts{p}.raw";
-	push(@lines, "$pre.split.log:$in_list");
-	push(@lines, "\t$fqs | \$(FERMI) splitfa - $pre $n_split 2> $pre.split.log\n");
-	&build_fmd(\@lines, $n_split, $pre, $opts{t}); # do not trim for the initial index
+	if (defined($opts{b})) {
+		push(@lines, "$pre.fmd:$in_list");
+		push(@lines, "\t$fqs | \$(FERMI) ropebwt -a bcr -btNf $pre.tmp - > \$@ 2> \$@.log", '');
+	} else {
+		push(@lines, "$pre.split.log:$in_list");
+		push(@lines, "\t$fqs | \$(FERMI) splitfa - $pre $n_split 2> $pre.split.log\n");
+		&build_fmd(\@lines, $n_split, $pre, $opts{t}); # do not trim for the initial index
+	}
 
 	push(@lines, "# Error correction");
 	push(@lines, "$opts{p}.ec.fq.gz:$opts{p}.raw.fmd");
 	push(@lines, "\t$fqs | \$(FERMI) correct -".(defined($opts{P})? 'p' : '')."t $opts{t} $opts{l} \$< - 2> \$@.log | gzip -1 > \$@\n");
 
 	push(@lines, "# Construct the FM-index for corrected sequences");
 	$pre = "$opts{p}.ec";
-	push(@lines, "$pre.split.log:$opts{p}.ec.fq.gz");
-	push(@lines, "\t\$(FERMI) fltuniq \$< 2> $opts{p}.fltuniq.log | \$(FERMI) splitfa - $pre $n_split 2> \$@\n");
-	&build_fmd(\@lines, $n_split, $pre, $opts{t});
+	if (defined($opts{b})) {
+		push(@lines, "$pre.fmd:$opts{p}.ec.fq.gz");
+		push(@lines, "\t\$(FERMI) fltuniq \$< 2> $opts{p}.fltuniq.log | \$(FERMI) ropebwt -a bcr -btf $pre.tmp - > \$@ 2> \$@.log", '');
+	} else {
+		push(@lines, "$pre.split.log:$opts{p}.ec.fq.gz");
+		push(@lines, "\t\$(FERMI) fltuniq \$< 2> $opts{p}.fltuniq.log | \$(FERMI) splitfa - $pre $n_split 2> \$@\n");
+		&build_fmd(\@lines, $n_split, $pre, $opts{t});
+	}
 
 	push(@lines, "# Generate unitigs");
 	if (defined($opts{P})) {