Commit 460f8534 authored by Ruqian Lyu's avatar Ruqian Lyu
Browse files

add cell barcode option

parent 02a1aa86
Pipeline #6948 passed with stages
in 1 minute and 10 seconds
......@@ -59,13 +59,13 @@ 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,bulkBam:bool): Table[string,allele_expr] =
barcodeTable:OrderedTableRef,minbsq:int,bulkBam:bool,barcodeTag:string): Table[string,allele_expr] =
var alleleCountTable = initTable[string,allele_expr]()
var rec_alt:char
var total_reads=0
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 cbt = tag[string](aln, barcodeTag)
var currentCB: string
if cbt.isNone:
# echo "no cb "
......@@ -127,7 +127,7 @@ proc countAllele(rec:Variant, ibam:Bam, maxTotalReads:int,
proc sscocaller(threads:int, vcff:string, barcodeFile:string,
bamfile:string, out_dir:string, mapq:int,
minbsq:int, mintotal:int, maxtotal:int, mindp:int, maxdp:int,
thetaREF:float, thetaALT:float, cmPmb:float,s_Chrs:seq): int =
thetaREF:float, thetaALT:float, cmPmb:float,s_Chrs:seq,barcodeTag:string): int =
var
ibam:Bam
......@@ -202,7 +202,8 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
minbsq=minbsq,
maxTotalReads = maxtotal,
minTotalReads = mintotal,
bulkBam = bulkBam)
bulkBam = bulkBam,
barcodeTag = barcodeTag)
if alleleCountTable.len==0: continue
var rec_alt:char
......@@ -404,7 +405,7 @@ proc sscocaller(threads:int, vcff:string, barcodeFile:string,
return 0
when(isMainModule):
let version = "0.2.0"
let version = "0.2.1"
var doc = format("""
$version
......@@ -424,6 +425,7 @@ Arguments:
Options:
-t --threads <threads> number of BAM decompression threads [default: 4]
-cb --cellbarcode <cellbarcode> the cell barcode tag, by default it is CB
-MQ --minMAPQ <mapq> Minimum MAPQ for read filtering [default: 20]
-BQ --baseq <baseq> base quality threshold for a base to be used for counting [default: 13]
-CHR --chrom <chrom> the selected chromsome (whole genome if not supplied,separate by comma if multiple chroms)
......@@ -459,9 +461,12 @@ Options:
thetaREF:float
thetaALT:float
cmPmb:float
barcodeTag="CB"
if($args["--threads"] != "nil"):
threads = parse_int($args["--threads"])
if($args["--cellbarcode"] != "nil"):
barcodeTag = $args["--barcodeTag"]
else: threads = 4
if($args["--minDP"] != "nil"):
mindp = parse_int($args["--minDP"])
......@@ -530,5 +535,5 @@ Options:
#var args = commandLineParams()
discard sscocaller(threads, vcff, barcodeFile, bamfile,
out_dir, mapq, minbsq, mintotal,
maxtotal, mindp, maxdp, thetaREF, thetaALT, cmPmb,s_Chrs)
maxtotal, mindp, maxdp, thetaREF, thetaALT, cmPmb,s_Chrs,barcodeTag)
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