From c7fdf3ec56fb536476db6f429fc19c5038de479c Mon Sep 17 00:00:00 2001 From: Jody Phelan Date: Tue, 26 Nov 2024 10:23:40 +0100 Subject: [PATCH] add crude mpileup varaint caller for barcode --- pathogenprofiler/bam.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/pathogenprofiler/bam.py b/pathogenprofiler/bam.py index 84c8bb3..3c55906 100644 --- a/pathogenprofiler/bam.py +++ b/pathogenprofiler/bam.py @@ -301,20 +301,26 @@ def get_bed_gt(self,bed_file: str,ref_file: str,caller: str,platform: str): cmd = "freebayes -f %(ref_file)s -t %(bed_file)s %(prefix)s.tmp.bam --haplotype-length -1 | bcftools query -f '%%CHROM\\t%%POS\\t%%REF\\t%%ALT[\\t%%GT\\t%%AD]\\n'" % vars(self) elif caller == "bcftools": cmd = "bcftools mpileup -f %(ref_file)s -T %(bed_file)s %(prefix)s.tmp.bam -BI -a AD | bcftools call -mv | bcftools query -f '%%CHROM\\t%%POS\\t%%REF\\t%%ALT[\\t%%GT\\t%%AD]\\n'" % vars(self) + elif caller == "mpileup": + cmd = "bcftools mpileup -f %(ref_file)s -T %(bed_file)s %(prefix)s.tmp.bam -BI -a AD | bcftools query -f '%%CHROM\\t%%POS\\t%%REF\\t%%ALT[\\t%%AD]\\n'" % vars(self) else: cmd = "freebayes -f %(ref_file)s -t %(bed_file)s %(prefix)s.tmp.bam --haplotype-length -1 | bcftools query -f '%%CHROM\\t%%POS\\t%%REF\\t%%ALT[\\t%%GT\\t%%AD]\\n'" % vars(self) for l in cmd_out(cmd): # Chromosome 4348079 0/0 51 - chrom, pos, ref, alt, gt, ad = l.rstrip().split() + row = l.strip().split() + if len(row)==6: + chrom, pos, ref, alt, gt, ad = l.rstrip().split() + elif len(row)==5: + chrom, pos, ref, alt, ad = l.rstrip().split() + gt = None + p = GenomePosition(chrom=chrom,pos=int(pos)) d = {} alts = alt.split(",") ad = [int(x) for x in ad.split(",")] - # if gt == "0/0": - # d[ref] = ad[0] - # elif gt == "./.": - if gt == "./.": + + if gt and gt == "./.": d[ref] = 0 else: genotypes = list([ref]+alts)