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

add bulk support

parent c17d69da
Pipeline #6676 passed with stages
in 36 seconds
......@@ -21,9 +21,6 @@ type
cref: int
calt: int
type
AlleleExpr = object
cRef: int
cAlt: int
# Nucleotide = enum
# A, C, T, G
ViState = enum
......@@ -42,16 +39,9 @@ type
## look up from current Sperm SNP index to SNPIndex
spermSnpIndexLookUp: Table[int,int]
SeqSpermViNodes = seq[SpermViNodes]
#var alexprSeq:seq[AlleleExpr]
proc inc_count_cref(a:allele_expr) = inc(a.cref)
proc inc_count_calt(a:allele_expr) = inc(a.calt)
proc reset_count(a:allele_expr) =
a.calt = 0
a.cref = 0
proc tostring(a:allele_expr, s:var string) {.inline.} =
s.set_len(0)
s.add("\t" & $a.cref & "," & $a.calt)
proc getEmission(thetaRef=0.1,thetaAlt=0.9,
cRef:int,cAlt:int): array[stateRef..stateAlt, float] =
......@@ -69,7 +59,7 @@ proc getTrans(pos1:int64,pos2:int64,cmPmb=0.1): float =
proc countAllele(rec:Variant, ibam:Bam, maxTotalReads:int,
minTotalReads:int,
chrom:string, mapq: int,
barcodeTable:OrderedTableRef,minbsq:int): Table[string,allele_expr] =
barcodeTable:OrderedTableRef,minbsq:int,bulkBam:bool): Table[string,allele_expr] =
var alleleCountTable = initTable[string,allele_expr]()
var rec_alt:char
var total_reads=0
......@@ -156,6 +146,7 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
kstr.m = 0
kstr.s = nil
var ithSperm=0
var bulkBam = false
## initiate the table with CB as has keys, allele counts (ref object) as elements
while hts_getline(hf, cint(10), addr kstr) > 0:
if $kstr.s[0] == "#":
......@@ -163,12 +154,14 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
var v = $kstr.s
discard barcodeTable.hasKeyOrPut(v, ithSperm)
ithSperm+=1
if barcodeTable.len == 1 and barcodeTable.hasKey("bulk"):
echo "running in bulk mode and all reads are regarded as from one sample/cell"
bulkBam = true
let initProb:array[stateRef..stateAlt, float]=[0.5,0.5]
## iterate through each selected chromosome
for chrom in s_Chrs:
# var scTable = initTable[string,SingleSpermRec]()
var snpIndex = 0
var spermIndex = 0
## number of non zeros
var nnsize = 0
var scSpermSeq:SeqSpermViNodes
......@@ -208,7 +201,8 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
barcodeTable=barcodeTable,
minbsq=minbsq,
maxTotalReads = maxtotal,
minTotalReads = mintotal)
minTotalReads = mintotal,
bulkBam = bulkBam)
if alleleCountTable.len==0: continue
var rec_alt:char
......@@ -410,7 +404,7 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
return 0
when(isMainModule):
let version = "0.1.0"
let version = "0.2.0"
var doc = format("""
$version
......@@ -444,7 +438,7 @@ Options:
-h --help show help
Examples
./sscocaller --threads 10 AAAGTAGCACGTCTCT-1.raw.bam AAAGTAGCACGTCTCT-1.raw.bam.dp3.alt.vcf.gz barcodeFile.tsv ./percell/ccsnp-
./sscocaller --threads 10 AAAGTAGCACGTCTCT-1.raw.bam AAAGTAGCACGTCTCT-1.raw.bam.dp3.alt.vcf.gz barcodeFile.tsv ./percell/ccsnp
""" % ["version", version])
let args = docopt(doc, version=version)
......
Markdown is supported
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