Skip to content

Commit

Permalink
fix: add decompress to manta vcf and annotation for normal ratio head…
Browse files Browse the repository at this point in the history
…er to Tonly samples
  • Loading branch information
elleira committed Nov 11, 2024
1 parent 56ebffa commit 969616b
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 30 deletions.
16 changes: 16 additions & 0 deletions workflow/Snakefile
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ ruleorder: filtering_bcftools_filter_include_region_snv_tn > bgzip_mutect
ruleorder: filtering_bcftools_filter_include_region_snv_t > bgzip_mutect
ruleorder: cnv_sv_manta_run_workflow_t > tabix_cnv
ruleorder: cnv_sv_manta_run_workflow_tn > tabix_cnv
ruleorder: filtering_bcftools_filter_include_region_sv > bgzip_cnv


aligner = config.get("aligner", None)
Expand Down Expand Up @@ -738,6 +739,21 @@ use rule tabix from misc as tabix_cnv with:
file="cnv_sv/[A-Za-z0-9-_./]+.vcf",


use rule bgzip from misc as decompress_manta with:
input:
"cnv_sv/manta_run_workflow_{analysis}/{sample}.ssa.include.{bed}.vcf.gz",
output:
"cnv_sv/manta_run_workflow_{analysis}/{sample}.ssa.include.{bed}.vcf",
params:
extra="-d"
log:
"cnv_sv/manta_run_workflow_{analysis}/{sample}.ssa.include.{bed}.vcf.log",
benchmark:
repeat("cnv_sv/manta_run_workflow_{analysis}/{sample}.ssa.include.{bed}.vcf.benchmark.tsv", config.get("bgzip", {}).get("benchmark_repeats", 1))
message:
"{rule}: decompress {input}"


use rule samtools_index from misc as misc_samtools_index with:
input:
cram="{file}.crumble.cram",
Expand Down
66 changes: 36 additions & 30 deletions workflow/scripts/annotate_normal_ratio.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,33 +19,39 @@
logging.debug("Creating header for " + snakemake.output.vcf)
vcf_out = VariantFile(snakemake.output.vcf, "w", header=new_header)

sample_tumor = [x for x in list(vcf_in.header.samples) if x.endswith("_T")][0]
sample_normal = [x for x in list(vcf_in.header.samples) if x.endswith("_N")][0]
logging.debug(f"{sample_tumor=}, {sample_normal=}")

logging.info("Starting to iterate through input.")
for record in vcf_in.fetch():
try:
n_freq = float(record.samples[sample_normal].get("AF")[0])
except TypeError:
n_freq = 0

try:
t_freq = float(record.samples[sample_tumor].get("AF")[0])
except TypeError:
t_freq = 0

if n_freq == 0:
logging.warning("Normal freq 0 or missing, setting ratio to 0")
ratio = 0
elif t_freq == 0:
logging.warning("Tumor freq 0 or missing, setting ratio to 100")
ratio = 100
else:
ratio = float(record.samples[sample_normal].get("AF")[0]) / float(record.samples[sample_tumor].get("AF")[0])

logging.debug(f"{record.samples[sample_normal].get('AF')[0]=}, {record.samples[sample_tumor].get('AF')[0]=}")
logging.debug(f"{ratio=}")
record.info["N_ratio"] = ratio

vcf_out.write(record)
if len(list(vcf_in.header.samples)) > 1:
logging.info("Two or more samples identified, proccessing as a TN vcf")
sample_tumor = [x for x in list(vcf_in.header.samples) if x.endswith("_T")][0]
sample_normal = [x for x in list(vcf_in.header.samples) if x.endswith("_N")][0]
logging.debug(f"{sample_tumor=}, {sample_normal=}")

logging.info("Starting to iterate through input.")
for record in vcf_in.fetch():
try:
n_freq = float(record.samples[sample_normal].get("AF")[0])
except TypeError:
n_freq = 0

try:
t_freq = float(record.samples[sample_tumor].get("AF")[0])
except TypeError:
t_freq = 0

if n_freq == 0:
logging.warning("Normal freq 0 or missing, setting ratio to 0")
ratio = 0
elif t_freq == 0:
logging.warning("Tumor freq 0 or missing, setting ratio to 100")
ratio = 100
else:
ratio = float(record.samples[sample_normal].get("AF")[0]) / float(record.samples[sample_tumor].get("AF")[0])

logging.debug(f"{record.samples[sample_normal].get('AF')[0]=}, {record.samples[sample_tumor].get('AF')[0]=}")
logging.debug(f"{ratio=}")
record.info["N_ratio"] = ratio

vcf_out.write(record)
else:
logging.info("Only one sample found in header processing as a T only vcf.")
for record in vcf_in.fetch():
vcf_out.write(record)

0 comments on commit 969616b

Please sign in to comment.