Merge bug fixes and cram compatibility

Author: Yuhe Cheng

src/DupCaller.egg-info/SOURCES.txt

@@ -20,4 +20,4 @@ src/subcommands/funcs/depth.py
 src/subcommands/funcs/indels.py
 src/subcommands/funcs/learn.py
 src/subcommands/funcs/misc.py
-src/subcommands/funcs/prob.py
+src/subcommands/funcs/prob.py

src/DupCaller.egg-info/dependency_links.txt

@@ -0,0 +1 @@
+

src/subcommands/Caller.py

@@ -37,12 +37,12 @@ def do_call(args):
         "threads": args.threads,
         # "amperr": args.amperrs,
         # "amperri": args.amperri,
-        "amperr_file": args.output + "/" + args.output + ".amp.tn.txt",
-        "amperri_file": args.output + "/" + args.output + ".amp.id.txt",
+        "amperr_file": args.output + ".amp.tn.txt",
+        "amperri_file": args.output + ".amp.id.txt",
         # "dmgerr": args.dmgerrs,
         # "dmgerri": args.dmgerri,
-        "dmgerr_file": args.output + "/" + args.output + ".dmg.tn.txt",
-        "dmgerri_file": args.output + "/" + args.output + ".dmg.id.txt",
+        "dmgerr_file": args.output + ".dmg.tn.txt",
+        "dmgerri_file": args.output + ".dmg.id.txt",
         "mutRate": args.mutRate,
         "pcutoff": args.threshold,
         "mapq": args.mapq,
@@ -133,7 +133,7 @@ def do_call(args):
         except OSError as e:
             if e.errno != errno.EEXIST:
                 raise
-    bamObject = BAM(args.bam, "rb")
+    bamObject = BAM(args.bam, "rb", args.reference)
 
     """
     Execulte variant calling
@@ -184,7 +184,7 @@ def do_call(args):
             # print(args.threads)
             # if args.normalBam:
             cutSites, chunkSize, contigs = splitBamRegions(
-                [args.bam], args.threads, contigs, args.windowSize
+                [args.bam], args.threads, contigs, args.windowSize, args.reference
             )
             # else:
             # cutSites, chunkSize, contigs = splitBamRegions(
@@ -284,17 +284,17 @@ def do_call(args):
             mismatch_dmg_profile, columns=["A", "T", "C", "G"], index=num2trinuc
         )
         # np.savetxt(params["output"] + "/" + args.output + ".amp.tn.txt",np.hstack([trinuc_cols[0:32],mismatch_profile]),delimiter="\t",header=" \tA\tT\tC\tG\n")
-        amp_tn_pd.to_csv(params["output"] + "/" + args.output + ".amp.tn.txt", sep="\t")
+        amp_tn_pd.to_csv(args.output + ".amp.tn.txt", sep="\t")
         np.savetxt(
-            params["output"] + "/" + args.output + ".amp.id.txt",
+            args.output + ".amp.id.txt",
             indelerr_profile,
             delimiter="\t",
             fmt="%d",
         )
-        dmg_tn_pd.to_csv(params["output"] + "/" + args.output + ".dmg.tn.txt", sep="\t")
+        dmg_tn_pd.to_csv(args.output + ".dmg.tn.txt", sep="\t")
         # np.savetxt(params["output"] + "/" + args.output + ".dmg.tn.txt",np.hstack([trinuc_cols,mismatch_dmg_profile]),delimiter="\t",header=" \tA\tT\tC\tG\n")
         np.savetxt(
-            params["output"] + "/" + args.output + ".dmg.id.txt",
+            args.output + ".dmg.id.txt",
             indelerr_dmg_profile,
             delimiter="\t",
             fmt="%d",
@@ -368,11 +368,11 @@ def do_call(args):
             print("....Splitting genomic regions for parallel execution.....")
             if args.normalBams:
                 cutSites, chunkSize, contigs = splitBamRegions(
-                    [args.bam], args.threads, contigs, args.windowSize
+                    [args.bam], args.threads, contigs, args.windowSize, args.reference
                 )
             else:
                 cutSites, chunkSize, contigs = splitBamRegions(
-                    [args.bam], args.threads, contigs, args.windowSize
+                    [args.bam], args.threads, contigs, args.windowSize, args.reference
                 )
             currentContigIndex = 0
             usedTime = (time.time() - startTime) / 60
@@ -491,7 +491,7 @@ def do_call(args):
         FPAll = sum(FPs, [])
         RPAll = sum(RPs, [])
 
-    tBam = BAM(args.bam, "rb")
+    tBam = BAM(args.bam, "rb", args.reference)
     contigs = tBam.references
     # print(contigs)
     chromDict = {contig: tBam.get_reference_length(contig) for contig in contigs}
@@ -518,11 +518,11 @@ def do_call(args):
     }
     filterDict = {"PASS": "All filter Passed"}
     vcfLines = createVcfStrings(chromDict, infoDict, formatDict, filterDict, mutsAll)
-    with open(args.output + "/" + args.output + "_snv.vcf", "w") as vcf:
+    with open(args.output + "_snv.vcf", "w") as vcf:
         vcf.write(vcfLines)
 
     vcfLines = createVcfStrings(chromDict, infoDict, formatDict, filterDict, indelsAll)
-    with open(args.output + "/" + args.output + "_indel.vcf", "w") as vcf:
+    with open(args.output + "_indel.vcf", "w") as vcf:
         vcf.write(vcfLines)
 
     burden_naive = muts_num / (coverage)
@@ -531,7 +531,7 @@ def do_call(args):
     pass_duprate = unique_read_num / pass_read_num
 
     with open(
-        params["output"] + "/" + args.output + "_duplex_group_stats.txt", "w"
+        args.output + "_duplex_group_stats.txt", "w"
     ) as f:
         f.write(
             "duplex_group_strand_composition\tduplex_group_number\t\
@@ -570,13 +570,13 @@ def do_call(args):
     trinuc_by_duplex_group = pd.DataFrame(duplex_read_num_trinuc)
     trinuc_by_duplex_group.insert(0, "", trinuc_list)
     trinuc_by_duplex_group.to_csv(
-        params["output"] + "/" + args.output + "_trinuc_by_duplex_group.txt",
+        args.output + "_trinuc_by_duplex_group.txt",
         sep="\t",
         index=False,
     )
 
     muts_by_group = np.loadtxt(
-        params["output"] + "/" + args.output + "_duplex_group_stats.txt",
+        args.output + "_duplex_group_stats.txt",
         skiprows=1,
         dtype=float,
         delimiter="\t",
@@ -610,7 +610,7 @@ def do_call(args):
     lgd2 = mpatches.Patch(color="blue", label="Naive")
     plt.legend(handles=[lgd1, lgd2])
     plt.savefig(
-        params["output"] + "/" + args.output + "_burden_by_duplex_group_size.png"
+        args.output + "_burden_by_duplex_group_size.png"
     )
     """
     if len(FPAll + RPAll) != 0:
