Merge from master.

GoodstudyChina · Jul 15, 2019 · f64beda · f64beda
2 parents e758d6e + 27bda6e
commit f64beda
Show file tree

Hide file tree

Showing 5 changed files with 19 additions and 7 deletions.
diff --git a/documentation/source/parameter-reference.rst b/documentation/source/parameter-reference.rst
@@ -574,6 +574,12 @@ trimReadsCoverage <integer=1>
   Minimum depth of evidence to retain bases.
 
 
+Trio binning Configuration
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. _hapUnknownFraction:
+hapUnknownFraction <float=0.05>
+   Fraction of unclassified bases to ignore for haplotype assemblies. If there are more than this fraction of unclassified bases, they are included in both haplotype assemblies.
 
 .. _grid-engine:
 

diff --git a/documentation/source/quick-start.rst b/documentation/source/quick-start.rst
@@ -184,7 +184,9 @@ The run will first bin the reads into the haplotypes (``ecoliTrio/haplotype/hapl
   --     4114 reads      6520294 bases written to haplotype file ./haplotype-unknown.fasta.gz.
 
 
-Next, the haplotypes are assembled in ``ecoliTrio/asm-haplotypeK12/asm-haplotypeK12.contigs.fasta`` and ``ecoliTrio/asm-haplotypeO157/asm-haplotypeO157.contigs.fasta``. As comparison, you can try co-assembling the datasets instead::
+Next, the haplotypes are assembled in ``ecoliTrio/asm-haplotypeK12/asm-haplotypeK12.contigs.fasta`` and ``ecoliTrio/asm-haplotypeO157/asm-haplotypeO157.contigs.fasta``. By default, if the unassigned bases are > 5% of the total, they are included in both haplotypes. This can be controlled with the :ref:`hapUnknownFraction <hapUnknownFraction>` option. 
+
+As comparison, you can try co-assembling the datasets instead::
 
  canu \
   -p asm -d ecoliHap \

diff --git a/src/pipelines/canu.pl b/src/pipelines/canu.pl
@@ -730,11 +730,15 @@ ($$)
 
     #  Ignore the unknown reads if there aren't that many.
 
-    my $withUnknown = (($totBases > 0) && ($hapBases{"unknown"} / $totBases < 0.02)) ? 0 : 1;
+    my $unknownFraction =  getGlobal("hapUnknownFraction");
+    my $withUnknown = (($totBases > 0) && ($hapBases{"unknown"} / $totBases < $unknownFraction)) ? 0 : 1;
 
     if ($withUnknown == 0) {
         print STDERR "--\n";
-        print STDERR "-- Too few bases in unassigned reads to care; don't use them in assemblies.\n";
+        print STDERR "-- Fewer than " . $unknownFraction*100 . " % of bases in unassigned reads; don't use them in assemblies.\n";
+    } else {
+        print STDERR "--\n";
+        print STDERR "-- More than " .  $unknownFraction*100 . " % of bases in unassigned reads; including them in assemblies.\n";
     }
 
     #  For each haplotype, emit a script to run the assembly.

diff --git a/src/pipelines/canu/Defaults.pm b/src/pipelines/canu/Defaults.pm
@@ -1006,9 +1006,10 @@ sub setDefaults () {
 
     #####  Haplotyping
 
-    setDefault("hapMemory",      undef,  "Amount of memory, in gigabytes, to use for haplotype assignment");
-    setDefault("hapThreads",     undef,  "Number of threads to use for haplotype assignment");
-    setDefault("hapConcurrency", undef,  "Unused, there is only one process");
+    setDefault("hapUnknownFraction", 0.05,   "Fraction of allowed unknown bases before they are included in the assembly, between 0-1; default 0.05");
+    setDefault("hapMemory",          undef,  "Amount of memory, in gigabytes, to use for haplotype assignment");
+    setDefault("hapThreads",         undef,  "Number of threads to use for haplotype assignment");
+    setDefault("hapConcurrency",     undef,  "Unused, there is only one process");
 
     #####  Overlap Based Trimming
 

diff --git a/src/utility/files-buffered-implementation.H b/src/utility/files-buffered-implementation.H
@@ -112,7 +112,6 @@ readBuffer::skipAhead(char stop, bool after) {
 
     //  If we hit the fake stop, fill the buffer again and continue.
     if (_bufferPos == _bufferLen) {
-      _filePos--;
       fillBuffer();
       continue;
     }