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"] + } +} + + 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}" + } +} 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 + } +} 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" + } +} 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 + } +} 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 + } +} 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" + } +} 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 + } +} 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 + } +} 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 + } +}