@@ -633,7 +633,7 @@ def do_call(args):
         clonal_num = 0
     """
 
-    with open(params["output"] + "/" + args.output + "_stats.txt", "w") as f:
+    with open(args.output + "_stats.txt", "w") as f:
         f.write(f"Number of Read Families\t{unique_read_num}\n")
         f.write(f"Number of Pass-filter Reads\t{pass_read_num}\n")
         f.write(f"Number of Effective Read Families\t{duplex_num}\n")

src/subcommands/Learn.py

@@ -68,7 +68,7 @@ def do_learn(args):
         except OSError as e:
             if e.errno != errno.EEXIST:
                 raise
-    bamObject = BAM(args.bam, "rb")
+    bamObject = BAM(args.bam, args.reference, "rb")
 
     """
     Execulte variant calling

src/subcommands/funcs/call.py

@@ -11,6 +11,7 @@
 from pysam import TabixFile as BED
 from pysam import VariantFile as VCF
 import pysam
+import os
 import re
 import h5py
 
@@ -207,8 +208,8 @@ def nums2str(nums, num2base="ATCG"):
     return "".join(bases)
 
 
-def bamIterateMultipleRegion(bam, regions):
-    bamObject = BAM(bam, "rb")
+def bamIterateMultipleRegion(bam, regions, ref):
+    bamObject = BAM(bam, "rb", ref)
     for region in regions:
         for rec in bamObject.fetch(*region):
             if len(region) >= 2:
@@ -223,7 +224,7 @@ def callBam(params, processNo):
     nbams = params["normalBams"]
     regions = params["regions"]
     if params["germline"]:
-        germline = VCF(params["germline"], "rb")
+        germline = VCF(params["germline"], mode = "rb")
     else:
         germline = None
     all_chroms = [_[0] for _ in regions]
@@ -235,7 +236,9 @@ def callBam(params, processNo):
     pcut = params["pcutoff"]
     isLearn = params.get("isLearn", False)
     nn = processNo
-    output = "tmp/" + params["output"] + "_" + str(nn)
+    output = os.path.join("tmp", params["output"].lstrip('/') + "_" + str(nn))
+    if not os.path.exists(os.path.dirname(output)):
+        os.makedirs(os.path.dirname(output))
     if params["noise"]:
         noise = BED(params["noise"])
     else:
