Log info and seed design improvement (#131)

* Update location of error message

* Correct spacing

* Use different shared mem location

* Make tmp dir

* Try newer python

* Test pinning dependency

* Move pinned version

* goldrush_path: expand reads skipped log

* goldrush_path: change variable name typo

* goldrush: change prefix

* spaced_seed: deterministic random space seed

* goldrush: only use seed preset if k and w are default

* goldrush: remove tab indent

* Update goldpolish to latest commit

* Remove debugging statements in azure-pipelines.yml and goldrush_test_demo.sh

* Update prefix in help docs

* goldrush: reverse seed preset logic regarding default k and w

* goldrush: remove s from help

* readme: update help to reflect goldrush

* Update p in goldrush help

---------

Co-authored-by: lcoombe <lauren.e.coombe@gmail.com>

Author: jwcodee

README.md

@@ -63,13 +63,12 @@ default 5 rounds)
         o                       occupancy of the miBF [0.1]
         x                       threshold for number of hits in miBF for a given frame to be considered assigned [10]
         h                       number of seed patterns to use [3]
-        s                       spaced seed design [1011011110110111101101]
         m                       minimum read length [20000]
         M                       maximum number of silver paths to generate [5]
         r                       ratio of full genome in golden path [0.9]
         P                       minimum average phred score for each read [15]
         d                       remove reads with greater or equal than d difference between average phred quality of first half and second half of the read [5]
-        p                       prefix to use for the output paths [w16_x10]
+        p                       prefix to use for the output paths [goldrush_asm]
 
         Tigmint-long options:
         span                    min number of spanning molecules [2]

azure-pipelines.yml

@@ -14,8 +14,8 @@ jobs:
     displayName: Create Anaconda environment
   - script: |
       source activate goldrush_CI
-      conda install --yes -c conda-forge mamba python=3.9
-      mamba install --yes -c conda-forge -c bioconda compilers meson gperftools sdsl-lite boost-cpp sparsehash btllib libdivsufsort minimap2 tigmint ntlink miller
+      conda install --yes -c conda-forge mamba python=3.10 
+      mamba install --yes -c conda-forge -c bioconda compilers meson gperftools sdsl-lite boost-cpp sparsehash btllib libdivsufsort minimap2 tigmint ntlink miller  
     displayName: Install dependencies
   - script: |
       source activate goldrush_CI
@@ -47,8 +47,8 @@ jobs:
     displayName: Create Anaconda environment
   - script: |
       source activate goldrush_CI
-      conda install --yes -c conda-forge mamba python=3.9
-      mamba install --yes -c conda-forge -c bioconda compilers meson gperftools sdsl-lite boost-cpp sparsehash btllib libdivsufsort minimap2 tigmint ntlink miller
+      conda install --yes -c conda-forge mamba python=3.10
+      mamba install --yes -c conda-forge -c bioconda compilers meson gperftools sdsl-lite boost-cpp sparsehash btllib libdivsufsort minimap2 tigmint ntlink miller 
     displayName: Install dependencies
   - script: |
       source activate goldrush_CI

bin/goldrush

@@ -68,7 +68,7 @@ h=3
 s=1011011110110111101101
 r=0.9
 M=5
-p=w$(w)_x$(x)
+p=goldrush_asm
 p1=$(p)_silver_path
 p2=$(p)_golden_path
 P=15
@@ -165,13 +165,12 @@ help:
 	@echo "	o			occupancy of the miBF [$(o)]"
 	@echo "	x			threshold for number of hits in miBF for a given frame to be considered assigned [$(x)]"
 	@echo "	h			number of seed patterns to use [$(h)]"
-	@echo "	s			spaced seed design [$(s)]"
 	@echo "	m			minimum read length [$(m)]"
 	@echo "	M			maximum number of silver paths to generate [$(M)]"
 	@echo "	r			ratio of full genome in golden path [$(r)]"
 	@echo "	P			minimum average phred score for each read [$(P)]"
 	@echo "	d			remove reads with greater or equal than d difference between average phred quality of first half and second half of the read [$(d)]"
-	@echo "	p			prefix to use for the output paths [w$(w)_x$(x)]"
+	@echo "	p			prefix to use for the output paths [$(p)]"
 	@echo ""
 	@echo "	Tigmint-long options:"
 	@echo "	span			min number of spanning molecules [$(span)]"
@@ -231,14 +230,26 @@ endif
 goldrush-path: $(p2).fa check-G check-reads clean
 
 $(p2).fa: $(p1)_all.fq
-	$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p2) -i $< -h $(h) -j $(t) -P $(P) -d $(d) -x$(x) -s $(s) -g $(G) -b $(b)  -m 0
+ifneq ($(k), 22)
+		$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p2) -i $< -h $(h) -j $(t) -P $(P) -d $(d) -x$(x) -g $(G) -b $(b)  -m 0
+else ifneq ($(w), 16)
+		$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p2) -i $< -h $(h) -j $(t) -P $(P) -d $(d) -x$(x) -g $(G) -b $(b)  -m 0
+else
+		$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p2) -i $< -h $(h) -j $(t) -P $(P) -d $(d) -x$(x) -s $(s) -g $(G) -b $(b)  -m 0
+endif
 	echo "Done GoldRush-Path! Golden path can be found in: $@"
 
 $(p1)_all.fq: $(p1)_$(M).fq
 	cat $(p1)_*.fq > $@
 
 $(p1)_$(M).fq: $(long_reads)
