Minor fix to prevent rflen from going beyond the true reference lengt…

…h in alignWithALTs_recur function
infphilo · Jun 7, 2017 · 1038abd · 1038abd
1 parent ff07bf6
commit 1038abd
Show file tree

Hide file tree

Showing 4 changed files with 16 additions and 8 deletions.
diff --git a/evaluation/tests/HLA_novel/README b/evaluation/tests/HLA_novel/README
@@ -2,8 +2,8 @@
 
    for i in `ls -l CP_80/LP*1.fq.gz | awk '{print $NF}' | cut -d'/' -f2 | cut -d'.' -f1`; do echo $i; ./hisatgenotype_locus.py --locus-list A --num-editdist 2 --assembly --reads CP_80/$i.extracted.1.fq.gz,CP_80/$i.extracted.2.fq.gz --assembly-base assembly_graph_$i 2> $i.tmp; done
 
-   hisatgenotype_scripts/hisat2_test_HLA_genotyping_CP.py -p 3 -v --num-editdist 2 --max-sample 400 --assembly --hla-list A CP_400 --out-dir test_A > cp_400_hla.txt
-   hisatgenotype_scripts/compare_HLA_Omixon.py cp_80_hla.txt ../CAAPA/omixon_caapa_hla.txt
+   hisatgenotype_locus_samples.py -p 3 -v --max-sample 400 --region-list hla.A --read-dir CP_400 --out-dir test_A > cp_400_A.txt
+   hisatgenotype_scripts/compare_HLA_Omixon.py cp_400_A.txt ../CAAPA/omixon_caapa_hla.txt
 
    hisatgenotype_locus.py --aligner-list hisat2.graph --locus-list A --num-editdist 2 --assembly -1 CP_80/LP6005106-DNA_H02.extracted.fq.1.gz -2 CP_80/LP6005106-DNA_H02.extracted.fq.2.gz --display-alleles A*23:17,A*02:02:01:01
 

diff --git a/hi_aligner.h b/hi_aligner.h
@@ -2621,6 +2621,14 @@ index_t GenomeHit<index_t>::alignWithALTs_recur(
     if(ht_llist.size() <= dep) ht_llist.expand();
     if(raw_refbufs.size() <= dep) raw_refbufs.expand();
     if(rfoff < -16) return 0;
+    size_t contig_len = ref.approxLen(tidx);
+    if(rfoff >= contig_len) return 0;
+    if(rfoff >= 0 && rfoff + rflen > contig_len) {
+        rflen = contig_len - rfoff;
+    } else if(rfoff < 0 && rflen > contig_len) {
+        rflen = contig_len;
+    }
+    if(rflen == 0) return 0;
     if(rfseq == NULL) {
         SStringExpandable<char>& raw_refbuf = raw_refbufs[dep];
         raw_refbuf.resize(rflen + 16 + 16);

diff --git a/hisat2.xcodeproj/project.pbxproj b/hisat2.xcodeproj/project.pbxproj
@@ -874,7 +874,7 @@
 					"MACOS=1",
 					POPCNT_CAPABILITY,
 					BOWTIE2,
-					"HISAT2_VERSION=\"\\\"2.0.6\\\"\"",
+					"HISAT2_VERSION=\"\\\"2.1.0\\\"\"",
 					"HISAT2_VERSION2=\"\\\"`cat VERSION`\\\"\"",
 					"BUILD_HOST=\"\\\"`hostname`\\\"\"",
 					"BUILD_TIME=\"\\\"`date\\\"\"",
@@ -887,7 +887,7 @@
 				GCC_WARN_UNINITIALIZED_AUTOS = YES_AGGRESSIVE;
 				GCC_WARN_UNUSED_FUNCTION = YES;
 				GCC_WARN_UNUSED_VARIABLE = YES;
-				MACOSX_DEPLOYMENT_TARGET = 10.11;
+				MACOSX_DEPLOYMENT_TARGET = 10.12;
 				MTL_ENABLE_DEBUG_INFO = YES;
 				ONLY_ACTIVE_ARCH = YES;
 				SDKROOT = macosx;
@@ -927,7 +927,7 @@
 					"MACOS=1",
 					POPCNT_CAPABILITY,
 					BOWTIE2,
-					"HISAT2_VERSION=\"\\\"2.0.6\\\"\"",
+					"HISAT2_VERSION=\"\\\"2.1.0\\\"\"",
 					"BUILD_HOST=\"\\\"`hostname`\\\"\"",
 					"BUILD_TIME=\"\\\"`date`\\\"\"",
 					"COMPILER_VERSION=\"\\\"`$(CXX) -v 2>&1 | tail -1`\\\"\"",
@@ -939,7 +939,7 @@
 				GCC_WARN_UNINITIALIZED_AUTOS = YES_AGGRESSIVE;
 				GCC_WARN_UNUSED_FUNCTION = YES;
 				GCC_WARN_UNUSED_VARIABLE = YES;
-				MACOSX_DEPLOYMENT_TARGET = 10.11;
+				MACOSX_DEPLOYMENT_TARGET = 10.12;
 				MTL_ENABLE_DEBUG_INFO = NO;
 				SDKROOT = macosx;
 			};

diff --git a/hisatgenotype_scripts/compare_HLA_Omixon.py b/hisatgenotype_scripts/compare_HLA_Omixon.py
@@ -91,13 +91,13 @@ def compare(hisatgenotype_fname, omixon_fname):
                         break
 
             # DK - for debugging purposes
-            # """
+            """
             if gene in ["A", "B", "C", "DQA1", "DQB1", "DRB1"] and num_match < 2:
                 print sample
                 print "\t", omixon_gene
                 print "\t", hisat_gene
                 # sys.exit(1)
-            # """
+            """
 
             assert num_match < len(count)
             count[num_match] += 1