add mafft_and_snp_annotated.wdl (#194)

* add mafft_and_snp_annotated.wdl

* rename workflow

* Update .dockstore.yml

* replace s/^1/$chr_name/g in annotate_vcf_snpeff

The vcf output by snp-sites usus "1" as the chromosome name for its output, while snpEff expects the chromosome name to match that of the reference. This renames the chr name in annotate_vcf_snpeff so snpEff succeeds. This will fail for multi-chr vcfs.

* remove txt output from merge_vcfs_and_annotate

Author: tomkinsc

.dockstore.yml

@@ -180,6 +180,11 @@ workflows:
    primaryDescriptorPath: /pipes/WDL/workflows/mafft_and_snp.wdl
    testParameterFiles:
     - empty.json
+ - name: mafft_and_snp_annotated
+   subclass: WDL
+   primaryDescriptorPath: /pipes/WDL/workflows/mafft_and_snp_annotated.wdl
+   testParameterFiles:
+    - empty.json
  - name: mafft_and_trim
    subclass: WDL
    primaryDescriptorPath: /pipes/WDL/workflows/mafft_and_trim.wdl

pipes/WDL/tasks/tasks_intrahost.wdl

@@ -150,32 +150,53 @@ task annotate_vcf_snpeff {
     fi
     echo "snpRefAccessions: $snpRefAccessions"
 
+    vcf_to_use=""
     if (file "${in_vcf}" | grep -q "gzip" ) ; then
       echo "${in_vcf} is already compressed"
+      vcf_to_use="${in_vcf}"
     else
       echo "${in_vcf} is not compressed; gzipping..."
       bgzip "${in_vcf}"
+      vcf_to_use="${in_vcf}.gz"
     fi
+
+    # renames the seq id using the first sequence in the alignment
+    ref_name=$(head -n1 "~{ref_fasta}" | sed -E 's/^>([^[:space:]]+).*$/\1/g')
+    ref_name_no_version=$(head -n1 "~{ref_fasta}" | sed -E 's/^>([^[:space:]\.]+).*$/\1/g')
+    # copy the input or just created gzipped vcf file
+    cp "$vcf_to_use" "temp.vcf.gz"
+    # ensure uncompressed
+    bgzip -d "temp.vcf.gz"
+    # rename chr field (first col) in vcf
+    cat "temp.vcf" | sed "s/^1/$ref_name_no_version/" > "temp2.vcf"
+    
+    # output the vcf, removing the reference sequence if present as a sample name
+    bgzip "temp2.vcf"
+    tabix -p vcf "temp2.vcf.gz"
+    bcftools index "temp2.vcf.gz"
+    bcftools view -s "^$ref_name" "temp2.vcf.gz" > "temp3.vcf"
+    rm "temp2.vcf.gz"
+    vcf_to_use="temp3.vcf"
+
+    # compress vcf
+    bgzip -c "$vcf_to_use" > "$vcf_to_use.gz"
+    vcf_to_use="$vcf_to_use.gz"
+
+    # index vcf
     echo "Creating vcf index"
-    tabix -p vcf "${in_vcf}"
-        
+    bcftools index "$vcf_to_use"
+    tabix -p vcf "$vcf_to_use"
+    
     interhost.py snpEff \
-        "${in_vcf}" \
+        "$vcf_to_use" \
         $snpRefAccessions \
         "${output_basename}.annot.vcf.gz" \
         ${'--emailAddress=' + emailAddress}
-
-    intrahost.py iSNV_table \
-        "${output_basename}.annot.vcf.gz" \
-        "${output_basename}.annot.txt.gz"
-
-    tabix -p vcf "${output_basename}.annot.vcf.gz"
   }
 
   output {
     File        annot_vcf_gz      = "${output_basename}.annot.vcf.gz"
     File        annot_vcf_gz_tbi  = "${output_basename}.annot.vcf.gz.tbi"
-    File        annot_txt_gz      = "${output_basename}.annot.txt.gz"
     String      viralngs_version  = read_string("VERSION")
   }
   runtime {

pipes/WDL/workflows/mafft_and_snp_annotated.wdl

@@ -0,0 +1,76 @@
+version 1.0
+
+import "../tasks/tasks_nextstrain.wdl" as nextstrain
+import "../tasks/tasks_intrahost.wdl" as intrahost
+
+workflow mafft_and_snp_annotated {
+    meta {
+        description: "Align assemblies with mafft and find SNPs with snp-sites."
+        author: "Broad Viral Genomics"
+        email:  "viral-ngs@broadinstitute.org"
+    }
+
+    input {
+        Array[File]     assembly_fastas
+        File            ref_fasta
+        Int             min_unambig_genome
+        Boolean         run_iqtree=false
+    }
+
+    parameter_meta {
+        assembly_fastas: {
+          description: "Set of assembled genomes to align and build trees. These must represent a single chromosome/segment of a genome only. Fastas may be one-sequence-per-individual or a concatenated multi-fasta (unaligned) or a mixture of the two. They may be compressed (gz, bz2, zst, lz4), uncompressed, or a mixture.",
+          patterns: ["*.fasta", "*.fa", "*.fasta.gz", "*.fasta.zst"]
+        }
+        ref_fasta: {
+          description: "A reference assembly (not included in assembly_fastas) to align assembly_fastas against. Typically from NCBI RefSeq or similar. Uncompressed.",
+          patterns: ["*.fasta", "*.fa"]
+        }
+        min_unambig_genome: {
+          description: "Minimum number of called bases in genome to pass prefilter."
+        }
+    }
+
+    call nextstrain.gzcat {
+        input:
+            infiles     = assembly_fastas,
+            output_name = "all_samples_combined_assembly.fasta.gz"
+    }
+    call nextstrain.filter_sequences_by_length {
+        input:
+            sequences_fasta = gzcat.combined,
+            min_non_N = min_unambig_genome
+    }
+    call nextstrain.mafft_one_chr as mafft {
+        input:
+            sequences = filter_sequences_by_length.filtered_fasta,
+            ref_fasta = ref_fasta,
+            basename  = "all_samples_aligned.fasta"
+    }
+    call nextstrain.snp_sites {
+        input:
+            msa_fasta = mafft.aligned_sequences
+    }
+
+    call intrahost.annotate_vcf_snpeff as annotate_vcf {
+        input:
+            ref_fasta = ref_fasta,
+            in_vcf    = snp_sites.snps_vcf
+    }
+
+    if(run_iqtree) {
+        call nextstrain.draft_augur_tree {
+            input:
+                msa_or_vcf = mafft.aligned_sequences
+        }
+    }
+
+    output {
+        File  combined_assemblies         = filter_sequences_by_length.filtered_fasta
+        File  multiple_alignment          = mafft.aligned_sequences
+        File  unmasked_snps               = snp_sites.snps_vcf
+        File  unmasked_snps_annotated     = annotate_vcf.annot_vcf_gz
+        File  unmasked_snps_annotated_tbi = annotate_vcf.annot_vcf_gz_tbi
+        File? ml_tree                     = draft_augur_tree.aligned_tree
+    }
+}

pipes/WDL/workflows/merge_vcfs_and_annotate.wdl

@@ -33,6 +33,5 @@ workflow merge_vcfs_and_annotate {
         File merged_vcf_gz_tbi       = merge_vcfs.merged_vcf_gz_tbi
         File merged_annot_vcf_gz     = annotate_vcf.annot_vcf_gz
         File merged_annot_vcf_gz_tbi = annotate_vcf.annot_vcf_gz_tbi
-        File merged_annot_txt_gz     = annotate_vcf.annot_txt_gz
     }
 }