From bd9f9ec0dce191c42a6a3348da17945f30fc8096 Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:55:42 -0400 Subject: [PATCH 01/10] Create add_DP.wdl --- wdl/add_DP.wdl | 96 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 wdl/add_DP.wdl diff --git a/wdl/add_DP.wdl b/wdl/add_DP.wdl new file mode 100644 index 0000000..84dbac2 --- /dev/null +++ b/wdl/add_DP.wdl @@ -0,0 +1,96 @@ +version 1.0 + +workflow AddDepthFields { + input { + Array[File] sharded_vcfs + Array[File] sharded_vcf_indices + String output_prefix + File? reference_fasta + File? reference_index + } + + scatter (i in range(length(sharded_vcfs))) { + call AddDepthToVCF { + input: + vcf = sharded_vcfs[i], + vcf_index = sharded_vcf_indices[i], + shard_name = basename(sharded_vcfs[i], ".vcf.gz"), + reference_fasta = reference_fasta, + reference_index = reference_index + } + } + + output { + Array[File] vcfs_with_depth = AddDepthToVCF.output_vcf + Array[File] vcf_indices = AddDepthToVCF.output_vcf_index + } +} + +task AddDepthToVCF { + input { + File vcf + File vcf_index + String shard_name + File? reference_fasta + File? reference_index + + Int disk_size = ceil(size(vcf, "GB") * 3) + 20 + Int mem_gb = 16 + Int cpu = 2 + } + + command <<< + set -euo pipefail + + # Install tabix + apt-get update && apt-get install -y tabix + + # Create Python script to run Hail + cat << 'EOF' > add_depth.py +import hail as hl +import sys + +input_vcf = sys.argv[1] +output_vcf = sys.argv[2] + +hl.init(log="/dev/null") + +# Import VCF with GRCh38 reference genome and allow missing array elements +mt = hl.import_vcf(input_vcf, force_bgz=True, reference_genome='GRCh38', array_elements_required=False) + +# Add FORMAT/DP = sum(AD) +mt = mt.annotate_entries( + DP = hl.sum(mt.AD) +) + +# Add INFO/DP = sum of sample DP values (add to info struct) +mt = mt.annotate_rows( + info = mt.info.annotate(DP = hl.agg.sum(mt.DP)) +) + +# Export final VCF +header = hl.get_vcf_metadata(input_vcf) +hl.export_vcf(mt, output_vcf, metadata=header) +EOF + + # Run hail script + python3 add_depth.py "~{vcf}" "~{shard_name}.with_depth.vcf.bgz" + + # Rename and index the VCF + mv "~{shard_name}.with_depth.vcf.bgz" "~{shard_name}.with_depth.vcf.gz" + tabix -p vcf "~{shard_name}.with_depth.vcf.gz" + >>> + + output { + File output_vcf = "~{shard_name}.with_depth.vcf.gz" + File output_vcf_index = "~{shard_name}.with_depth.vcf.gz.tbi" + } + + runtime { + docker: "us.gcr.io/talkowski-sv-gnomad/shineren:hail-0.2.134" + memory: "~{mem_gb} GB" + disks: "local-disk ~{disk_size} HDD" + cpu: cpu + preemptible: 3 + } +} From c597abae57639277055af74f646b7742a1d18eb3 Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:56:32 -0400 Subject: [PATCH 02/10] Add WDL workflow for BAM to FASTQ conversion --- wdl/bam_to_fastq.wdl | 76 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 wdl/bam_to_fastq.wdl diff --git a/wdl/bam_to_fastq.wdl b/wdl/bam_to_fastq.wdl new file mode 100644 index 0000000..67461e9 --- /dev/null +++ b/wdl/bam_to_fastq.wdl @@ -0,0 +1,76 @@ +version 1.0 + +workflow BamToFastq { + input { + File input_reads + String output_prefix + + File? ref_fasta + File? ref_fasta_fai + Int bam_to_fastq_disk_gb + } + + call SamtoolsBamToFastq { + input: + input_reads = input_reads, + output_prefix = output_prefix, + ref_fasta = ref_fasta, + ref_fasta_fai = ref_fasta_fai, + bam_to_fastq_disk_gb = bam_to_fastq_disk_gb + } + + output { + File fastq_1 = SamtoolsBamToFastq.fastq_1 + File fastq_2 = SamtoolsBamToFastq.fastq_2 + } +} + +task SamtoolsBamToFastq { + input { + File input_reads + String output_prefix + File? ref_fasta + File? ref_fasta_fai + Int bam_to_fastq_disk_gb + } + + Boolean is_cram = sub(input_reads, ".*\\.", "") == "cram" + Boolean has_ref = defined(ref_fasta) + + command <<< + set -euo pipefail + + # 1. Determine if a reference is required (for CRAMs) + REF_ARGS="" + if ~{true="true" false="false" is_cram} ; then + if ! ~{true="true" false="false" has_ref} ; then + echo "Error: Reference FASTA is required for CRAM files." >&2 + exit 1 + fi + REF_ARGS="--reference ~{ref_fasta}" + fi + + # 2. Collate alignments by queryname + samtools collate -@ 4 $REF_ARGS -o collated.bam -T tmp_collate "~{input_reads}" + + # 3. Convert to paired compressed FASTQs + samtools fastq -@ 4 \ + -1 "~{output_prefix}_1.fastq.gz" \ + -2 "~{output_prefix}_2.fastq.gz" \ + -0 /dev/null -s /dev/null \ + collated.bam + >>> + + runtime { + docker: "staphb/samtools:1.17" + cpu: 4 + memory: "16 GB" + disks: "local-disk " + bam_to_fastq_disk_gb + " SSD" + preemptible: 1 + } + + output { + File fastq_1 = "~{output_prefix}_1.fastq.gz" + File fastq_2 = "~{output_prefix}_2.fastq.gz" + } +} From aff4ac12c0db78aa119c4fcb77a7f34d19a4540b Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:58:33 -0400 Subject: [PATCH 03/10] Add WDL workflow for computing het SNV AB stats with CNV validation using small variant VCF This WDL workflow computes heterozygosity allele balance (AB) statistics for samples, with optional validation against CNVs. It includes tasks for both genome-wide analysis and CNV-specific analysis, producing multiple output files. --- wdl/compute_het_AB_stats | 529 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 529 insertions(+) create mode 100644 wdl/compute_het_AB_stats diff --git a/wdl/compute_het_AB_stats b/wdl/compute_het_AB_stats new file mode 100644 index 0000000..7edda73 --- /dev/null +++ b/wdl/compute_het_AB_stats @@ -0,0 +1,529 @@ +version 1.0 + +workflow compute_het_ab_stats { + + input { + File vcf + File vcf_index + Int min_gq = 20 + Int min_ad_total = 10 + String docker + Int mem_gb + Int cpu = 4 + + # Optional: BED file of SVs/CNVs to validate with het AB + File? cnv_bed + } + + # Original genome-wide task (unchanged) + call ComputeHetAB { + input: + vcf = vcf, + vcf_index = vcf_index, + min_gq = min_gq, + min_ad_total = min_ad_total, + docker = docker, + mem_gb = mem_gb, + cpu = cpu + } + + # Optional CNV validation task + if (defined(cnv_bed)) { + call ComputeCNVHetAB { + input: + vcf = vcf, + vcf_index = vcf_index, + cnv_bed = select_first([cnv_bed]), + min_gq = min_gq, + min_ad_total = min_ad_total, + docker = docker, + mem_gb = mem_gb, + cpu = cpu + } + } + + output { + File het_ab_stats = ComputeHetAB.het_ab_stats + File? cnv_het_ab_stats = ComputeCNVHetAB.cnv_het_ab_stats + File? cnv_het_ab_raw = ComputeCNVHetAB.cnv_het_ab_raw + File? cnv_het_ab_variant_summary = ComputeCNVHetAB.cnv_het_ab_variant_summary + } + + meta { + author : "Generated WDL" + description : "Compute het AB statistics per sample per chrom for triploidy detection; optionally validate CNVs via het AB" + } +} + +# ── Original genome-wide task ───────────────────────────────────────────────── +task ComputeHetAB { + input { + File vcf + File vcf_index + Int min_gq + Int min_ad_total + String docker + Int mem_gb + Int cpu + Int disk_gb = ceil(size(vcf, "GiB") * 2) + 500 + } + + command <<< + set -euo pipefail + + bcftools view \ + --exclude-type other \ + --include 'FILTER="PASS" || FILTER="."' \ + --threads ~{cpu} \ + ~{vcf} \ + | bcftools query \ + -f '[%CHROM\t%SAMPLE\t%GT\t%GQ\t%AD\n]' \ + -o raw_genotypes.tsv + + echo "Total lines extracted: $(wc -l < raw_genotypes.tsv)" >&2 + head -3 raw_genotypes.tsv >&2 + + python3 -c " +import re, sys +from collections import defaultdict + +min_gq = ~{min_gq} +min_ad_total = ~{min_ad_total} + +DIPLOID_LO = 0.40 +DIPLOID_HI = 0.50 +TRIP_LO1 = 0.28 +TRIP_HI1 = 0.38 + +HET_RE = re.compile(r'^(0[/|][1-9][0-9]*|[1-9][0-9]*[/|]0)$') + +ab_vals = defaultdict(list) + +with open('raw_genotypes.tsv') as fh: + for line in fh: + parts = line.rstrip('\n').split('\t') + if len(parts) < 5: + continue + chrom, sample, gt, gq_s, ad_s = parts + + if not HET_RE.match(gt): + continue + + try: + if float(gq_s) < min_gq: + continue + except ValueError: + continue + + try: + ad = ad_s.split(',') + ref_ad = float(ad[0]) + alt_ad = float(ad[1]) + total = ref_ad + alt_ad + if total < min_ad_total: + continue + ab = min(ref_ad, alt_ad) / total + except (ValueError, IndexError): + continue + + ab_vals[(sample, chrom)].append(ab) + +with open('het_ab_stats.tsv', 'w') as out: + out.write('sample\tchrom\tn_het\tmean_ab\tmedian_ab\t' + 'frac_diploid\tfrac_trip\ttriploid_score\n') + for (sample, chrom), vals in sorted(ab_vals.items()): + n = len(vals) + if n == 0: + continue + mean_ab = sum(vals) / n + svals = sorted(vals) + median_ab = svals[n // 2] if n % 2 else (svals[n//2 - 1] + svals[n//2]) / 2 + + frac_diploid = sum(DIPLOID_LO <= v <= DIPLOID_HI for v in vals) / n + frac_trip = sum(TRIP_LO1 <= v <= TRIP_HI1 for v in vals) / n + triploid_score = frac_trip - frac_diploid + + out.write(f'{sample}\t{chrom}\t{n}\t{mean_ab:.4f}\t{median_ab:.4f}\t' + f'{frac_diploid:.4f}\t{frac_trip:.4f}\t{triploid_score:.4f}\n') + +print(f'Done. Wrote stats for {len(ab_vals)} sample-chrom combinations.', file=sys.stderr) +" + >>> + + output { + File het_ab_stats = "het_ab_stats.tsv" + } + + runtime { + docker : docker + memory : "~{mem_gb} GiB" + cpu : cpu + disks : "local-disk ~{disk_gb} HDD" + preemptible : 2 + } +} + +# ── CNV validation task ─────────────────────────────────────────────────────── +task ComputeCNVHetAB { + input { + File vcf + File vcf_index + File cnv_bed + Int min_gq + Int min_ad_total + String docker + Int mem_gb + Int cpu + Int disk_gb = ceil(size(vcf, "GiB") * 2) + 500 + } + + command <<< + set -euo pipefail + + # ── Step 1: extract unique samples from BED ─────────────────────────── + tail -n +2 ~{cnv_bed} \ + | awk '{print $6}' \ + | sort -u \ + > samples_of_interest.txt + + echo "Samples to query: $(cat samples_of_interest.txt | tr '\n' ' ')" >&2 + + # ── Step 2: build merged regions BED for bcftools query ─────────────── + tail -n +2 ~{cnv_bed} \ + | awk 'BEGIN{OFS="\t"} {print $1, $2+1, $3}' \ + | sort -k1,1 -k2,2n \ + | bedtools merge -i stdin \ + > merged_regions.bed + + echo "Merged query regions:" >&2 + cat merged_regions.bed >&2 + + # ── Step 3: extract genotypes for samples of interest in CNV regions ── + bcftools view \ + --exclude-type other \ + --include 'FILTER="PASS" || FILTER="."' \ + --samples-file samples_of_interest.txt \ + --regions-file merged_regions.bed \ + --threads ~{cpu} \ + ~{vcf} \ + | bcftools query \ + -f '[%CHROM\t%POS\t%SAMPLE\t%GT\t%GQ\t%AD\n]' \ + -o raw_cnv_genotypes.tsv + + echo "Genotype lines extracted: $(wc -l < raw_cnv_genotypes.tsv)" >&2 + + # ── Step 4: match SNVs to CNV intervals, compute AB, write outputs ──── + python3 -c " +import re, sys +from collections import defaultdict, OrderedDict + +min_gq = ~{min_gq} +min_ad_total = ~{min_ad_total} + +DIPLOID_LO = 0.40 +DIPLOID_HI = 0.50 +TRIP_LO = 0.28 +TRIP_HI = 0.38 + +HET_RE = re.compile(r'^(0[/|][1-9][0-9]*|[1-9][0-9]*[/|]0)$') +HOM_RE = re.compile(r'^([0-9]+)[/|]\1$') + +# ── Helper: nan-safe float formatter ───────────────────────────────────────── +def fmt(x): + return f'{x:.4f}' if x == x else 'NA' + +# ── Helper: compute interpretation call and stats for one sample-variant ────── +def interp_call(svtype, vals, n_het, n_hom, n_other): + n = len(vals) + if n == 0: + return 'NO_HET_SNVs', float('nan'), float('nan'), float('nan'), float('nan') + mean_ab = sum(vals) / n + svals = sorted(vals) + median_ab = svals[n // 2] if n % 2 else (svals[n//2-1] + svals[n//2]) / 2 + frac_dip = sum(DIPLOID_LO <= x <= DIPLOID_HI for x in vals) / n + frac_dup = sum(TRIP_LO <= x <= TRIP_HI for x in vals) / n + + if svtype == 'DUP': + if frac_dup >= 0.30 and median_ab < 0.42: + call = 'SUPPORTS_DUP' + elif frac_dip >= 0.50: + call = 'DIPLOID_AB_NO_DUP_SIGNAL' + else: + call = 'AMBIGUOUS' + elif svtype == 'DEL': + total_gts = n_het + n_hom + n_other + hom_rate = n_hom / total_gts if total_gts > 0 else 0 + if hom_rate >= 0.50: + call = 'SUPPORTS_DEL_HOM_DROPOUT' + elif n_het < 5: + call = 'TOO_FEW_HETS_FOR_DEL' + elif frac_dip >= 0.50: + call = 'DIPLOID_AB_NO_DEL_SIGNAL' + else: + call = 'AMBIGUOUS' + else: + call = 'UNKNOWN_SVTYPE' + + return call, mean_ab, median_ab, frac_dip, frac_dup + +# ── Load CNV BED ────────────────────────────────────────────────────────────── +# Flexible schema: +# 7-col: chrom start end name svtype sample CPX_INTERVALS +# 8-col: chrom start end variant_name svtype sample batch CN +# try/except on int conversion silently skips header and malformed rows +variants = [] +with open('~{cnv_bed}') as fh: + for line in fh: + line = line.rstrip('\n') + if not line: + continue + parts = line.split('\t') + if len(parts) < 6: + continue + try: + start = int(parts[1]) + end = int(parts[2]) + except ValueError: + continue + chrom = parts[0] + vname = parts[3] + svtype = parts[4].upper() + sample = parts[5] + batch = parts[6] if len(parts) > 6 else 'NA' + cn = parts[7] if len(parts) > 7 else 'NA' + variants.append({ + 'chrom' : chrom, + 'start' : start, + 'end' : end, + 'name' : vname, + 'svtype' : svtype, + 'sample' : sample, + 'batch' : batch, + 'cn' : cn, + }) + +print(f'Loaded {len(variants)} CNV records', file=sys.stderr) + +# Index variants by (chrom, sample) for fast overlap lookup +var_index = defaultdict(list) +for v in variants: + var_index[(v['chrom'], v['sample'])].append(v) + +# Per-(variant, sample) accumulators +raw_rows = [] +agg_data = defaultdict(lambda: {'ab_vals': [], 'n_hom': 0, 'n_het': 0, 'n_other': 0}) + +# Metadata map: (variant_name, sample) -> variant dict +var_meta = {} +for v in variants: + key = (v['name'], v['sample']) + var_meta[key] = v + +# ── Stream genotype file ────────────────────────────────────────────────────── +with open('raw_cnv_genotypes.tsv') as fh: + for line in fh: + parts = line.rstrip('\n').split('\t') + if len(parts) < 6: + continue + chrom, pos_s, sample, gt, gq_s, ad_s = parts + pos = int(pos_s) + + overlapping = [ + v for v in var_index.get((chrom, sample), []) + if v['start'] <= pos <= v['end'] + ] + if not overlapping: + continue + + try: + if float(gq_s) < min_gq: + continue + except ValueError: + continue + + try: + ad = ad_s.split(',') + ref_ad = float(ad[0]) + alt_ad = float(ad[1]) + total = ref_ad + alt_ad + except (ValueError, IndexError): + continue + + is_het = bool(HET_RE.match(gt)) + is_hom = bool(HOM_RE.match(gt)) + + for v in overlapping: + key = (v['name'], v['sample']) + + if is_het: + if total < min_ad_total: + continue + ab = min(ref_ad, alt_ad) / total + agg_data[key]['ab_vals'].append(ab) + agg_data[key]['n_het'] += 1 + raw_rows.append('\t'.join([ + v['name'], v['svtype'], v['sample'], v['batch'], v['cn'], + chrom, pos_s, gt, gq_s, + f'{ref_ad:.0f}', f'{alt_ad:.0f}', f'{ab:.4f}' + ])) + elif is_hom: + agg_data[key]['n_hom'] += 1 + else: + agg_data[key]['n_other'] += 1 + +# ── Output 1: raw het SNV calls ─────────────────────────────────────────────── +with open('cnv_het_ab_raw.tsv', 'w') as out: + out.write('variant_name\tsvtype\tsample\tbatch\tcn\t' + 'chrom\tpos\tgt\tgq\tref_ad\talt_ad\tab\n') + for row in raw_rows: + out.write(row + '\n') + +print(f'Raw file: {len(raw_rows)} het SNV rows', file=sys.stderr) + +# ── Output 2: per-(variant, sample) aggregated stats ───────────────────────── +with open('cnv_het_ab_stats.tsv', 'w') as out: + out.write( + 'variant_name\tsvtype\tsample\tbatch\tcn\tchrom\tstart\tend\t' + 'n_het\tn_hom\tn_other\t' + 'mean_ab\tmedian_ab\t' + 'frac_diploid\tfrac_dup_trip\t' + 'ab_interpretation\n' + ) + for key in sorted(var_meta.keys()): + v = var_meta[key] + data = agg_data[key] + vals = data['ab_vals'] + + call, mean_ab, median_ab, frac_dip, frac_dup = interp_call( + v['svtype'], vals, + data['n_het'], data['n_hom'], data['n_other'] + ) + + out.write( + f\"{v['name']}\t{v['svtype']}\t{v['sample']}\t{v['batch']}\t{v['cn']}\t\" + f\"{v['chrom']}\t{v['start']}\t{v['end']}\t\" + f\"{data['n_het']}\t{data['n_hom']}\t{data['n_other']}\t\" + f\"{fmt(mean_ab)}\t{fmt(median_ab)}\t\" + f\"{fmt(frac_dip)}\t{fmt(frac_dup)}\t\" + f\"{call}\n\" + ) + +print('Done writing per-sample stats.', file=sys.stderr) + +# ── Output 3: per-variant summary across all samples ───────────────────────── +# Unique variant names in BED order, preserving first-seen metadata +seen_vnames = OrderedDict() +for v in variants: + if v['name'] not in seen_vnames: + seen_vnames[v['name']] = v + +# Samples per variant +samples_per_variant = defaultdict(list) +for v in variants: + samples_per_variant[v['name']].append(v['sample']) + +with open('cnv_het_ab_variant_summary.tsv', 'w') as out: + out.write( + 'variant_name\tsvtype\tchrom\tstart\tend\t' + 'n_samples_in_bed\t' + 'n_samples_with_hets\t' + 'n_supports\t' + 'n_diploid_no_signal\t' + 'n_ambiguous\t' + 'n_no_het_snvs\t' + 'median_ab_across_samples\t' + 'mean_ab_across_samples\t' + 'frac_samples_supporting\t' + 'variant_verdict\n' + ) + + for vname, vmeta in seen_vnames.items(): + svtype = vmeta['svtype'] + bed_samples = samples_per_variant[vname] + n_in_bed = len(bed_samples) + + sample_median_abs = [] + n_supports = 0 + n_diploid = 0 + n_ambiguous = 0 + n_no_hets = 0 + n_with_hets = 0 + + for samp in bed_samples: + key = (vname, samp) + data = agg_data[key] + vals = data['ab_vals'] + + call, _, median_ab, _, _ = interp_call( + svtype, vals, + data['n_het'], data['n_hom'], data['n_other'] + ) + + if call in ('SUPPORTS_DUP', 'SUPPORTS_DEL_HOM_DROPOUT'): + n_supports += 1 + elif call in ('DIPLOID_AB_NO_DUP_SIGNAL', 'DIPLOID_AB_NO_DEL_SIGNAL'): + n_diploid += 1 + elif call in ('NO_HET_SNVs', 'TOO_FEW_HETS_FOR_DEL'): + n_no_hets += 1 + else: + n_ambiguous += 1 + + if len(vals) > 0: + n_with_hets += 1 + sample_median_abs.append(median_ab) + + # Median and mean of per-sample median ABs + if sample_median_abs: + sma = sorted(sample_median_abs) + nm = len(sma) + cross_median = sma[nm//2] if nm % 2 else (sma[nm//2-1] + sma[nm//2]) / 2 + cross_mean = sum(sma) / nm + else: + cross_median = float('nan') + cross_mean = float('nan') + + # Fraction of evaluable samples supporting the call + frac_support = (n_supports / n_with_hets) if n_with_hets > 0 else float('nan') + + # Variant-level verdict + if n_with_hets == 0: + verdict = 'NO_EVALUABLE_SAMPLES' + elif frac_support >= 0.50: + verdict = 'LIKELY_REAL' + elif frac_support <= 0.15: + verdict = 'LIKELY_ARTIFACT' + else: + verdict = 'UNCERTAIN' + + out.write( + f'{vname}\t{svtype}\t{vmeta[\"chrom\"]}\t{vmeta[\"start\"]}\t{vmeta[\"end\"]}\t' + f'{n_in_bed}\t' + f'{n_with_hets}\t' + f'{n_supports}\t' + f'{n_diploid}\t' + f'{n_ambiguous}\t' + f'{n_no_hets}\t' + f'{fmt(cross_median)}\t' + f'{fmt(cross_mean)}\t' + f'{fmt(frac_support)}\t' + f'{verdict}\n' + ) + +print('Done writing variant summary.', file=sys.stderr) +" + >>> + + output { + File cnv_het_ab_raw = "cnv_het_ab_raw.tsv" + File cnv_het_ab_stats = "cnv_het_ab_stats.tsv" + File cnv_het_ab_variant_summary = "cnv_het_ab_variant_summary.tsv" + } + + runtime { + docker : docker + memory : "~{mem_gb} GiB" + cpu : cpu + disks : "local-disk ~{disk_gb} HDD" + preemptible : 2 + } +} From d4ef9459b42ab35d11111dc0c861034dda93a1bf Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 16:59:18 -0400 Subject: [PATCH 04/10] Add workflow to extract read group information from BAM/CRAM --- wdl/extract_read_group_info | 145 ++++++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) create mode 100644 wdl/extract_read_group_info diff --git a/wdl/extract_read_group_info b/wdl/extract_read_group_info new file mode 100644 index 0000000..360829c --- /dev/null +++ b/wdl/extract_read_group_info @@ -0,0 +1,145 @@ +version 1.0 + +workflow ExtractReadGroup { + meta { + description: "Extract all read group information from a BAM or CRAM file" + } + + input { + File input_bam + File? ref_fasta + File? ref_fasta_fai + + Int cpu = 2 + Int mem_gb = 4 + Int disk_gb = 50 + } + + call GetReadGroups { + input: + input_bam = input_bam, + ref_fasta = ref_fasta, + ref_fasta_fai = ref_fasta_fai, + cpu = cpu, + mem_gb = mem_gb, + disk_gb = disk_gb + } + + output { + File read_group_file = GetReadGroups.read_group_file + Array[String] read_groups = GetReadGroups.read_groups + + Array[String] rg_ids = GetReadGroups.rg_ids + Array[String] rg_lbs = GetReadGroups.rg_lbs + Array[String] rg_pus = GetReadGroups.rg_pus + Array[String] rg_sms = GetReadGroups.rg_sms + Array[String] rg_pls = GetReadGroups.rg_pls + Array[String] rg_cns = GetReadGroups.rg_cns + Array[String] rg_pms = GetReadGroups.rg_pms + Array[String] rg_dts = GetReadGroups.rg_dts + } +} + +task GetReadGroups { + input { + File input_bam + File? ref_fasta + File? ref_fasta_fai + Int cpu = 2 + Int mem_gb = 4 + Int disk_gb = 50 + } + + Boolean has_ref = defined(ref_fasta) + + command <<< + set -euo pipefail + + REF_ARGS="" + if ~{true="true" false="false" has_ref} ; then + REF_ARGS="--reference ~{ref_fasta}" + fi + + # Extract all @RG lines from header + samtools view -H $REF_ARGS ~{input_bam} \ + | grep "^@RG" \ + > read_group.txt + + if [ ! -s read_group.txt ]; then + echo "Error: no read group found in ~{input_bam}" >&2 + exit 1 + fi + + echo "Found $(wc -l < read_group.txt) read group(s):" + cat read_group.txt + + # Helper: extract a tag from a given RG line + extract_tag() { + local line="$1" + local tag="$2" + echo "$line" | grep -oP "(?<=\t${tag}:)[^\t]+" || true + } + + # Clear output files + > rg_escaped.tsv + > rg_ids.txt + > rg_lbs.txt + > rg_pus.txt + > rg_sms.txt + > rg_pls.txt + > rg_cns.txt + > rg_pms.txt + > rg_dts.txt + + # Iterate over every RG line + while IFS= read -r rg_line; do + extract_tag "$rg_line" "ID" >> rg_ids.txt + extract_tag "$rg_line" "LB" >> rg_lbs.txt + extract_tag "$rg_line" "PU" >> rg_pus.txt + extract_tag "$rg_line" "SM" >> rg_sms.txt + extract_tag "$rg_line" "PL" >> rg_pls.txt + extract_tag "$rg_line" "CN" >> rg_cns.txt + extract_tag "$rg_line" "PM" >> rg_pms.txt + extract_tag "$rg_line" "DT" >> rg_dts.txt + + # Escaped full RG string (tabs -> \t) for use as bwa-mem2 -R argument + echo "$rg_line" \ + | tr '\t' '\n' \ + | paste -sd '\t' - \ + | sed 's/\t/\\t/g' \ + >> rg_escaped.tsv + done < read_group.txt + + echo "--- Parsed RG fields ---" + echo "IDs: $(cat rg_ids.txt | tr '\n' ' ')" + echo "LBs: $(cat rg_lbs.txt | tr '\n' ' ')" + echo "PUs: $(cat rg_pus.txt | tr '\n' ' ')" + echo "SMs: $(cat rg_sms.txt | tr '\n' ' ')" + echo "PLs: $(cat rg_pls.txt | tr '\n' ' ')" + echo "CNs: $(cat rg_cns.txt | tr '\n' ' ')" + echo "PMs: $(cat rg_pms.txt | tr '\n' ' ')" + echo "DTs: $(cat rg_dts.txt | tr '\n' ' ')" + >>> + + output { + File read_group_file = "read_group.txt" + Array[String] read_groups = read_lines("rg_escaped.tsv") + + Array[String] rg_ids = read_lines("rg_ids.txt") + Array[String] rg_lbs = read_lines("rg_lbs.txt") + Array[String] rg_pus = read_lines("rg_pus.txt") + Array[String] rg_sms = read_lines("rg_sms.txt") + Array[String] rg_pls = read_lines("rg_pls.txt") + Array[String] rg_cns = read_lines("rg_cns.txt") + Array[String] rg_pms = read_lines("rg_pms.txt") + Array[String] rg_dts = read_lines("rg_dts.txt") + } + + runtime { + docker: "staphb/samtools:1.17" + cpu: cpu + memory: mem_gb + " GB" + disks: "local-disk " + disk_gb + " HDD" + preemptible: 1 + } +} From c8f5e28d5427038601a470beb1122fe669a41ad2 Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:00:55 -0400 Subject: [PATCH 05/10] Add WDL workflow for filtering SVs in VCF files for association analysis --- wdl/Filter_SVs_for_Association_Analysis.wdl | 461 ++++++++++++++++++++ 1 file changed, 461 insertions(+) create mode 100644 wdl/Filter_SVs_for_Association_Analysis.wdl diff --git a/wdl/Filter_SVs_for_Association_Analysis.wdl b/wdl/Filter_SVs_for_Association_Analysis.wdl new file mode 100644 index 0000000..8126dd0 --- /dev/null +++ b/wdl/Filter_SVs_for_Association_Analysis.wdl @@ -0,0 +1,461 @@ +version 1.0 + +workflow VcfToBed { + + input { + File vcf_file + String cohort_prefix + String variant_interpretation_docker + File outlier_samples_file + File ped_file + + Array[String] svtypes = ["DEL", "DUP", "CNV"] + Float gnomad_af_threshold = 0.01 + Int min_proband_ac = 1 + Int max_proband_ac = 5 + Boolean filter_predicted_lof = true + Boolean filter_predicted_intragenic_exon_dup = true + Boolean filter_predicted_copy_gain = true + Boolean filter_predicted_partial_exon_dup = true + } + + call FilterAndSplitVcf { + input: + vcf_file = vcf_file, + cohort_prefix = cohort_prefix, + variant_interpretation_docker = variant_interpretation_docker, + outlier_samples_file = outlier_samples_file, + ped_file = ped_file, + svtypes = svtypes, + gnomad_af_threshold = gnomad_af_threshold, + min_proband_ac = min_proband_ac, + max_proband_ac = max_proband_ac, + filter_predicted_lof = filter_predicted_lof, + filter_predicted_intragenic_exon_dup = filter_predicted_intragenic_exon_dup, + filter_predicted_copy_gain = filter_predicted_copy_gain, + filter_predicted_partial_exon_dup = filter_predicted_partial_exon_dup + } + + output { + # Outlier-exclusive outputs (main PROBAND_AC branch) + File outlier_vcf = FilterAndSplitVcf.outlier_vcf + File outlier_vcf_tbi = FilterAndSplitVcf.outlier_vcf_tbi + File outlier_bed = FilterAndSplitVcf.outlier_bed + + # Non-outlier outputs (main PROBAND_AC branch, outlier samples fully excluded) + File nonoutlier_vcf = FilterAndSplitVcf.nonoutlier_vcf + File nonoutlier_vcf_tbi = FilterAndSplitVcf.nonoutlier_vcf_tbi + File nonoutlier_bed = FilterAndSplitVcf.nonoutlier_bed + + # All-samples BED (main PROBAND_AC branch) + File all_samples_bed = FilterAndSplitVcf.all_samples_bed + + # Variant ID lists (main branch) + File ids_outlier_exclusive = FilterAndSplitVcf.ids_outlier_exclusive + File ids_in_nonoutliers = FilterAndSplitVcf.ids_in_nonoutliers + + # Parent-only outputs: PROBAND_AC=0 branch, split by outlier status + File no_proband_carriers_outlier_vcf = FilterAndSplitVcf.no_proband_carriers_outlier_vcf + File no_proband_carriers_outlier_vcf_tbi = FilterAndSplitVcf.no_proband_carriers_outlier_vcf_tbi + File no_proband_carriers_outlier_bed = FilterAndSplitVcf.no_proband_carriers_outlier_bed + File no_proband_carriers_nonoutlier_vcf = FilterAndSplitVcf.no_proband_carriers_nonoutlier_vcf + File no_proband_carriers_nonoutlier_vcf_tbi = FilterAndSplitVcf.no_proband_carriers_nonoutlier_vcf_tbi + File no_proband_carriers_nonoutlier_bed = FilterAndSplitVcf.no_proband_carriers_nonoutlier_bed + File ids_no_proband_carriers_outlier = FilterAndSplitVcf.ids_no_proband_carriers_outlier + File ids_no_proband_carriers_nonoutlier = FilterAndSplitVcf.ids_no_proband_carriers_nonoutlier + File proband_ids = FilterAndSplitVcf.proband_ids + } +} + +task FilterAndSplitVcf { + + input { + File vcf_file + String cohort_prefix + String variant_interpretation_docker + File outlier_samples_file + File ped_file + Array[String] svtypes = ["DEL", "DUP", "CNV"] + Float gnomad_af_threshold = 0.01 + Int min_proband_ac = 1 + Int max_proband_ac = 5 + Boolean filter_predicted_lof = true + Boolean filter_predicted_intragenic_exon_dup = true + Boolean filter_predicted_copy_gain = true + Boolean filter_predicted_partial_exon_dup = true + Int mem_gb = 8 + Int cpu_cores = 2 + Int preemptible_tries = 3 + Int max_retries = 1 + Int boot_disk_gb = 10 + } + + Float input_size = size(vcf_file, "GB") + Int disk_gb = ceil(10.0 + input_size * 10.0) + + runtime { + memory: "~{mem_gb} GB" + disks: "local-disk ~{disk_gb} HDD" + cpu: cpu_cores + preemptible: preemptible_tries + maxRetries: max_retries + docker: variant_interpretation_docker + bootDiskSizeGb: boot_disk_gb + } + + command <<< + set -eou pipefail + + VCF="~{vcf_file}" + OUTLIERS="~{outlier_samples_file}" + PED="~{ped_file}" + PREFIX="~{cohort_prefix}" + + # --------------------------------------------------------------- + # Step 0: Extract proband IDs from PED file + # PED format: family_id, sample_id, father_id, mother_id, + # sex, affected (2=affected/proband) + # Primary: affected status = 2 + # Fallback: samples with both parents listed (col3/col4 != 0) + # --------------------------------------------------------------- + + awk '$6 == 2 {print $2}' "$PED" | sort -u > proband_ids.txt + + if [ ! -s proband_ids.txt ]; then + echo "WARNING: No affected=2 samples found in PED; falling back to samples with both parents listed" + awk '$3 != "0" && $3 != "" && $4 != "0" && $4 != "" {print $2}' "$PED" | sort -u > proband_ids.txt + fi + + echo "Probands identified from PED: $(wc -l < proband_ids.txt)" + + # --------------------------------------------------------------- + # Step 1: Restrict to requested SVTYPEs (all FILTER statuses) + # --------------------------------------------------------------- + SVTYPE_ARRAY=(~{sep=" " svtypes}) + + SVTYPE_CLAUSES=() + for ST in "${SVTYPE_ARRAY[@]}"; do + SVTYPE_CLAUSES+=("SVTYPE=\"${ST}\"") + done + SVTYPE_EXPR=$(printf '%s | ' "${SVTYPE_CLAUSES[@]}") + SVTYPE_EXPR="(${SVTYPE_EXPR% | })" + + bcftools view \ + -f "" \ + -i "${SVTYPE_EXPR}" \ + "$VCF" \ + -Oz -o svtype_filtered.vcf.gz + + bcftools index -t svtype_filtered.vcf.gz + + # --------------------------------------------------------------- + # Step 2: Build shared functional + AF expression + # PROBAND_AC thresholds are applied separately per branch + # --------------------------------------------------------------- + FUNC_CLAUSES=() + + if [ "~{filter_predicted_lof}" = "true" ]; then + FUNC_CLAUSES+=('INFO/PREDICTED_LOF!="."') + fi + if [ "~{filter_predicted_intragenic_exon_dup}" = "true" ]; then + FUNC_CLAUSES+=('INFO/PREDICTED_INTRAGENIC_EXON_DUP!="."') + fi + if [ "~{filter_predicted_copy_gain}" = "true" ]; then + FUNC_CLAUSES+=('INFO/PREDICTED_COPY_GAIN!="."') + fi + if [ "~{filter_predicted_partial_exon_dup}" = "true" ]; then + FUNC_CLAUSES+=('INFO/PREDICTED_PARTIAL_EXON_DUP!="."') + fi + + if [ ${#FUNC_CLAUSES[@]} -eq 0 ]; then + FUNC_EXPR="1=1" + else + FUNC_EXPR=$(printf '%s | ' "${FUNC_CLAUSES[@]}") + FUNC_EXPR="(${FUNC_EXPR% | })" + fi + + AF_EXPR="INFO/gnomad_v4.1_sv_AF<~{gnomad_af_threshold}" + + # Main branch: variants with proband carriers (PROBAND_AC min-max) + FILTER_EXPR="${FUNC_EXPR} & INFO/PROBAND_AC>=~{min_proband_ac} & INFO/PROBAND_AC<=~{max_proband_ac} & ${AF_EXPR}" + + # Parent-only branch: variants with zero proband carriers (PROBAND_AC=0) + PARENT_FILTER_EXPR="${FUNC_EXPR} & INFO/PROBAND_AC=0 & ${AF_EXPR}" + + echo "SVTYPE filter: ${SVTYPE_EXPR}" + echo "Main variant filter: ${FILTER_EXPR}" + echo "Parent variant filter: ${PARENT_FILTER_EXPR}" + + # --------------------------------------------------------------- + # Step 3a: Main branch — apply PROBAND_AC min-max filter + # --------------------------------------------------------------- + bcftools view \ + -f "" \ + -i "${FILTER_EXPR}" \ + svtype_filtered.vcf.gz \ + -Oz -o filtered.vcf.gz + + bcftools index -t filtered.vcf.gz + + echo "Main branch variants after filtering: $(bcftools view -H filtered.vcf.gz | wc -l)" + + # --------------------------------------------------------------- + # Step 3b: Parent-only branch — apply PROBAND_AC=0 filter + # --------------------------------------------------------------- + bcftools view \ + -f "" \ + -i "${PARENT_FILTER_EXPR}" \ + svtype_filtered.vcf.gz \ + -Oz -o filtered_parent_only.vcf.gz + + bcftools index -t filtered_parent_only.vcf.gz + + echo "Parent-only branch variants after filtering: $(bcftools view -H filtered_parent_only.vcf.gz | wc -l)" + + # --------------------------------------------------------------- + # Step 4: Main branch — determine outlier-exclusive vs non-outlier + # variant IDs from filtered.vcf.gz + # --------------------------------------------------------------- + + # Variant IDs carried by at least one outlier sample + bcftools view -f "" -S "$OUTLIERS" filtered.vcf.gz \ + | bcftools view -f "" -i 'GT="alt"' \ + | bcftools query -f '%ID\n' \ + | sort -u > ids_in_outliers.txt + + # Variant IDs carried by at least one non-outlier sample + bcftools view -f "" -S ^"$OUTLIERS" filtered.vcf.gz \ + | bcftools view -f "" -i 'GT="alt"' \ + | bcftools query -f '%ID\n' \ + | sort -u > ids_in_nonoutliers.txt + + # Outlier-exclusive: in outliers but zero non-outlier carriers + comm -23 ids_in_outliers.txt ids_in_nonoutliers.txt > ids_outlier_exclusive.txt + + echo "Variants carried by outlier samples: $(wc -l < ids_in_outliers.txt)" + echo "Variants carried by non-outlier samples: $(wc -l < ids_in_nonoutliers.txt)" + echo "Variants exclusive to outlier samples: $(wc -l < ids_outlier_exclusive.txt)" + + # --------------------------------------------------------------- + # Step 5: Main branch — generate outlier-exclusive and non-outlier VCFs + # --------------------------------------------------------------- + + # Outlier-exclusive VCF — variants with zero non-outlier carriers + bcftools view -f "" \ + -i "ID=@ids_outlier_exclusive.txt" \ + -S "$OUTLIERS" \ + filtered.vcf.gz \ + -Oz -o "${PREFIX}.outlier_samples.vcf.gz" + bcftools index -t "${PREFIX}.outlier_samples.vcf.gz" + + # Non-outlier VCF — variants carried by at least one non-outlier, + # samples column restricted to non-outliers only + bcftools view -f "" \ + -i "ID=@ids_in_nonoutliers.txt" \ + -S ^"$OUTLIERS" \ + filtered.vcf.gz \ + -Oz -o "${PREFIX}.nonoutlier_samples.vcf.gz" + bcftools index -t "${PREFIX}.nonoutlier_samples.vcf.gz" + + # --------------------------------------------------------------- + # Step 6: Parent-only branch — split by outlier status + # Uses filtered_parent_only.vcf.gz (PROBAND_AC=0) + # --------------------------------------------------------------- + + # Parent-only variant IDs carried by at least one outlier sample + bcftools view -f "" -S "$OUTLIERS" filtered_parent_only.vcf.gz \ + | bcftools view -f "" -i 'GT="alt"' \ + | bcftools query -f '%ID\n' \ + | sort -u > ids_parent_only_in_outliers.txt + + # Parent-only variant IDs carried by at least one non-outlier sample + bcftools view -f "" -S ^"$OUTLIERS" filtered_parent_only.vcf.gz \ + | bcftools view -f "" -i 'GT="alt"' \ + | bcftools query -f '%ID\n' \ + | sort -u > ids_parent_only_in_nonoutliers.txt + + echo "Parent-only variants in outlier samples: $(wc -l < ids_parent_only_in_outliers.txt)" + echo "Parent-only variants in non-outlier samples: $(wc -l < ids_parent_only_in_nonoutliers.txt)" + + # Outlier parent-only VCF — samples restricted to outliers only + bcftools view -f "" \ + -i "ID=@ids_parent_only_in_outliers.txt" \ + -S "$OUTLIERS" \ + filtered_parent_only.vcf.gz \ + -Oz -o "${PREFIX}.no_proband_carriers_outlier.vcf.gz" + bcftools index -t "${PREFIX}.no_proband_carriers_outlier.vcf.gz" + + # Nonoutlier parent-only VCF — samples restricted to non-outliers only + bcftools view -f "" \ + -i "ID=@ids_parent_only_in_nonoutliers.txt" \ + -S ^"$OUTLIERS" \ + filtered_parent_only.vcf.gz \ + -Oz -o "${PREFIX}.no_proband_carriers_nonoutlier.vcf.gz" + bcftools index -t "${PREFIX}.no_proband_carriers_nonoutlier.vcf.gz" + + # --------------------------------------------------------------- + # Step 7: Convert all VCFs to BED + # --------------------------------------------------------------- + svtk vcf2bed --include-filters -i ALL \ + filtered.vcf.gz "${PREFIX}.all_samples.bed.gz" + + svtk vcf2bed --include-filters -i ALL \ + "${PREFIX}.outlier_samples.vcf.gz" "${PREFIX}.outlier_samples.bed.gz" + + svtk vcf2bed --include-filters -i ALL \ + "${PREFIX}.nonoutlier_samples.vcf.gz" "${PREFIX}.nonoutlier_samples.bed.gz" + + svtk vcf2bed --include-filters -i ALL \ + "${PREFIX}.no_proband_carriers_outlier.vcf.gz" "${PREFIX}.no_proband_carriers_outlier.bed.gz" + + svtk vcf2bed --include-filters -i ALL \ + "${PREFIX}.no_proband_carriers_nonoutlier.vcf.gz" "${PREFIX}.no_proband_carriers_nonoutlier.bed.gz" + + # --------------------------------------------------------------- + # Step 8: Build per-subset carrier maps and fix samples column in BEDs + # - Main outlier BED: outlier carriers only + # - Main nonoutlier BED: non-outlier carriers only + # - Parent-only outlier BED: outlier carriers only (from parent-only VCF) + # - Parent-only nonoutlier BED: non-outlier carriers only (from parent-only VCF) + # --------------------------------------------------------------- + + # Main branch: outlier-only carrier map + bcftools view -f "" -S "$OUTLIERS" filtered.vcf.gz \ + | bcftools query -f '[%ID\t%SAMPLE\n]' -i 'GT="alt"' \ + | awk '{carriers[$1] = (carriers[$1] == "" ? $2 : carriers[$1] "," $2)} + END {for (id in carriers) print id"\t"carriers[id]}' \ + | sort > outlier_carrier_map.txt + + echo "Main outlier carrier map: $(wc -l < outlier_carrier_map.txt) variants" + + # Main branch: non-outlier-only carrier map + bcftools view -f "" -S ^"$OUTLIERS" filtered.vcf.gz \ + | bcftools query -f '[%ID\t%SAMPLE\n]' -i 'GT="alt"' \ + | awk '{carriers[$1] = (carriers[$1] == "" ? $2 : carriers[$1] "," $2)} + END {for (id in carriers) print id"\t"carriers[id]}' \ + | sort > nonoutlier_carrier_map.txt + + echo "Main nonoutlier carrier map: $(wc -l < nonoutlier_carrier_map.txt) variants" + + # Parent-only branch: outlier-only carrier map + bcftools view -f "" -S "$OUTLIERS" filtered_parent_only.vcf.gz \ + | bcftools query -f '[%ID\t%SAMPLE\n]' -i 'GT="alt"' \ + | awk '{carriers[$1] = (carriers[$1] == "" ? $2 : carriers[$1] "," $2)} + END {for (id in carriers) print id"\t"carriers[id]}' \ + | sort > parent_outlier_carrier_map.txt + + echo "Parent-only outlier carrier map: $(wc -l < parent_outlier_carrier_map.txt) variants" + + # Parent-only branch: non-outlier-only carrier map + bcftools view -f "" -S ^"$OUTLIERS" filtered_parent_only.vcf.gz \ + | bcftools query -f '[%ID\t%SAMPLE\n]' -i 'GT="alt"' \ + | awk '{carriers[$1] = (carriers[$1] == "" ? $2 : carriers[$1] "," $2)} + END {for (id in carriers) print id"\t"carriers[id]}' \ + | sort > parent_nonoutlier_carrier_map.txt + + echo "Parent-only nonoutlier carrier map: $(wc -l < parent_nonoutlier_carrier_map.txt) variants" + + # Fix BED samples columns using the appropriate carrier map per BED + python3 << 'PYEOF' +import gzip +import csv + +def load_carrier_map(path): + carrier_map = {} + with open(path) as f: + for line in f: + parts = line.strip().split("\t") + if len(parts) == 2: + carrier_map[parts[0]] = parts[1] + print(f"Loaded carrier map from {path}: {len(carrier_map)} variants") + return carrier_map + +def fix_bed(bed_in, bed_out, carrier_map): + n_fixed = 0 + n_total = 0 + with gzip.open(bed_in, 'rt') as fin, gzip.open(bed_out, 'wt') as fout: + reader = csv.reader(fin, delimiter='\t') + writer = csv.writer(fout, delimiter='\t', lineterminator='\n') + for i, row in enumerate(reader): + if i == 0: + writer.writerow(row) + continue + n_total += 1 + variant_id = row[3] # name column (0-indexed) + if variant_id in carrier_map: + row[5] = carrier_map[variant_id] # replace samples column + n_fixed += 1 + writer.writerow(row) + print(f"{bed_in}: fixed {n_fixed}/{n_total} rows") + +outlier_map = load_carrier_map("outlier_carrier_map.txt") +nonoutlier_map = load_carrier_map("nonoutlier_carrier_map.txt") +parent_outlier_map = load_carrier_map("parent_outlier_carrier_map.txt") +parent_nonoutlier_map = load_carrier_map("parent_nonoutlier_carrier_map.txt") + +fix_bed( + "~{cohort_prefix}.outlier_samples.bed.gz", + "~{cohort_prefix}.outlier_samples.fixed.bed.gz", + outlier_map +) +fix_bed( + "~{cohort_prefix}.nonoutlier_samples.bed.gz", + "~{cohort_prefix}.nonoutlier_samples.fixed.bed.gz", + nonoutlier_map +) +fix_bed( + "~{cohort_prefix}.no_proband_carriers_outlier.bed.gz", + "~{cohort_prefix}.no_proband_carriers_outlier.fixed.bed.gz", + parent_outlier_map +) +fix_bed( + "~{cohort_prefix}.no_proband_carriers_nonoutlier.bed.gz", + "~{cohort_prefix}.no_proband_carriers_nonoutlier.fixed.bed.gz", + parent_nonoutlier_map +) + +print("Done.") +PYEOF + + # Replace originals with fixed versions + mv "${PREFIX}.outlier_samples.fixed.bed.gz" "${PREFIX}.outlier_samples.bed.gz" + mv "${PREFIX}.nonoutlier_samples.fixed.bed.gz" "${PREFIX}.nonoutlier_samples.bed.gz" + mv "${PREFIX}.no_proband_carriers_outlier.fixed.bed.gz" "${PREFIX}.no_proband_carriers_outlier.bed.gz" + mv "${PREFIX}.no_proband_carriers_nonoutlier.fixed.bed.gz" "${PREFIX}.no_proband_carriers_nonoutlier.bed.gz" + + # Rename parent-only ID lists to final output names + mv ids_parent_only_in_outliers.txt ids_no_proband_carriers_outlier.txt + mv ids_parent_only_in_nonoutliers.txt ids_no_proband_carriers_nonoutlier.txt + + >>> + + output { + # Outlier-exclusive outputs (main PROBAND_AC branch) + File outlier_vcf = "~{cohort_prefix}.outlier_samples.vcf.gz" + File outlier_vcf_tbi = "~{cohort_prefix}.outlier_samples.vcf.gz.tbi" + File outlier_bed = "~{cohort_prefix}.outlier_samples.bed.gz" + + # Non-outlier outputs (main PROBAND_AC branch, outlier samples fully excluded) + File nonoutlier_vcf = "~{cohort_prefix}.nonoutlier_samples.vcf.gz" + File nonoutlier_vcf_tbi = "~{cohort_prefix}.nonoutlier_samples.vcf.gz.tbi" + File nonoutlier_bed = "~{cohort_prefix}.nonoutlier_samples.bed.gz" + + # All-samples BED (main PROBAND_AC branch) + File all_samples_bed = "~{cohort_prefix}.all_samples.bed.gz" + + # Variant ID lists (main branch) + File ids_outlier_exclusive = "ids_outlier_exclusive.txt" + File ids_in_nonoutliers = "ids_in_nonoutliers.txt" + + # Parent-only outputs: PROBAND_AC=0 branch, split by outlier status + File no_proband_carriers_outlier_vcf = "~{cohort_prefix}.no_proband_carriers_outlier.vcf.gz" + File no_proband_carriers_outlier_vcf_tbi = "~{cohort_prefix}.no_proband_carriers_outlier.vcf.gz.tbi" + File no_proband_carriers_outlier_bed = "~{cohort_prefix}.no_proband_carriers_outlier.bed.gz" + File no_proband_carriers_nonoutlier_vcf = "~{cohort_prefix}.no_proband_carriers_nonoutlier.vcf.gz" + File no_proband_carriers_nonoutlier_vcf_tbi = "~{cohort_prefix}.no_proband_carriers_nonoutlier.vcf.gz.tbi" + File no_proband_carriers_nonoutlier_bed = "~{cohort_prefix}.no_proband_carriers_nonoutlier.bed.gz" + File ids_no_proband_carriers_outlier = "ids_no_proband_carriers_outlier.txt" + File ids_no_proband_carriers_nonoutlier = "ids_no_proband_carriers_nonoutlier.txt" + File proband_ids = "proband_ids.txt" + } +} From 1efdd5d64c3b9bf49705f7df0a9f4c632299a5da Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:01:39 -0400 Subject: [PATCH 06/10] Add WDL workflow for creating Manta Tloc FOFN --- wdl/Create_Manta_Tloc_FOFN.wdl | 42 ++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 wdl/Create_Manta_Tloc_FOFN.wdl diff --git a/wdl/Create_Manta_Tloc_FOFN.wdl b/wdl/Create_Manta_Tloc_FOFN.wdl new file mode 100644 index 0000000..5b71a11 --- /dev/null +++ b/wdl/Create_Manta_Tloc_FOFN.wdl @@ -0,0 +1,42 @@ +version 1.0 + +workflow CreateMantaTlocFOFN { + input { + Array[String] manta_tloc_vcfs + String output_prefix + } + + call WriteFOFN { + input: + vcfs = manta_tloc_vcfs, + output_prefix = output_prefix + } + + output { + File manta_tloc_fofn = WriteFOFN.fofn + } +} + +task WriteFOFN { + input { + Array[String] vcfs + String output_prefix + } + + command <<< + set -euo pipefail + cat ~{write_lines(vcfs)} | sed 's|/mnt/disks/cromwell_root/|gs://|g' > ~{output_prefix}.manta_tloc_vcfs.fofn + >>> + + output { + File fofn = "~{output_prefix}.manta_tloc_vcfs.fofn" + } + + runtime { + cpu: 1 + memory: "1 GiB" + disks: "local-disk 10 HDD" + docker: "ubuntu:latest" + preemptible: 3 + } +} From dd33b0f9e4d2d886101fbfdfd4338bd56a52af75 Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:02:53 -0400 Subject: [PATCH 07/10] Add WDL workflow for counting HQ heterozygous SNP sites from small variant VCF This WDL script defines a workflow to count high-quality heterozygous SNV/Indel sites from VCF files, including input parameters and a task for processing the data. --- wdl/count_HQ_het_sites.wdl | 150 +++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) create mode 100644 wdl/count_HQ_het_sites.wdl diff --git a/wdl/count_HQ_het_sites.wdl b/wdl/count_HQ_het_sites.wdl new file mode 100644 index 0000000..65bc65b --- /dev/null +++ b/wdl/count_HQ_het_sites.wdl @@ -0,0 +1,150 @@ +version 1.0 + +workflow count_hq_het_sites { + + input { + File vcf + File vcf_index + Int min_gq = 20 + Int min_ad_total = 10 + Float min_ab = 0.2 + Float max_ab = 0.8 + String docker + Int mem_gb = 16 + Int cpu = 4 + } + + call CountHQHet { + input: + vcf = vcf, + vcf_index = vcf_index, + min_gq = min_gq, + min_ad_total = min_ad_total, + min_ab = min_ab, + max_ab = max_ab, + docker = docker, + mem_gb = mem_gb, + cpu = cpu + } + + output { + File hq_het_counts = CountHQHet.hq_het_counts + } + + meta { + author : "Generated WDL" + description : "Count HQ heterozygous SNV/Indel sites per chromosome per sample" + } +} + +task CountHQHet { + input { + File vcf + File vcf_index + Int min_gq + Int min_ad_total + Float min_ab + Float max_ab + String docker + Int mem_gb + Int cpu + Int disk_gb = ceil(size(vcf, "GiB") * 2) + 500 + } + + command <<< + set -euo pipefail + + # Debug: show what FILTER values actually look like to bcftools + echo "Sample FILTER values seen by bcftools:" >&2 + bcftools query -f '%FILTER\n' ~{vcf} | sort | uniq -c | sort -rn | head -10 >&2 + + # Keep only PASS sites (handles both VQSR-filtered VCFs with "PASS" + # and unfiltered GATK VCFs with ".") + bcftools view \ + --exclude-type other \ + --include 'FILTER="PASS" || FILTER="."' \ + --threads ~{cpu} \ + ~{vcf} \ + | bcftools query \ + -f '[%CHROM\t%SAMPLE\t%GT\t%GQ\t%AD\n]' \ + -o raw_genotypes.tsv + + echo "Total lines extracted: $(wc -l < raw_genotypes.tsv)" >&2 + echo "First 5 lines:" >&2 + head -5 raw_genotypes.tsv >&2 + + python3 -c " +import re, sys +from collections import defaultdict + +min_gq = ~{min_gq} +min_ad_total = ~{min_ad_total} +min_ab = ~{min_ab} +max_ab = ~{max_ab} + +HET_RE = re.compile(r'^(0[/|][1-9][0-9]*|[1-9][0-9]*[/|]0)$') + +counts = defaultdict(int) +skipped_gt = skipped_gq = skipped_ad = passed = 0 + +with open('raw_genotypes.tsv') as fh: + for line in fh: + parts = line.rstrip('\n').split('\t') + if len(parts) < 5: + continue + chrom, sample, gt, gq_s, ad_s = parts + + if not HET_RE.match(gt): + skipped_gt += 1 + continue + + try: + if float(gq_s) < min_gq: + skipped_gq += 1 + continue + except ValueError: + skipped_gq += 1 + continue + + try: + ad = ad_s.split(',') + ref_ad = float(ad[0]) + alt_ad = float(ad[1]) + total = ref_ad + alt_ad + if total < min_ad_total: + skipped_ad += 1 + continue + if not (min_ab <= alt_ad / total <= max_ab): + skipped_ad += 1 + continue + except (ValueError, IndexError): + skipped_ad += 1 + continue + + counts[(sample, chrom)] += 1 + passed += 1 + +print(f'Passed={passed} Skipped_GT={skipped_gt} Skipped_GQ={skipped_gq} Skipped_AD={skipped_ad}', file=sys.stderr) + +with open('hq_het_counts.tsv', 'w') as out: + out.write('sample\tchrom\tn_hq_het\n') + for (sample, chrom), n in sorted(counts.items()): + out.write(f'{sample}\t{chrom}\t{n}\n') + +print(f'Done. Wrote {len(counts)} sample-chrom combinations.', file=sys.stderr) +" + + >>> + + output { + File hq_het_counts = "hq_het_counts.tsv" + } + + runtime { + docker : docker + memory : "~{mem_gb} GiB" + cpu : cpu + disks : "local-disk ~{disk_gb} HDD" + preemptible : 2 + } +} From 06d7416b8c3a60c4a9423c04a00c8a476aa8ccb6 Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:03:59 -0400 Subject: [PATCH 08/10] Add WGS metrics collection workflow in WDL --- wdl/College_Multiple_Metrics_WGS.wdl | 351 +++++++++++++++++++++++++++ 1 file changed, 351 insertions(+) create mode 100644 wdl/College_Multiple_Metrics_WGS.wdl diff --git a/wdl/College_Multiple_Metrics_WGS.wdl b/wdl/College_Multiple_Metrics_WGS.wdl new file mode 100644 index 0000000..b024dcd --- /dev/null +++ b/wdl/College_Multiple_Metrics_WGS.wdl @@ -0,0 +1,351 @@ +workflow CollectMultipleMetricsWgs { + File input_bam + File input_bam_index + String base_name = basename(input_bam, ".bam") + + File ref_fasta + File ref_fasta_index + File ref_dict + File wgs_coverage_interval_list + + # average read length in the file: Picard's default 150bp + Int? read_length_input + Int read_length = select_first([read_length_input, 150]) + + # optional inputs for runtime parameters + Int? preemptible_attempts + Int? max_retries + Int? disk_pad + + # fixed runtime parameters + String? docker_override + String docker = select_first([docker_override, "us.gcr.io/broad-gotc-prod/genomes-in-the-cloud:2.3.2-1510681135"]) + Int preemptible_tries = select_first([preemptible_attempts, 3]) + Int max_retry_param = select_first([max_retries, 2]) + Int disk_size = ceil(bam_size + ref_size) + select_first([disk_pad, 10]) + + # disk size parameters + Float bam_size = size(input_bam, "GB") + size(input_bam_index, "GB") + Float ref_size = size(ref_fasta, "GB") + size(ref_fasta_index, "GB") + size(ref_dict, "GB") + + call CollectReadgroupBamQualityMetrics { + input: + input_bam = input_bam, + input_bam_index = input_bam_index, + output_prefix = base_name + ".readgroup", + ref_dict = ref_dict, + ref_fasta = ref_fasta, + ref_fasta_index = ref_fasta_index, + docker = docker, + preemptible_tries = preemptible_tries, + max_retries = max_retry_param, + disk_size = disk_size + } + + call CollectAggregationMetrics { + input: + input_bam = input_bam, + input_bam_index = input_bam_index, + output_prefix = base_name, + ref_dict = ref_dict, + ref_fasta = ref_fasta, + ref_fasta_index = ref_fasta_index, + docker = docker, + preemptible_tries = preemptible_tries, + max_retries = max_retry_param, + disk_size = disk_size + } + + call CollectWgsMetrics { + input: + input_bam = input_bam, + input_bam_index = input_bam_index, + output_prefix = base_name + ".wgs_metrics", + ref_dict = ref_dict, + ref_fasta = ref_fasta, + ref_fasta_index = ref_fasta_index, + wgs_coverage_interval_list = wgs_coverage_interval_list, + read_length = read_length, + docker = docker, + preemptible_tries = preemptible_tries, + max_retries = max_retry_param, + disk_size = disk_size + } + + call CollectRawWgsMetrics { + input: + input_bam = input_bam, + input_bam_index = input_bam_index, + output_prefix = base_name + ".raw_wgs_metrics", + ref_dict = ref_dict, + ref_fasta = ref_fasta, + ref_fasta_index = ref_fasta_index, + wgs_coverage_interval_list = wgs_coverage_interval_list, + read_length = read_length, + docker = docker, + preemptible_tries = preemptible_tries, + max_retries = max_retry_param, + disk_size = disk_size + } + + call ValidateSamFile { + input: + input_bam = input_bam, + input_bam_index = input_bam_index, + output_prefix = base_name + ".validation_report", + ref_dict = ref_dict, + ref_fasta = ref_fasta, + ref_fasta_index = ref_fasta_index, + docker = docker, + preemptible_tries = preemptible_tries, + max_retries = max_retry_param, + disk_size = disk_size + } + + output { + File readgroup_alignment_summary_metrics = CollectReadgroupBamQualityMetrics.alignment_summary_metrics + File readgroup_gc_bias_detail_metrics = CollectReadgroupBamQualityMetrics.gc_bias_detail_metrics + File readgroup_gc_bias_pdf = CollectReadgroupBamQualityMetrics.gc_bias_pdf + File readgroup_gc_bias_summary_metrics = CollectReadgroupBamQualityMetrics.gc_bias_summary_metrics + + File alignment_summary_metrics = CollectAggregationMetrics.alignment_summary_metrics + File bait_bias_detail_metrics = CollectAggregationMetrics.bait_bias_detail_metrics + File bait_bias_summary_metrics = CollectAggregationMetrics.bait_bias_summary_metrics + File gc_bias_detail_metrics = CollectAggregationMetrics.gc_bias_detail_metrics + File gc_bias_pdf = CollectAggregationMetrics.gc_bias_pdf + File gc_bias_summary_metrics = CollectAggregationMetrics.gc_bias_summary_metrics + File insert_size_histogram_pdf = CollectAggregationMetrics.insert_size_histogram_pdf + File insert_size_metrics = CollectAggregationMetrics.insert_size_metrics + File pre_adapter_detail_metrics = CollectAggregationMetrics.pre_adapter_detail_metrics + File pre_adapter_summary_metrics = CollectAggregationMetrics.pre_adapter_summary_metrics + File quality_distribution_pdf = CollectAggregationMetrics.quality_distribution_pdf + File quality_distribution_metrics = CollectAggregationMetrics.quality_distribution_metrics + + File wgs_metrics = CollectWgsMetrics.metrics + File raw_wgs_metrics = CollectRawWgsMetrics.metrics + File validate_sam_report = ValidateSamFile.report + } +} + + + +# Collect alignment summary and GC bias quality metrics +task CollectReadgroupBamQualityMetrics { + File input_bam + File input_bam_index + String output_prefix + File ref_dict + File ref_fasta + File ref_fasta_index + + String docker + Int preemptible_tries + Int max_retries + Int disk_size + + command { + java -Xms5000m -jar /usr/gitc/picard.jar \ + CollectMultipleMetrics \ + INPUT=${input_bam} \ + REFERENCE_SEQUENCE=${ref_fasta} \ + OUTPUT=${output_prefix} \ + ASSUME_SORTED=true \ + PROGRAM="null" \ + PROGRAM="CollectAlignmentSummaryMetrics" \ + PROGRAM="CollectGcBiasMetrics" \ + METRIC_ACCUMULATION_LEVEL="null" \ + METRIC_ACCUMULATION_LEVEL="READ_GROUP" + } + + runtime { + docker: docker + memory: "7 GB" + disks: "local-disk " + disk_size + " HDD" + preemptible: preemptible_tries + maxRetries: max_retries + } + output { + File alignment_summary_metrics = "${output_prefix}.alignment_summary_metrics" + File gc_bias_detail_metrics = "${output_prefix}.gc_bias.detail_metrics" + File gc_bias_pdf = "${output_prefix}.gc_bias.pdf" + File gc_bias_summary_metrics = "${output_prefix}.gc_bias.summary_metrics" + } +} + +# Collect quality metrics from the aggregated bam +task CollectAggregationMetrics { + File input_bam + File input_bam_index + String output_prefix + File ref_dict + File ref_fasta + File ref_fasta_index + + String docker + Int preemptible_tries + Int max_retries + Int disk_size + + command { + java -Xms5000m -jar /usr/gitc/picard.jar \ + CollectMultipleMetrics \ + INPUT=${input_bam} \ + REFERENCE_SEQUENCE=${ref_fasta} \ + OUTPUT=${output_prefix} \ + ASSUME_SORTED=true \ + PROGRAM="null" \ + PROGRAM="CollectAlignmentSummaryMetrics" \ + PROGRAM="CollectInsertSizeMetrics" \ + PROGRAM="CollectSequencingArtifactMetrics" \ + PROGRAM="CollectGcBiasMetrics" \ + PROGRAM="QualityScoreDistribution" \ + METRIC_ACCUMULATION_LEVEL="null" \ + METRIC_ACCUMULATION_LEVEL="SAMPLE" \ + METRIC_ACCUMULATION_LEVEL="LIBRARY" + + touch ${output_prefix}.insert_size_metrics + touch ${output_prefix}.insert_size_histogram.pdf + } + runtime { + docker: docker + memory: "7 GB" + disks: "local-disk " + disk_size + " HDD" + preemptible: preemptible_tries + maxRetries: max_retries + } + output { + File alignment_summary_metrics = "${output_prefix}.alignment_summary_metrics" + File bait_bias_detail_metrics = "${output_prefix}.bait_bias_detail_metrics" + File bait_bias_summary_metrics = "${output_prefix}.bait_bias_summary_metrics" + File gc_bias_detail_metrics = "${output_prefix}.gc_bias.detail_metrics" + File gc_bias_pdf = "${output_prefix}.gc_bias.pdf" + File gc_bias_summary_metrics = "${output_prefix}.gc_bias.summary_metrics" + File insert_size_histogram_pdf = "${output_prefix}.insert_size_histogram.pdf" + File insert_size_metrics = "${output_prefix}.insert_size_metrics" + File pre_adapter_detail_metrics = "${output_prefix}.pre_adapter_detail_metrics" + File pre_adapter_summary_metrics = "${output_prefix}.pre_adapter_summary_metrics" + File quality_distribution_pdf = "${output_prefix}.quality_distribution.pdf" + File quality_distribution_metrics = "${output_prefix}.quality_distribution_metrics" + } +} + +# Note these tasks will break if the read lengths in the bam are greater than 250. +task CollectWgsMetrics { + File input_bam + File input_bam_index + String output_prefix + File ref_dict + File ref_fasta + File ref_fasta_index + + File wgs_coverage_interval_list + Int read_length + + String docker + Int preemptible_tries + Int max_retries + Int disk_size + + command { + java -Xms2000m -jar /usr/gitc/picard.jar \ + CollectWgsMetrics \ + INPUT=${input_bam} \ + VALIDATION_STRINGENCY=SILENT \ + REFERENCE_SEQUENCE=${ref_fasta} \ + INCLUDE_BQ_HISTOGRAM=true \ + INTERVALS=${wgs_coverage_interval_list} \ + OUTPUT=${output_prefix} \ + USE_FAST_ALGORITHM=true \ + READ_LENGTH=${read_length} + } + runtime { + docker: docker + preemptible: preemptible_tries + memory: "3 GB" + disks: "local-disk " + disk_size + " HDD" + maxRetries: max_retries + } + output { + File metrics = "${output_prefix}" + } +} + +# Collect raw WGS metrics (commonly used QC thresholds) +task CollectRawWgsMetrics { + File input_bam + File input_bam_index + String output_prefix + File ref_dict + File ref_fasta + File ref_fasta_index + + File wgs_coverage_interval_list + Int read_length + + String docker + Int preemptible_tries + Int max_retries + Int disk_size + + command { + java -Xms2000m -jar /usr/gitc/picard.jar \ + CollectRawWgsMetrics \ + INPUT=${input_bam} \ + VALIDATION_STRINGENCY=SILENT \ + REFERENCE_SEQUENCE=${ref_fasta} \ + INCLUDE_BQ_HISTOGRAM=true \ + INTERVALS=${wgs_coverage_interval_list} \ + OUTPUT=${output_prefix} \ + USE_FAST_ALGORITHM=true \ + READ_LENGTH=${read_length} + } + runtime { + docker: docker + preemptible: preemptible_tries + memory: "30 GB" + disks: "local-disk " + disk_size + " HDD" + maxRetries: max_retries + } + output { + File metrics = "${output_prefix}" + } +} + + +task ValidateSamFile { + File input_bam + File input_bam_index + String output_prefix + File ref_dict + File ref_fasta + File ref_fasta_index + + String docker + Int preemptible_tries + Int max_retries + Int disk_size + Int mem_multiplier = 1 + Int mem_gb = 8 * mem_multiplier + + command { + java -jar /usr/gitc/picard.jar \ + ValidateSamFile \ + INPUT=${input_bam} \ + OUTPUT=${output_prefix} \ + REFERENCE_SEQUENCE=${ref_fasta} \ + IGNORE="MISSING_TAG_NM" \ + IGNORE="INVALID_VERSION_NUMBER" \ + MODE=VERBOSE \ + IS_BISULFITE_SEQUENCED=false + } + runtime { + docker: docker + preemptible: preemptible_tries + memory: "${mem_gb} GB" + disks: "local-disk " + disk_size + " HDD" + maxRetries: max_retries + } + output { + File report = "${output_prefix}" + } +} From 754d3bb6a3a5be170f52f4ee5886f4ff90a2ca16 Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:05:16 -0400 Subject: [PATCH 09/10] Additional QC metrics for WGS samples --- ...Collect_Dup_Contamination_gVCF_metrics.wdl | 261 ++++++++++++++++++ 1 file changed, 261 insertions(+) create mode 100644 wdl/Collect_Dup_Contamination_gVCF_metrics.wdl diff --git a/wdl/Collect_Dup_Contamination_gVCF_metrics.wdl b/wdl/Collect_Dup_Contamination_gVCF_metrics.wdl new file mode 100644 index 0000000..82e1d91 --- /dev/null +++ b/wdl/Collect_Dup_Contamination_gVCF_metrics.wdl @@ -0,0 +1,261 @@ +version 1.0 + +## Portions Copyright Broad Institute, 2018 +## +## This WDL pipeline implements QC in human whole-genome or exome/targeted sequencing data. +## +## Requirements/expectations +## - Human paired-end sequencing data in aligned BAM or CRAM format +## - Input BAM/CRAM files must additionally comply with the following requirements: +## - - files must pass validation by ValidateSamFile +## - - reads are provided in query-sorted order +## - - all reads must have an RG tag +## - Reference genome must be Hg38 with ALT contigs +## +## Runtime parameters are optimized for Broad's Google Cloud Platform implementation. +## For program versions, see docker containers. +## +## LICENSING : +## This script is released under the WDL open source code license (BSD-3). +## Full license text at https://github.com/openwdl/wdl/blob/master/LICENSE +## Note however that the programs it calls may be subject to different licenses. +## Users are responsible for checking that they are authorized to run all programs before running this script. +## - [Picard](https://broadinstitute.github.io/picard/) +## - [VerifyBamID2](https://github.com/Griffan/VerifyBamID) + +# Git URL import +#import "tasks/Qc.wdl" as QC + +# WORKFLOW DEFINITION +workflow SingleSampleQc { + input { + File input_bam + File input_bam_index + File input_gvcf + File input_gvcf_index + File dbsnp_vcf + File dbsnp_vcf_index + File ref_cache + File ref_dict + File ref_fasta + File ref_fasta_index + String base_name + Int preemptible_tries + File coverage_interval_list + File contamination_sites_ud + File contamination_sites_bed + File contamination_sites_mu + File evaluation_interval_list + Boolean is_wgs + Boolean? is_outlier_data + Boolean is_gvcf = true + + File evaluation_thresholds + } + + # Not overridable: + Int read_length = 250 + + # Generate a BAM or CRAM index + # Estimate level of cross-sample contamination + call CheckContamination { + input: + input_bam = input_bam, + input_bam_index = input_bam_index, + contamination_sites_ud = contamination_sites_ud, + contamination_sites_bed = contamination_sites_bed, + contamination_sites_mu = contamination_sites_mu, + ref_fasta = ref_fasta, + ref_fasta_index = ref_fasta_index, + output_prefix = base_name + ".verify_bam_id", + preemptible_tries = preemptible_tries, + } + + # Calculate the duplication rate since MarkDuplicates was already performed + call CollectDuplicateMetrics { + input: + input_bam = input_bam, + input_bam_index = input_bam_index, + output_bam_prefix = base_name, + ref_dict = ref_dict, + ref_fasta = ref_fasta, + ref_fasta_index = ref_fasta_index, + preemptible_tries = preemptible_tries + } + +call CollectVariantCallingMetrics { + input: + input_vcf = input_gvcf, + input_vcf_index = input_gvcf_index, + metrics_basename = base_name, + dbsnp_vcf = dbsnp_vcf, + dbsnp_vcf_index = dbsnp_vcf_index, + ref_dict = ref_dict, + evaluation_interval_list = evaluation_interval_list, + is_gvcf = is_gvcf, + preemptible_tries = preemptible_tries + } + + # Outputs that will be retained when execution is complete + output { + + + File selfSM = CheckContamination.selfSM + Float contamination = CheckContamination.contamination + + File duplication_metrics_file = CollectDuplicateMetrics.duplication_metrics_file + String percent_duplication = CollectDuplicateMetrics.percent_duplication + + File vcf_summary_metrics = CollectVariantCallingMetrics.summary_metrics + File vcf_detail_metrics = CollectVariantCallingMetrics.detail_metrics + + } +} +task CheckContamination { + input { + File input_bam + File input_bam_index + File contamination_sites_ud + File contamination_sites_bed + File contamination_sites_mu + File ref_fasta + File ref_fasta_index + String output_prefix + Int preemptible_tries + Boolean disable_sanity_check = false + } + + Int disk_size = ceil(size(input_bam, "GiB") + size(ref_fasta, "GiB")) + 30 + + command <<< + set -e + + # creates a ~{output_prefix}.selfSM file, a TSV file with 2 rows, 19 columns. + # First row are the keys (e.g., SEQ_SM, RG, FREEMIX), second row are the associated values + /usr/gitc/VerifyBamID \ + --Verbose \ + --NumPC 4 \ + --Output ~{output_prefix} \ + --BamFile ~{input_bam} \ + --Reference ~{ref_fasta} \ + --UDPath ~{contamination_sites_ud} \ + --MeanPath ~{contamination_sites_mu} \ + --BedPath ~{contamination_sites_bed} \ + ~{true="--DisableSanityCheck" false="" disable_sanity_check} \ + 1>/dev/null + + # used to read from the selfSM file and calculate contamination, which gets printed out + python3 <>> + runtime { + preemptible: preemptible_tries + memory: "4 GiB" + disks: "local-disk " + disk_size + " HDD" + docker: "us.gcr.io/broad-gotc-prod/verify-bam-id:c1cba76e979904eb69c31520a0d7f5be63c72253-1553018888" + cpu: "2" + } + output { + File selfSM = "~{output_prefix}.selfSM" + Float contamination = read_float(stdout()) + Map[String, String] metrics = { "FREEMIX": read_string(stdout()) } + } +} + +task CollectVariantCallingMetrics { + input { + File input_vcf + File input_vcf_index + String metrics_basename + File dbsnp_vcf + File dbsnp_vcf_index + File ref_dict + File evaluation_interval_list + Boolean is_gvcf = true + Int preemptible_tries + } + + Int disk_size = ceil(size(input_vcf, "GiB") + size(dbsnp_vcf, "GiB")) + 20 + + command { + java -Xms2000m -Xmx2500m -jar /usr/picard/picard.jar \ + CollectVariantCallingMetrics \ + INPUT=~{input_vcf} \ + OUTPUT=~{metrics_basename} \ + DBSNP=~{dbsnp_vcf} \ + SEQUENCE_DICTIONARY=~{ref_dict} \ + TARGET_INTERVALS=~{evaluation_interval_list} \ + ~{true="GVCF_INPUT=true" false="" is_gvcf} + } + runtime { + docker: "us.gcr.io/broad-gotc-prod/picard-cloud:2.26.10" + preemptible: preemptible_tries + memory: "3000 MiB" + disks: "local-disk " + disk_size + " HDD" + } + output { + File summary_metrics = "~{metrics_basename}.variant_calling_summary_metrics" + File detail_metrics = "~{metrics_basename}.variant_calling_detail_metrics" + } +} + +task CollectDuplicateMetrics { + input { + File input_bam + File input_bam_index + String output_bam_prefix + File ref_dict + File ref_fasta + File ref_fasta_index + Int preemptible_tries + } + + Float ref_size = size(ref_fasta, "GiB") + size(ref_fasta_index, "GiB") + size(ref_dict, "GiB") + Int disk_size = ceil(size(input_bam, "GiB") + ref_size) + 20 + + String duplication_metric_object_file = '~{output_bam_prefix}.duplication_metrics.metrics_only' + + command <<< + java -Xms5000m -jar /usr/picard/picard.jar \ + CollectDuplicateMetrics \ + METRICS_FILE=~{output_bam_prefix}.duplication_metrics \ + INPUT=~{input_bam} \ + ASSUME_SORTED=true \ + REFERENCE_SEQUENCE=~{ref_fasta} + + grep -v '#' '~{output_bam_prefix}.duplication_metrics' | grep '.\+' | perl -E 'my ($keys, $values) = <>; chomp $keys; chomp $values; my @k = split("\t", $keys); my @v = split("\t", $values); for(0..$#k) { say join("\t", $k[$_], $v[$_]); }' > ~{duplication_metric_object_file} + + >>> + runtime { + docker: "us.gcr.io/broad-gotc-prod/picard-cloud:2.21.7" + memory: "7 GiB" + disks: "local-disk " + disk_size + " HDD" + preemptible: preemptible_tries + } + output { + File duplication_metrics_file = "~{output_bam_prefix}.duplication_metrics" + Map[String, String] duplication_metrics = read_map(duplication_metric_object_file) + String percent_duplication = duplication_metrics["PERCENT_DUPLICATION"] + } +} + + From e062313cd7dcdf181f4ac42b4e032b6d81c80cfd Mon Sep 17 00:00:00 2001 From: knathamuni <144271815+knathamuni@users.noreply.github.com> Date: Tue, 23 Jun 2026 17:06:04 -0400 Subject: [PATCH 10/10] Add WDL workflow for converting paired FASTQ to BAM --- wdl/Paired_Fastq_to_ubam.wdl | 205 +++++++++++++++++++++++++++++++++++ 1 file changed, 205 insertions(+) create mode 100644 wdl/Paired_Fastq_to_ubam.wdl diff --git a/wdl/Paired_Fastq_to_ubam.wdl b/wdl/Paired_Fastq_to_ubam.wdl new file mode 100644 index 0000000..155465e --- /dev/null +++ b/wdl/Paired_Fastq_to_ubam.wdl @@ -0,0 +1,205 @@ +version 1.0 + +workflow ConvertPairedFastQsToUnmappedBamWf { + + input { + String sample_name + File fastq_1 + File fastq_2 + Array[String] readgroup_names + Array[String] library_names + Array[String] platform_units + Array[String] run_dates + Array[String] platform_names + Array[String] sequencing_centers + + File ref_dict + + String samtools_docker = "staphb/samtools:1.17" + String gatk_docker = "broadinstitute/gatk:4.5.0.0" + + Int additional_disk_space_gb = 50 + Int machine_mem_gb = 7 + Int preemptible_attempts = 3 + } + + call PairedFastQsToUnmappedBAM { + input: + sample_name = sample_name, + fastq_1 = fastq_1, + fastq_2 = fastq_2, + readgroup_names = readgroup_names, + library_names = library_names, + platform_units = platform_units, + run_dates = run_dates, + platform_names = platform_names, + sequencing_centers = sequencing_centers, + ref_dict = ref_dict, + docker = samtools_docker, + additional_disk_space_gb = additional_disk_space_gb, + machine_mem_gb = machine_mem_gb, + preemptible_attempts = preemptible_attempts + } + + call SortSam { + input: + input_bam = PairedFastQsToUnmappedBAM.output_unmapped_bam, + sample_name = sample_name, + docker = gatk_docker, + machine_mem_gb = machine_mem_gb, + additional_disk_space_gb = additional_disk_space_gb, + preemptible_attempts = preemptible_attempts + } + + output { + File output_unmapped_bam = SortSam.sorted_bam + } +} + +task PairedFastQsToUnmappedBAM { + input { + String sample_name + File fastq_1 + File fastq_2 + Array[String] readgroup_names + Array[String] library_names + Array[String] platform_units + Array[String] run_dates + Array[String] platform_names + Array[String] sequencing_centers + + File ref_dict + + Int additional_disk_space_gb = 50 + Int machine_mem_gb = 7 + Int preemptible_attempts = 3 + String docker + } + + Int disk_space_gb = ceil((size(fastq_1, "GB") + size(fastq_2, "GB")) * 4) + additional_disk_space_gb + + command <<< + set -euo pipefail + + # Read arrays from files to avoid shell quoting issues + mapfile -t RG_NAMES < ~{write_lines(readgroup_names)} + mapfile -t LIB_NAMES < ~{write_lines(library_names)} + mapfile -t PLAT_UNITS < ~{write_lines(platform_units)} + mapfile -t RUN_DATES < ~{write_lines(run_dates)} + mapfile -t PLAT_NAMES < ~{write_lines(platform_names)} + mapfile -t SEQ_CENTERS < ~{write_lines(sequencing_centers)} + + n_rg=${#RG_NAMES[@]} + + # Validate all arrays are the same length + for arr_name in LIB_NAMES PLAT_UNITS RUN_DATES PLAT_NAMES SEQ_CENTERS; do + declare -n arr="$arr_name" + n=${#arr[@]} + if [ "$n" -ne "$n_rg" ]; then + echo "ERROR: $arr_name length ($n) does not match readgroup_names length ($n_rg)" >&2 + exit 1 + fi + done + + # Build header with @SQ lines from ref_dict. + # Picard 2.26.x MergeBamAlignment requires the uBAM sequence dictionary + # to match the reference — an empty dictionary causes a merge failure. + printf "@HD\tVN:1.6\tSO:queryname\n" > new_header.sam + grep "^@SQ" ~{ref_dict} >> new_header.sam + + for (( i=0; i> new_header.sam + done + + echo "Header to be applied:" + cat new_header.sam + echo "RG lines: $(grep -c "^@RG" new_header.sam)" + echo "SQ lines: $(grep -c "^@SQ" new_header.sam)" + + FIRST_RG=$(printf "ID:%s\tSM:%s\tLB:%s\tPU:%s\tDT:%s\tPL:%s\tCN:%s" \ + "${RG_NAMES[0]}" \ + "~{sample_name}" \ + "${LIB_NAMES[0]}" \ + "${PLAT_UNITS[0]}" \ + "${RUN_DATES[0]}" \ + "${PLAT_NAMES[0]}" \ + "${SEQ_CENTERS[0]}") + + samtools import \ + -1 ~{fastq_1} \ + -2 ~{fastq_2} \ + -r "$FIRST_RG" \ + -O BAM \ + -o raw.unmapped.bam + + samtools reheader new_header.sam raw.unmapped.bam > ~{sample_name}.unmapped.bam + + echo "Done. Verifying output BAM header:" + samtools view -H ~{sample_name}.unmapped.bam | grep "^@RG" + echo "SQ lines: $(samtools view -H ~{sample_name}.unmapped.bam | grep -c "^@SQ")" + echo "Flagstat:" + samtools flagstat ~{sample_name}.unmapped.bam + >>> + + output { + File output_unmapped_bam = "~{sample_name}.unmapped.bam" + } + + runtime { + docker: docker + memory: machine_mem_gb + " GB" + disks: "local-disk " + disk_space_gb + " HDD" + preemptible: preemptible_attempts + maxRetries: preemptible_attempts + } +} + +task SortSam { + input { + File input_bam + String sample_name + String docker + Int machine_mem_gb = 7 + Int additional_disk_space_gb = 50 + Int preemptible_attempts = 3 + } + + Int command_mem_gb = machine_mem_gb - 1 + Int disk_space_gb = ceil(size(input_bam, "GB") * 4) + additional_disk_space_gb + + command <<< + set -euo pipefail + + mkdir -p tmp + + gatk --java-options "-Xmx~{command_mem_gb}g -Djava.io.tmpdir=./tmp" \ + SortSam \ + -I ~{input_bam} \ + -O ~{sample_name}.query_sorted.unmapped.bam \ + --SORT_ORDER queryname \ + --TMP_DIR ./tmp + + echo "Sort complete. Flagstat:" + samtools flagstat ~{sample_name}.query_sorted.unmapped.bam + >>> + + output { + File sorted_bam = "~{sample_name}.query_sorted.unmapped.bam" + } + + runtime { + docker: docker + memory: machine_mem_gb + " GB" + disks: "local-disk " + disk_space_gb + " HDD" + preemptible: preemptible_attempts + maxRetries: preemptible_attempts + } +}