From c17d69dae6402c8d67bb6c374b182b697ea65f04 Mon Sep 17 00:00:00 2001
From: rlyu <rlyu@svi.edu.au>
Date: Wed, 26 May 2021 16:33:53 +1000
Subject: [PATCH] add bulk sample support
---
src/sscocaller.nim | 25 ++++++++++++++-----------
1 file changed, 14 insertions(+), 11 deletions(-)
diff --git a/src/sscocaller.nim b/src/sscocaller.nim
index 9976b7a..50201b1 100755
--- a/src/sscocaller.nim
+++ b/src/sscocaller.nim
@@ -76,14 +76,19 @@ proc countAllele(rec:Variant, ibam:Bam, maxTotalReads:int,
rec_alt = rec.ALT[0][0]
for aln in ibam.query(chrom = chrom,start = rec.POS.cint-1, stop = rec.POS.cint):
var cbt = tag[string](aln, "CB")
+ var currentCB: string
if cbt.isNone:
# echo "no cb "
- continue
+ if not bulkBam:
+ continue
+ else:
+ currentCB = "bulk"
+ else:
+ currentCB = cbt.get
if aln.flag.unmapped or aln.mapping_quality.cint < mapq or aln.flag.dup: continue
## skip unmapped, duplicated, mapq low reads or aln.flag.secondary or aln.flag.supplementary
## if not aln.flag.proper_pair: continue
-
- if not barcodeTable.hasKey(cbt.get):
+ if not barcodeTable.hasKey(currentCB):
continue
var
off = aln.start.int
@@ -110,21 +115,21 @@ proc countAllele(rec:Variant, ibam:Bam, maxTotalReads:int,
if aln.base_quality_at(qoff - over-1).cint < minbsq:
continue
total_reads+=1
- if alleleCountTable.hasKey(cbt.get):
+ if alleleCountTable.hasKey(currentCB):
if aln.base_at(qoff - over-1) == rec.REF[0]:
- alleleCountTable[cbt.get].inc_count_cref
+ alleleCountTable[currentCB].inc_count_cref
continue
if aln.base_at(qoff - over-1) == rec_alt:
- alleleCountTable[cbt.get].inc_count_calt
+ alleleCountTable[currentCB].inc_count_calt
continue
else:
var new_snp = allele_expr(cref:0,calt:0)
- alleleCountTable[cbt.get] = new_snp
+ alleleCountTable[currentCB] = new_snp
if aln.base_at(qoff - over-1) == rec.REF[0]:
- alleleCountTable[cbt.get].inc_count_cref
+ alleleCountTable[currentCB].inc_count_cref
continue
if aln.base_at(qoff - over-1) == rec_alt:
- alleleCountTable[cbt.get].inc_count_calt
+ alleleCountTable[currentCB].inc_count_calt
continue
if total_reads > maxTotalReads or total_reads < minTotalReads:
return initTable[string,allele_expr]()
@@ -179,7 +184,6 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
outFileSNPanno = openFileStream(out_dir & chrom & "_snpAnnot.txt", fmWrite)
outFileTotalCountMtx = openFileStream(out_dir & chrom & "_totalCount.mtx", fmWrite)
outFileAltCountMtx = openFileStream(out_dir & chrom & "_altCount.mtx", fmWrite)
- # outHeaderMtx = openFileStream(out_dir & chrom & "_header.mtx", fmWrite)
outFileVStateMtx = openFileStream(out_dir & chrom & "_vi.mtx", fmWrite)
viSegmentInfo = openFileStream(out_dir & chrom & "_viSegInfo.txt", fmWrite)
#echo strm.readLine()
@@ -217,7 +221,6 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
# ac.tostring(acs)
## mindp, maxdp, they are values per cell
if (ac.cref+ac.calt) <= mindp or (ac.cref+ac.calt) >= maxdp: continue
-
var ithSperm = barcodeTable[bc]
## write to mtx Ref count
outFileTotalCountMtx.writeLine($snpIndex & " " & $(ithSperm+1) & " " & $(ac.cRef+ac.cAlt))
--
GitLab