-	$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p1) -i $<  -h $(h) -j $(t) -x$(x) -P $(P) -d $(d) -s $(s) -g $(G) -b $(b) -r $(r) --silver_path -M $(M) -m $(m)
+ifneq ($(k), 22)
+		$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p1) -i $<  -h $(h) -j $(t) -x$(x) -P $(P) -d $(d) -g $(G) -b $(b) -r $(r) --silver_path -M $(M) -m $(m)
+else ifneq ($(w), 16)
+		$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p1) -i $<  -h $(h) -j $(t) -x$(x) -P $(P) -d $(d) -g $(G) -b $(b) -r $(r) --silver_path -M $(M) -m $(m)
+else
+		$(time) goldrush-path  -k $(k) -w $(w) -t $(tile) -u $(u) -a $(a) -o $(o) -p $(p1) -i $<  -h $(h) -j $(t) -x$(x) -P $(P) -d $(d) -s $(s) -g $(G) -b $(b) -r $(r) --silver_path -M $(M) -m $(m)
+endif
 
 %.racon-polished.fa: %.fa.$(long_reads).sam %.fa
 	$(time) racon -u -t$(t) $(long_reads) $^ > $@

goldrush_path/goldrush_path.cpp

@@ -140,13 +140,16 @@ fill_bit_vector(const std::string& input_file,
   }
   size_t num_reads = 0;
   size_t num_passed_reads = 0;
-  size_t size_num_reads_skipped_by_phred = 0;
-  size_t size_num_reads_skipped_by_delta = 0;
+  size_t num_reads_skipped_by_phred = 0;
+  size_t num_reads_skipped_by_delta = 0;
+  size_t num_reads_skipped_by_length = 0;
 #pragma omp parallel
   for (const auto record : reader) {
 #pragma omp atomic
     ++num_reads;
     if (record.seq.size() < min_seq_len) {
+#pragma omp atomic
+      ++num_reads_skipped_by_length;
       continue;
     }
     const auto phred_stat = calc_phred_average(record.qual);
@@ -155,11 +158,11 @@ fill_bit_vector(const std::string& input_file,
       if (opt::verbose) {
         if (phred_stat.first < opt::phred_min) {
 #pragma omp atomic
-          ++size_num_reads_skipped_by_phred;
+          ++num_reads_skipped_by_phred;
         }
         if (phred_stat.second >= opt::phred_delta) {
 #pragma omp atomic
-          ++size_num_reads_skipped_by_delta;
+          ++num_reads_skipped_by_delta;
         }
       }
 #pragma omp critical
@@ -183,9 +186,15 @@ fill_bit_vector(const std::string& input_file,
               << floor((double)(num_reads - num_passed_reads) / num_reads)
               << "\n"
               << "num_reads_skipped_by_phred: "
-              << size_num_reads_skipped_by_phred << "\n"
+              << num_reads_skipped_by_phred << "\n"
               << "num_reads_skipped_by_delta: "
-              << size_num_reads_skipped_by_delta << "\n"
+              << num_reads_skipped_by_delta << "\n"
+              << "num_reads_skipped_by_length: "
+              << num_reads_skipped_by_length << "\n"
+              << "Total reads skipped: "
+              << num_reads_skipped_by_phred +
+                   num_reads_skipped_by_delta +
+                   num_reads_skipped_by_length
               << std::endl;
   }
 
@@ -1083,6 +1092,14 @@ main(int argc, char** argv)
                    opt::min_length,
                    filter_out_reads);
     }
+
+  }
+  if (opt::silver_path && opt::max_paths > curr_path) {
+    std::cerr << "WARNING: Expected " << std::to_string(opt::max_paths) 
+              << " silver paths, but only " << std::to_string(curr_path)
+              <<" generated.\n" << "Possible reasons include:\n"
+              << "\t- Input reads sorted by chromosome/position\n"
+              << "\t- Genome size set too large\n";
   }
 
   std::cerr << "assigned" << std::endl;

goldrush_path/spaced_seeds.cpp

@@ -16,7 +16,7 @@ make_seed_pattern(const std::string& seed_preset,
   std::string right_seed_str;
 
   if (seed_preset == "") {
-    srand(time(NULL));
+    srand(123);
     // seed generation
     std::cerr << "Designing base symmetrical spaced seed"
               << "\n"

subprojects/goldpolish

@@ -7,7 +7,7 @@ curl -L --output test_reads.fq https://www.bcgsc.ca/downloads/btl/goldrush/test/
 
 # Run this demo to test your GoldRush installation
 echo "Launching GoldRush"
-goldrush run reads=test_reads G=1e6 t=4 p=goldrush_test -B
+goldrush run reads=test_reads G=1e6 t=4 p=goldrush_test  -B
 
 l50=$(abyss-fac goldrush_test_golden_path.goldpolish-polished.span2.dist500.tigmint.fa.k40.w250.ntLink-5rounds.fa |awk '{print $3}' |tail -n1)