fix bug

Author: Yuhe Cheng

README.md

@@ -16,7 +16,7 @@ The complete DupCaller pipeline also requires the following tools for data prepr
 The tool uses pip for installing scripts and prerequisites. To install DupCaller, simply clone this repository and install via pip:
 
 ```bash
-git clone https://github.com/AlexandrovLab/DupCaller.git
+git clone https://github.com/YuheCheng62/DupCaller.git
 cd DupCaller
 pip install .
 ```

hs_err_pid74663.log

@@ -418,6 +418,20 @@
         help="contigs to consider for trinucleotide calculation",
         default=["chr" + str(_) for _ in range(1, 23, 1)] + ["chrX"],
     )
+    estimate_parser.add_argument(
+        "-c",
+        "--clonal",
+        type=bool,
+        help="If True, mutations detected in more than one molecule will be considered as clonal mutations",
+        default=False,
+    )
+    estimate_parser.add_argument(
+        "-d",
+        "--dilute",
+        type=bool,
+        help="Set to true when sample and matched normal are from the same starting DNA material",
+        default=False,
+    )
 
     summarize_parser = subparsers.add_parser(
         "summarize", help="Summarize results from multiple samples"

src/DupCaller.py

@@ -6,7 +6,7 @@
 import pandas as pd
 from matplotlib import pyplot as plt
 from pysam import VariantFile as VCF
-from scipy.stats import chi2
+from scipy.stats import chi2, barnard_exact
 
 
 def calculate_ref_trinuc(args):
@@ -168,10 +168,24 @@ def do_estimate(args):
     revcomp = {"A": "T", "T": "A", "C": "G", "G": "C"}
     for nn, duplex_no in enumerate(trinuc_by_rf.columns):
         duplex_no_dict[duplex_no] = nn
+
+    if args.dilute:
+        vcf_out = VCF(args.prefix + "/" + sample + "_flt.vcf", "w", header=vcf.header)
     for rec in vcf.fetch():
         F1R2 = rec.info["F1R2"]
         F2R1 = rec.info["F2R1"]
-        duplex_no = str(min(F1R2, F2R1)) + "+" + str(max(F1R2, F2R1))
+        TAC = rec.samples["TUMOR"]["AC"]
+        if TAC > 1 and args.dilute:
+            TDP = rec.samples["TUMOR"]["DP"]
+            NAC = rec.samples["NORMAL"]["AC"]
+            NDP = rec.samples["NORMAL"]["DP"]
+            barnard_p = barnard_exact([[TAC, TDP - TAC], [NAC, NDP - NAC]]).pvalue
+            if barnard_p <= 0.05:
+                continue
+        if args.dilute:
+            vcf_out.write(rec)
+        # duplex_no = str(min(F1R2, F2R1)) + "+" + str(max(F1R2, F2R1))
+        duplex_no = str(F1R2) + "+" + str(F2R1)
         ref = rec.ref
         if ref == "C" or ref == "T":
             trinuc = rec.info["TN"]

src/subcommands/Estimate.py

@@ -367,6 +367,7 @@ def callBam(params, processNo):
     dmgmat = dmgmat / dmgmat.sum(axis=1, keepdims=True)
     dmgmat_min_error = dmgmat.min(axis=1, keepdims=True)
     dmgmat = np.concatenate([dmgmat, dmgmat_min_error], axis=1)
+    # dmgmat += 1e-9
 
     params["dmgmat_top"] = dmgmat
     params["trinuc2num_dict"] = trinuc2num
@@ -380,6 +381,7 @@ def callBam(params, processNo):
             ]
     dmgmat_rev_min_error = dmgmat_rev.min(axis=1, keepdims=True)
     dmgmat_rev = np.concatenate([dmgmat_rev, dmgmat_rev_min_error], axis=1)
+    # dmgmat_rev += 1e-9
     params["dmgmat_rev_top"] = dmgmat_rev
 
     dmgmat_b = np.vstack((dmgmat[32:64, [1, 0, 3, 2]], dmgmat[:32, [1, 0, 3, 2]]))
@@ -390,6 +392,8 @@ def callBam(params, processNo):
     dmgmat_rev_b_min_error = dmgmat_rev_b.min(axis=1, keepdims=True)
     dmgmat_b = np.concatenate([dmgmat_b, dmgmat_b_min_error], axis=1)
     dmgmat_rev_b = np.concatenate([dmgmat_rev_b, dmgmat_rev_b_min_error], axis=1)
+    # dmgmat_b += 1e-9
+    # dmgmat_rev_b += 1e-9
     params["dmgmat_bot"] = dmgmat_b
     params["dmgmat_rev_bot"] = dmgmat_rev_b
 
@@ -402,14 +406,15 @@ def callBam(params, processNo):
             current_row = dmgmat_indel[nn, :]
             current_row[current_row == 0] = dmgmat_indel[nn - 1, :][current_row == 0]
             dmgmat_indel[nn, :] = current_row
+        # dmgmat_indel += 1e-9
         params["dmgmat_indel"] = dmgmat_indel
     else:
         params["dmgmat_indel"] = np.ones([40, 11]) * params["dmgerri"]
         dmgmat_indel = params["dmgmat_indel"]
 
     dmgmat_indel_rev = np.fliplr(dmgmat_indel)
+    # dmgmat_indel_rev += 1e-9
     params["dmgmat_indel_rev"] = dmgmat_indel_rev
-
     params["dmgmat_indel_mean"] = np.mean(dmgmat_indel, axis=1)
     params["dmgmat_indel_rev_mean"] = np.mean(dmgmat_indel_rev, axis=1)
 
@@ -936,6 +941,7 @@ def callBam(params, processNo):
                     # else:
                     ### Calculate genotype probability
                     if not any(indel_bool) or isLearn:
+                        # if 1:
                         if isLearn:
                             (
                                 mismatch_now,
@@ -1429,6 +1435,7 @@ def callBam(params, processNo):
                     indel_dict[indel_str] = 1
             # else:
             ### Calculate genotype probability
+            # if 1:
             if not any(indel_bool) or isLearn:  # or len(indels_pass) == 0:
                 if isLearn:
                     (