feat: make sample_mixup_configurable

config/config.yaml

@@ -323,6 +323,9 @@ report_fusions:
 report_gene_fuse:
   min_unique_reads: 6
 
+sample_mixup_check:
+  match_cutoff: 0.7
+
 samtools_merge_bam:
   extra: "-c -p"
 

workflow/rules/bcftools.smk

@@ -41,6 +41,6 @@ rule bcftools_id_snps:
         "{input.bam} "
         "| bcftools call "
         "--skip-variants indels "
-        "-m -O v "
+        "-c -O v "
         "-o {output.vcf}) "
         "&> {log}"

workflow/rules/sample_mixup_check.smk

@@ -22,6 +22,7 @@ rule sample_mixup_check:
         mixup_report="qc/sample_mixup_check/sample_mixup_check.tsv",
     params:
         extra=config.get("sample_mixup_check", {}).get("extra", ""),
+        match_cutoff=config.get("sample_mixup_check", {}).get("match_cutoff", "0.7"),
     log:
         "twist_solid/sample_mixup_check/sample_mixup_check.tsv.log",
     benchmark:

workflow/schemas/config.schema.yaml

@@ -635,6 +635,9 @@ properties:
       extra:
         type: string
         description: parameters that should be forwarded
+      match_cutoff:
+        type: number
+        description: minimum fraction of matching SNP calls for an RNA and DNA matched sample
 
 required:
   - samples

workflow/scripts/sample_mixup_check.py

@@ -4,6 +4,7 @@
 vcf_dna_filenames = snakemake.input.id_snp_vcf_dna
 vcf_rna_filenames = snakemake.input.id_snp_vcf_rna
 report = open(snakemake.output.mixup_report, "w")
+match_cutoff = snakemake.params.match_cutoff
 
 vcf_dict_dna = {}
 dna_samples = {}
@@ -64,9 +65,9 @@ def read_vcf(vcf_filename, vcf_dict, samples):
         if rna_samples[rna_sample][dna_sample] > best_gt_match:
             best_dna_sample = dna_sample
             best_gt_match = rna_samples[rna_sample][dna_sample]
-    p_match = round(best_gt_match * 100 / 42.0, 1)
+    p_match = round(best_gt_match / 42.0, 1)
     report.write(f"{rna_sample}\t{best_dna_sample}\t{best_gt_match}\t{p_match}%\t")
-    if p_match > 80:
+    if p_match > match_cutoff:
         report.write(f"yes\n")
     else:
         report.write(f"no\n")