carrier GQ

dellytools · Oct 10, 2023 · a3b15a6 · a3b15a6
1 parent 4436481
commit a3b15a6
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 12 deletions.
diff --git a/src/delly.cpp b/src/delly.cpp
@@ -43,19 +43,19 @@ displayUsage() {
   std::cerr << "Long-read SV calling:" << std::endl;
   std::cerr << "    lr           long-read SV discovery" << std::endl;
   std::cerr << std::endl;
-  //std::cerr << "Pan-genome based SV calling (work-in-progress):" << std::endl;
-  //std::cerr << "    pg           pan-genome SV discovery" << std::endl;
-  //std::cerr << std::endl;
-  //std::cerr << "Assembly-based SV calling (work-in-progress):" << std::endl;
-  //std::cerr << "    asm          assembly SV site discovery" << std::endl;
-  //std::cerr << std::endl;
+  std::cerr << "Pan-genome based SV calling (work-in-progress):" << std::endl;
+  std::cerr << "    pg           pan-genome SV discovery" << std::endl;
+  std::cerr << std::endl;
+  std::cerr << "Assembly-based SV calling (work-in-progress):" << std::endl;
+  std::cerr << "    asm          assembly SV site discovery" << std::endl;
+  std::cerr << std::endl;
   std::cerr << "Copy-number variant calling:" << std::endl;
   std::cerr << "    cnv          discover and genotype copy-number variants" << std::endl;
   std::cerr << "    classify     classify somatic or germline copy-number variants" << std::endl;
-  //std::cerr << std::endl;
-  //std::cerr << "Multi-sample VCF operations:" << std::endl;
-  //std::cerr << "    markdup      mark duplicate SV sites based on SV allele and GT concordance" << std::endl;
-  //std::cerr << "    compvcf      compare multi-sample VCF file to a ground truth VCF file" << std::endl;
+  std::cerr << std::endl;
+  std::cerr << "Multi-sample VCF operations:" << std::endl;
+  std::cerr << "    markdup      mark duplicate SV sites based on SV allele and GT concordance" << std::endl;
+  std::cerr << "    compvcf      compare multi-sample VCF file to a ground truth VCF file" << std::endl;
   //std::cerr << "Deprecated:" << std::endl;
   //std::cerr << "    dpe          double paired-end signatures" << std::endl;
   //std::cerr << "    chimera      ONT chimera flagging" << std::endl;

diff --git a/src/markdup.h b/src/markdup.h
@@ -380,15 +380,27 @@ namespace torali
 
 	// Find better SV
 	double gqsum1 = 0;
+	uint32_t gqn1 = 0;
 	double gqsum2 = 0;
+	uint32_t gqn2 = 0;
 	for(uint32_t k = 0; k < allsv[i].gq.size(); ++k) {
-	  gqsum1 += allsv[i].gq[k];
-	  gqsum2 += allsv[j].gq[k];
+	  if (allsv[i].gt[k] != 0) {
+	    gqsum1 += allsv[i].gq[k];
+	    ++gqn1;
+	  }
+	  if (allsv[j].gt[k] != 0) {
+	    gqsum2 += allsv[j].gq[k];
+	    ++gqn2;
+	  }
 	}
+	// Normalize
 	if ((gqsum1 == 0) && (gqsum2 == 0)) {
 	  // Use QUAL
 	  gqsum1 = allsv[i].qual;
 	  gqsum2 = allsv[j].qual;
+	} else {
+	  gqsum1 /= (double) gqn1;
+	  gqsum2 /= (double) gqn2;
 	}
 	double sharedperc = _sharedCarriers(allsv[i].gt, allsv[j].gt);
 	if (sharedperc < c.sharedcarrier) continue;