@@ -252,6 +255,7 @@ def callBam(params, processNo):
         feature_beds = [BED(_) for _ in params["feature_files"]]
     else:
         feature_beds = None
+
     base2num = {"A": 0, "T": 1, "C": 2, "G": 3}
     num2base = "ATCG"
     muts = []
@@ -302,7 +306,7 @@ def callBam(params, processNo):
     params["trinuc_convert"] = trinuc_convert_np
     params["trinuc2num_dict"] = trinuc2num
     params["num2trinuc_list"] = num2trinuc
-
+    print(params["amperr_file"])
     ### Load amp error matrix
     if not params["amperr_file"]:
         prob_amp = params["amperr"]
@@ -417,16 +421,16 @@ def callBam(params, processNo):
     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)
-
     # Initialize
+
     total_coverage = 0
     total_coverage_indel = 0
     starttime = time.time()
-    tumorBam = BAM(bam, "rb")
+    tumorBam = BAM(bam, "rb", params.get("reference"))
     if nbams:
         normalBams = list()
         for nbam in nbams:
-            normalBams.append(BAM(nbam, "rb"))
+            normalBams.append(BAM(nbam, "rb", params.get("reference")))
     else:
         normalBams = None
     currentStart = -1
@@ -457,7 +461,7 @@ def callBam(params, processNo):
         )
         read_blacklist = set()
         rec_num = 0
-        for rec, region in bamIterateMultipleRegion(bam, regions):
+        for rec, region in bamIterateMultipleRegion(bam, regions, params.get("reference")):
             rec_num += 1
             if rec.query_name in read_blacklist or rec.is_unmapped:
                 continue
@@ -477,7 +481,7 @@ def callBam(params, processNo):
         print(
             f"Process {str(processNo)}: finished screening highly damaged reads in {currentTime: .2f} minutes. Blacklisted {len(read_blacklist)} ({percent_blocked: .2f}%) possible highly damaged read and started variant calling."
         )
-    for rec, region in bamIterateMultipleRegion(bam, regions):
+    for rec, region in bamIterateMultipleRegion(bam, regions, params.get("reference")):
         recCount += 1
         if recCount == currentCheckPoint:
             currentTime = (time.time() - starttime) / 60

src/subcommands/funcs/misc.py

@@ -48,13 +48,13 @@ def createVcfStrings(chromDict, infoDict, formatDict, filterDict, recs):
     return "\n".join(lines) + "\n"
 
 
-def bamReadCount(bam, chrom, start, end):
-    count = getAlignmentObject(bam, "rb").count(chrom, start, end)
+def bamReadCount(bam, chrom, start, end, ref):
+    count = getAlignmentObject(bam, "rb", ref).count(chrom, start, end)
     return count
 
 
-def splitBamRegions(bams, num, contigs, step):
-    bamObject = getAlignmentObject(bams[0], "rb")
+def splitBamRegions(bams, num, contigs, step, ref):
+    bamObject = getAlignmentObject(bams[0], "rb", ref)
     contigs_set = set(contigs)
     contigs_sorted = [_ for _ in bamObject.references if _ in contigs]
 
@@ -88,6 +88,7 @@ def splitBamRegions(bams, num, contigs, step):
                         contig,
                         current_start + k * step,
                         current_start + k * step + step,
+                        ref,
                     )
                     for k in range(num)
                 ]
@@ -97,7 +98,7 @@ def splitBamRegions(bams, num, contigs, step):
                 current_start += num * step
                 current_window += num
             current_arguments = [
-                (bam, contig, current_start + k * step, current_start + k * step + step)
+                (bam, contig, current_start + k * step, current_start + k * step + step, ref)
                 for k in range(window_nums_cumulative[c] - current_window)
             ]
             pool = Pool()
@@ -126,11 +127,11 @@ def splitBamRegions(bams, num, contigs, step):
 
 
 ## Allow reading either bam or cram as bamObject
-def getAlignmentObject(bam, refpath=None):
+def getAlignmentObject(bam, mode, refpath=None):
     if bam.endswith(".bam"):
-        bamObject = BAM(bam, "rb")
+        bamObject = BAM(bam, mode)
     elif bam.endswith(".cram"):
-        bamObject = BAM(bam, "rc", reference_filename=refpath)
+        bamObject = BAM(bam, mode, reference_filename=refpath)
     else:
         raise NameError(f"{bam} should have .bam or .cram as file extension")
     return bamObject