Commit c17d69da authored by Ruqian Lyu's avatar Ruqian Lyu
Browse files

add bulk sample support

parent d755de85
Pipeline #6674 failed with stages
in 33 seconds
......@@ -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))
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment