Enable detailed inspection of mapping.

stela2502 · Nov 14, 2023 · 0a34044 · 0a34044
1 parent 272d4fb
commit 0a34044
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 23 deletions.
diff --git a/idenitfy_index_failing_reads.sh b/idenitfy_index_failing_reads.sh
@@ -1,10 +1,18 @@
 target/release/quantify_rhapsody -r testData/cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f testData/cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o testData/BD_results/Rustody_S1 -s mouse  -e testData/2276_20220531_chang_to_rpl36a_amplicons.fasta -a testData/MyAbSeqPanel.fasta -m 200 -v v2.96 > testData/BD_results/Rustody_S1/mapped_reads.txt
-target/release/quantify_rhapsody -r testData/cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f testData/cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o testData/BD_results/Rustody_S1 -s mouse  -e testData/2276_20220531_chang_to_rpl36a_amplicons.fasta -a testData/MyAbSeqPanel.fasta -m 200 -v v2.96 > testData/BD_results/Rustody_S1/mapped_reads.txt
+target/release/quantify_rhapsody -r testData/cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f testData/cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o testData/BD_results/Rustody_S1 -s mouse  -e testData/addOn.fa -i testData/mapperTest/index -a testData/MyAbSeqPanel.fasta -m 200 -v v2.96 > testData/BD_results/Rustody_S1_index/mapped_reads.txt
 
 Rscript idenitfy_index_failing_reads.R
 
 zgrep -f R1_ids_not_in_index.txt testData/cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -A3 | sed '/--/d' | gzip > testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz
 zgrep -f R2_ids_not_in_index.txt testData/cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -A3 | sed '/--/d' | gzip > testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz
 
 
-target/release/quantify_rhapsody_multi -r testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o testData/BD_results/Rustody_S1 -s mouse  -e testData/2276_20220531_chang_to_rpl36a_amplicons.fasta -a testData/MyAbSeqPanel.fasta -m 200 -v v2.96 
+target/release/quantify_rhapsody_multi -r testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o testData/BD_results/Rustody_S1_failed_in_indexing -s mouse  -e testData/2276_20220531_chang_to_rpl36a_amplicons.fasta -a testData/MyAbSeqPanel.fasta -m 200 -v v2.96 
+target/release/quantify_rhapsody_multi -r testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o testData/BD_results/Rustody_S1_index_failed -s mouse  -e testData/addOn.fa -i testData/mapperTest/index -a testData/MyAbSeqPanel.fasta -m 200 -v v2.96
+
+# there are still indices that can not be mapped using the muti program?!
+
+target/release/quantify_rhapsody -r testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R1_001.fastq.gz -f testData/failing_in_index_cells.1.Rhapsody_SV_index1_S1_R2_001.fastq.gz -o testData/BD_results/Rustody_S1_failed_in_indexing -s mouse  -e testData/2276_20220531_chang_to_rpl36a_amplicons.fasta -a testData/MyAbSeqPanel.fasta -m 200 -v v2.96  > testData/BD_results/Rustody_S1_failed_in_indexing/look_into_these_indexes.txt
+
+
+# damn these are all mapping using the not multi program?! How so??
diff --git a/this/src/fast_mapper.rs b/this/src/fast_mapper.rs
@@ -356,21 +356,21 @@ impl FastMapper{
     fn get_best_gene( &self, genes:&HashMap::<(usize, usize), usize>, ret: &mut Vec::<usize> ) -> bool{
         ret.clear();
 
-        // let mut report = false;
-        // if genes.len() > 2{
-        //     report = true;
-        //     eprintln!("Lots of genes matched here?: {genes:?}");
-        //     for  (gene_id, _level) in genes.keys(){
-        //         eprint!(" {}", self.names_store[*gene_id] );
-        //     }
-        //     eprint!("\n");
-        // }
+        let mut report = false;
+        if genes.len() > 2{
+            report = true;
+            eprintln!("Lots of genes matched here?: {genes:?}");
+            for  (gene_id, _level) in genes.keys(){
+                eprint!(" {}", self.names_store[*gene_id] );
+            }
+            eprint!("\n");
+        }
         if genes.len() == 1 {
             if let Some((key, _)) = genes.iter().next() {
                ret.push(key.0.clone());
-               // if report {
-               //  eprintln!("1 This was selected as good: {} or {}\n",key.0, self.names_store[key.0] )
-               // }
+               if report {
+                eprintln!("1 This was selected as good: {} or {}\n",key.0, self.names_store[key.0] )
+               }
                return true
             }
         }
@@ -399,24 +399,24 @@ impl FastMapper{
             }
             if good.len() ==1 {
                 ret.push( good[0] );
-                // if report {
-                //     eprintln!("2 This was selected as good: {} or {}\n",good[0], self.names_store[good[0]] )
-                // }
+                if report {
+                    eprintln!("2 This was selected as good: {} or {}\n",good[0], self.names_store[good[0]] )
+                }
                 return true
             }
             if genes_with_best_matches == 1{
                 ret.push( best_gene );
-                // if report {
-                //     eprintln!("3 This was selected as good: {} or {}\n",best_gene, self.names_store[best_gene] )
-                // }
+                if report {
+                    eprintln!("3 This was selected as good: {} or {}\n",best_gene, self.names_store[best_gene] )
+                }
                 return true
             }
 
 
         }
-        // if report {
-        //     eprintln!("! No good gene identified!\n" )
-        // }
+        if report {
+            eprintln!("! No good gene identified!\n" )
+        }
         return false
 
     }