Run cnmops on chromosome arms separately (#925)

Author: VJalili

src/sv-pipeline/scripts/cytoband_to_arms.py

@@ -0,0 +1,59 @@
+import gzip
+import argparse
+import re
+from collections import defaultdict
+from dataclasses import dataclass
+
+
+@dataclass
+class ArmData:
+    p_end: int = 0
+    q_start: int = None
+    q_end: int = 0
+
+
+def natural_sort_key(s):
+    return [int(text) if text.isdigit() else text.lower() for text in re.split("([0-9]+)", s)]
+
+
+def main(input_file, output_file):
+    arms = defaultdict(ArmData)
+
+    with gzip.open(input_file, "rt") as f:
+        for line in f:
+            if line.startswith("#") or not line.strip():
+                continue
+
+            parts = line.strip().split("\t")
+            chrom, start, end, band = parts[0], int(parts[1]), int(parts[2]), parts[3]
+            arm_type = band[0].lower()
+
+            if arm_type == "p":
+                if end > arms[chrom].p_end:
+                    arms[chrom].p_end = end
+            elif arm_type == "q":
+                # Convert 0-based BED start to 1-based GATK start
+                one_based_start = start + 1
+                if arms[chrom].q_start is None:
+                    arms[chrom].q_start = one_based_start
+                arms[chrom].q_end = end
+
+    with open(output_file, "w") as out:
+        for chrom in sorted(arms.keys(), key=natural_sort_key):
+            a = arms[chrom]
+
+            # P-arm starts at 1, Q-arm starts at the first q-band + 1
+            p_range = f"1-{a.p_end}" if a.p_end > 0 else "0"
+            q_range = f"{a.q_start}-{a.q_end}" if a.q_start is not None else "0"
+
+            # Output format: chrom \t p-interval \t q-interval
+            out.write(f"{chrom}\t{p_range}\t{q_range}\n")
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-i", "--input", required=True)
+    parser.add_argument("-o", "--output", required=True)
+    args = parser.parse_args()
+
+    main(args.input, args.output)

src/sv_shell/cnmops.sh

@@ -19,21 +19,23 @@ JVM_MAX_MEM=$(getJavaMem MemTotal)
 echo "JVM memory: $JVM_MAX_MEM"
 
 
-function CNSampleNormal() {
+function CNSampleNormalArm() {
   local _chr=$1
   local _mode=$2
   local _r=$3
+  local _interval=${4:-${_chr}}
 
   echo "----------- Starting CN Sample Normal -------------"
   echo "chr: ${_chr}"
+  echo "interval: ${_interval}"
   echo "mode: ${_mode}"
   echo "r: ${_r}"
   echo "---------------------------------------------------"
 
   java "-Xmx${JVM_MAX_MEM}" -jar /opt/gatk.jar PrintSVEvidence \
     --sequence-dictionary "${ref_dict}" \
     --evidence-file "${bincov_matrix}" \
-    -L "${_chr}" \
+    -L "${_interval}" \
     -O "${_chr}.RD.txt"
 
   if [ "${_mode}" == "normal" ]; then
@@ -50,23 +52,55 @@ function CNSampleNormal() {
   EMPTY_OUTPUT_ERROR="No CNV regions in result object. Rerun cn.mops with different parameters!"
   set +e
   echo "Starting to run cnMOPS_workflow"
-  bash /opt/WGD/bin/cnMOPS_workflow.sh -S "${exclude_list}" -x "${exclude_list}" -r "${_r}" -o . -M "${_chr}.${_mode}.RD.txt" </dev/null wor2>&1 | tee cnmops.out
+  bash /opt/WGD/bin/cnMOPS_workflow.sh -S "${exclude_list}" -x "${exclude_list}" -r "${_r}" -o . -M "${_chr}.${_mode}.RD.txt" </dev/null 2>&1 | tee cnmops.out
   echo "Finished running cnMOPS_workflow"
   RC=$?
   set -e
-  if [ ! $RC -eq 0 ]; then
-    if grep -q "$EMPTY_OUTPUT_ERROR" "cnmops.out"; then
-      touch calls/cnMOPS.cnMOPS.gff
-    else
-      echo "cnMOPS_workflow.sh returned a non-zero code that was not due to an empty call file."
-      exit $RC
-    fi
+  if grep -q "$EMPTY_OUTPUT_ERROR" "cnmops.out"; then
+    echo "No CNV regions detected, creating empty GFF."
+    touch calls/cnMOPS.cnMOPS.gff
+  elif [ ! $RC -eq 0 ]; then
+    echo "cnMOPS_workflow.sh returned a non-zero code that was not due to an empty call file."
+    exit $RC
+  elif [ ! -f calls/cnMOPS.cnMOPS.gff ]; then
+    echo "cnMOPS_workflow.sh succeeded but did not produce calls/cnMOPS.cnMOPS.gff. Creating empty file."
+    touch calls/cnMOPS.cnMOPS.gff
   fi
 
   echo "----------- Finished CN Sample Normal -------------"
 }
 
 
+function CNSampleNormal() {
+  local _chr=$1
+  local _mode=$2
+  local _r=$3
+  local _parent_dir
+  _parent_dir=$(pwd)
+  local _arm_gffs=()
+
+  local _arm_idx=0
+  for _arm_range in $(awk -v chr="${_chr}" '$1 == chr {for(i=2;i<=NF;i++) print $i}' "${chr_arms_file}"); do
+    _arm_idx=$((_arm_idx + 1))
+    local _arm_dir="${_parent_dir}/arm_${_arm_idx}"
+    mkdir -p "${_arm_dir}"
+    cd "${_arm_dir}"
+
+    CNSampleNormalArm "${_chr}_arm${_arm_idx}" "${_mode}" "${_r}" "${_chr}:${_arm_range}"
+    _arm_gffs+=("${_arm_dir}/calls/cnMOPS.cnMOPS.gff")
+
+    cd "${_parent_dir}"
+  done
+
+  mkdir -p calls
+  if [ ${#_arm_gffs[@]} -gt 0 ]; then
+    cat "${_arm_gffs[@]}" > calls/cnMOPS.cnMOPS.gff
+  else
+    touch calls/cnMOPS.cnMOPS.gff
+  fi
+}
+
+
 function CleanCNMops() {
   local _sample_list=$1
   local _cnmops_gff=$2
@@ -149,25 +183,32 @@ working_dir="$(realpath ${working_dir})"
 cd "${working_dir}"
 echo "cnMOPS Working directory: ${working_dir}"
 
-batch=($(jq -r '.batch' "$input_json"))
-allo_file=($(jq -r '.allo_file' "$input_json"))
-chrom_file=($(jq -r '.chrom_file' "$input_json"))
-exclude_list=($(jq -r '.exclude_list' "$input_json"))
-ped_file=($(jq -r '.ped_file' "$input_json"))
-r1=($(jq -r '.r1' "$input_json"))
-r2=($(jq -r '.r2' "$input_json"))
-ref_dict=($(jq -r '.ref_dict' "$input_json"))
-bincov_matrix=($(jq -r '.bincov_matrix' "$input_json"))
-stitch_and_clean_large_events=($(jq -r '.stitch_and_clean_large_events' "$input_json"))
-min_size=($(jq -r '.min_size' "$input_json"))
-prefix=($(jq -r '.prefix' "$input_json"))
+batch=$(jq -r '.batch' "$input_json")
+allo_file=$(jq -r '.allo_file' "$input_json")
+chrom_file=$(jq -r '.chrom_file' "$input_json")
+exclude_list=$(jq -r '.exclude_list' "$input_json")
+ped_file=$(jq -r '.ped_file' "$input_json")
+r1=$(jq -r '.r1' "$input_json")
+r2=$(jq -r '.r2' "$input_json")
+ref_dict=$(jq -r '.ref_dict' "$input_json")
+bincov_matrix=$(jq -r '.bincov_matrix' "$input_json")
+stitch_and_clean_large_events=$(jq -r '.stitch_and_clean_large_events' "$input_json")
+min_size=$(jq -r '.min_size' "$input_json")
+prefix=$(jq -r '.prefix' "$input_json")
+cytobands=$(jq -r '.cytobands' "$input_json")
 
 
 # -------------------------------------------------------
 # ======================= Command =======================
 # -------------------------------------------------------
 
 
+# Get chr arms from cytobands
+# ---------------------------
+chr_arms_file=$(realpath "chr_arms.bed")
+python /opt/sv-pipeline/scripts/cytoband_to_arms.py -i "${cytobands}" -o "${chr_arms_file}"
+
+
 male_r1_gff=()
 male_r2_gff=()
 allos=($(awk '{print $1}' "${allo_file}"))

src/sv_shell/gather_batch_evidence.sh

@@ -186,6 +186,7 @@ cnmops_inputs_json="$(realpath "${output_dir}/cnmops_inputs.json")"
 cnmops_outputs_json="$(realpath "${output_dir}/cnmops_outputs.json")"
 
 jq -n \
+  --slurpfile inputs "${input_json}" \
   --arg r1 "3" \
   --arg r2 "10" \
   --arg batch "${batch}" \
@@ -213,7 +214,8 @@ jq -n \
       "allo_file": $allo_file,
       "ref_dict": $ref_dict,
       "prefix": $prefix,
-      "stitch_and_clean_large_events": $stitch_and_clean_large_events
+      "stitch_and_clean_large_events": $stitch_and_clean_large_events,
+      "cytobands": $inputs[0].cytobands
   }' > "${cnmops_inputs_json}"
 
 bash /opt/sv_shell/cnmops.sh "${cnmops_inputs_json}" "${cnmops_outputs_json}"
@@ -229,6 +231,7 @@ cnmops_large_inputs_json="$(realpath "${output_dir}/cnmops_large_inputs.json")"
 cnmops_large_outputs_json="$(realpath "${output_dir}/cnmops_large_outputs.json")"
 
 jq -n \
+  --slurpfile inputs "${input_json}" \
   --arg r1 "1000" \
   --arg r2 "100" \
   --arg batch "${batch}" \
@@ -256,7 +259,8 @@ jq -n \
       "allo_file": $allo_file,
       "ref_dict": $ref_dict,
       "prefix": $prefix,
-      "stitch_and_clean_large_events": $stitch_and_clean_large_events
+      "stitch_and_clean_large_events": $stitch_and_clean_large_events,
+      "cytobands": $inputs[0].cytobands
   }' > "${cnmops_large_inputs_json}"
 
 bash /opt/sv_shell/cnmops.sh "${cnmops_large_inputs_json}" "${cnmops_large_outputs_json}"

src/sv_shell/sample_inputs/cnmops.json

@@ -2,14 +2,15 @@
   "batch": "NA12878",
   "r1": "3",
   "r2": "10",
-  "bincov_matrix": "/inputs/NA12878.RD.txt.gz",
+  "bincov_matrix": "/intermediate/NA12878.RD.txt.gz",
   "chrom_file": "/inputs/autosome.fai",
-  "ped_file": "/inputs/combined_ped_file.ped",
+  "ped_file": "/intermediate/combined_ped_file.ped",
   "exclude_list": "/inputs/GRCh38_Nmask.bed",
   "allo_file": "/inputs/allosome.fai",
   "ref_dict": "/inputs/Homo_sapiens_assembly38.dict",
   "samples": [ "NA12878", "HG00096", "HG00129", "HG00140", "HG00150", "HG00187", "HG00239", "HG00277", "HG00288", "HG00337", "HG00349", "HG00375", "HG00410", "HG00457", "HG00557", "HG00599", "HG00625", "HG00701", "HG00740", "HG00844", "HG01060", "HG01085", "HG01112", "HG01275", "HG01325", "HG01344", "HG01356", "HG01384", "HG01393", "HG01396", "HG01474", "HG01507", "HG01572", "HG01607", "HG01709", "HG01747", "HG01790", "HG01794", "HG01799", "HG01861", "HG01874", "HG01880", "HG01885", "HG01958", "HG01982", "HG02002", "HG02010", "HG02019", "HG02020", "HG02069", "HG02085", "HG02186", "HG02221", "HG02235", "HG02272", "HG02275", "HG02299", "HG02332", "HG02367", "HG02374", "HG02489", "HG02490", "HG02491", "HG02586", "HG02588", "HG02611", "HG02620", "HG02642", "HG02648", "HG02658", "HG02855", "HG02953", "HG03007", "HG03009", "HG03085", "HG03099", "HG03100", "HG03111", "HG03369", "HG03370", "HG03436", "HG03449", "HG03472", "HG03476", "HG03556", "HG03604", "HG03649", "HG03684", "HG03694", "HG03709", "HG03722", "HG03727", "HG03744", "HG03756", "HG03789", "HG03850", "HG03864", "HG03872", "HG03888", "HG04118", "HG04158", "HG04161", "HG04183", "NA06984", "NA10847", "NA11894", "NA12340", "NA12489", "NA12872", "NA18499", "NA18507", "NA18530", "NA18539", "NA18549", "NA18553", "NA18560", "NA18638", "NA18923", "NA18941", "NA18945", "NA18956", "NA18995", "NA19001", "NA19035", "NA19062", "NA19102", "NA19143", "NA19184", "NA19350", "NA19351", "NA19377", "NA19443", "NA19449", "NA19661", "NA19678", "NA19679", "NA19684", "NA19746", "NA19795", "NA19818", "NA19913", "NA20126", "NA20320", "NA20321", "NA20346", "NA20509", "NA20510", "NA20522", "NA20752", "NA20764", "NA20802", "NA20845", "NA20869", "NA20895", "NA21102", "NA21122", "NA21133" ],
   "prefix": "header",
   "stitch_and_clean_large_events": false,
-  "min_size": 1000000
+  "min_size": 1000000,
+  "cytobands": "/inputs/cytobands_hg38.bed.gz"
 }

src/sv_shell/sample_inputs/gather_batch_evidence.json

@@ -520,5 +520,6 @@
   "wham_vcfs": [
     "/inputs/NA12878.wham.vcf.gz"
   ],
-  "min_svsize": 50
+  "min_svsize": 50,
+  "cytobands": "/inputs/cytobands_hg38.bed.gz"
 }

src/sv_shell/single_sample_pipeline.sh

@@ -319,7 +319,8 @@ jq -n \
       run_module_metrics: $inputs[0].run_batchevidence_metrics,
       median_cov_mem_gb_per_sample: $inputs[0].median_cov_mem_gb_per_sample,
       ref_panel_median_cov: $inputs[0].ref_panel_median_cov,
-      sample_median_cov: $eqc_outputs[0].bincov_median
+      sample_median_cov: $eqc_outputs[0].bincov_median,
+      "cytobands": $inputs[0].cytobands,
   }' > "${gather_batch_evidence_inputs_json_filename}"
 
 bash /opt/sv_shell/gather_batch_evidence.sh \