Merge pull request #800 from broadinstitute/dp-genbank

genbank and other new WDL workflows

ncbi.py

@@ -381,6 +381,7 @@ def fasta2fsa(infname, outdir, biosample=None):
 def make_structured_comment_file(cmt_fname, name=None, seq_tech=None, coverage=None):
     with open(cmt_fname, 'wt') as outf:
         outf.write("StructuredCommentPrefix\t##Genome-Assembly-Data-START##\n")
+        # note: the <tool name> v. <version name> format is required by NCBI, don't remove the " v. "
         outf.write("Assembly Method\tgithub.com/broadinstitute/viral-ngs v. {}\n".format(util.version.get_version()))
         if name:
             outf.write("Assembly Name\t{}\n".format(name))
@@ -448,7 +449,7 @@ def prep_genbank_files(templateFile, fasta_files, annotDir,
                             outf.write(inf.readline())
                             for line in inf:
                                 row = line.rstrip('\n').split('\t')
-                                if row[0] == sample_base:
+                                if row[0] == sample_base or row[0] == sample:
                                     row[0] = sample
                                     outf.write('\t'.join(row) + '\n')
 

pipes/WDL/workflows/coverage_table.wdl

@@ -0,0 +1,5 @@
+import "reports.wdl" as reports
+
+workflow coverage_table {
+  call reports.coverage_report
+}

pipes/WDL/workflows/genbank.wdl

@@ -0,0 +1,29 @@
+import "interhost.wdl" as interhost
+import "ncbi.wdl" as ncbi
+
+workflow genbank {
+
+  File          reference_fasta
+  Array[File]+  assemblies_fasta     # one per genome
+  Array[File]+  ref_annotations_tbl  # one per chromosome
+
+  call interhost.multi_align_mafft_ref as mafft {
+    input:
+      reference_fasta = reference_fasta,
+      assemblies_fasta = assemblies_fasta
+  }
+
+  call ncbi.annot_transfer as annot {
+    input:
+      multi_aln_fasta = mafft.alignments_by_chr,
+      reference_fasta = reference_fasta,
+      reference_feature_table = ref_annotations_tbl
+  }
+
+  call ncbi.prepare_genbank as prep_genbank {
+    input:
+      assemblies_fasta = assemblies_fasta,
+      annotations_tbl = annot.transferred_feature_tables
+  }
+
+}

pipes/WDL/workflows/mafft.wdl

@@ -0,0 +1,5 @@
+import "interhost.wdl" as interhost
+
+workflow mafft {
+  call interhost.multi_align_mafft
+}

pipes/WDL/workflows/merge_bams.wdl

@@ -0,0 +1,5 @@
+import "demux.wdl" as demux
+
+workflow merge_bams {
+  call demux.merge_and_reheader_bams
+}

pipes/WDL/workflows/tasks/demux.wdl

@@ -107,3 +107,38 @@ task illumina_demux {
     preemptible: 0  # this is the very first operation before scatter, so let's get it done quickly & reliably
   }
 }
+
+task merge_and_reheader_bams {
+  Array[File]+  in_bams
+  File?         reheader_table # tsv with 3 cols: field, old value, new value
+  String        out_basename
+
+  command {
+    set -ex -o pipefail
+
+    if [ ${length(in_bams)} -gt 1 ]; then
+      read_utils.py merge_bams ${sep=' ' in_bams} merged.bam --loglevel DEBUG
+    else
+      echo "Skipping merge, only one input file"
+      ln -s ${select_first(in_bams)} merged.bam
+    fi    
+
+    if [ -f ${reheader_table} ]; then
+      read_utils.py merged.bam ${reheader_table} ${out_basename}.bam --loglevel DEBUG
+    else
+      echo "Skipping reheader, no mapping table specified"
+      ln -s merged.bam ${out_basename}.bam
+    fi
+  }
+
+  output {
+    File  out_bam = "${out_basename}.bam"
+  }
+
+  runtime {
+    docker: "quay.io/broadinstitute/viral-ngs"
+    memory: "2000 MB"
+    cpu: 2
+    dx_instance_type: "mem1_ssd2_x4"
+  }
+}

pipes/WDL/workflows/tasks/interhost.wdl

@@ -1,92 +1,18 @@
-task ref_guided_consensus {
-  String sample
-
-  File referenceGenome #fasta
-  File inBam
-  # TODO: input settings must compute chr_names
-  Array[String] chrNames # sampleName-1, sampleName-2, ...
-
-  Int? threads
-  Int? minCoverage
-
-  String? novoalignOptions
-
-  command {
-    assembly.py refine_assembly \
-      "${referenceGenome}" \
-      "${inBam}" \
-      "${sample}.fasta" \
-      --outBam "${sample}.realigned.bam" \
-      --outVcf "${sample}.vcf.gz" \
-      "${'--min_coverage' + minCoverage || '3'}" \
-      "${'--novo_params' + novoalignOptions || '-r Random -l 40 -g 40 -x 20 -t 100'}" \
-      --keep_all_reads \
-      --chr_names "${sep=' ' chrNames+}" \
-      "${'--threads' + threads}"
-  }
-
-  output {
-    File assembly = "${sample}.fasta"
-    File realignedBam = "${sample}.realigned.bam"
-    File variantCalls = "${sample}.vcf.gz"
-  }
-  runtime {
-    memory: "4 GB"
-    docker: "quay.io/broadinstitute/viral-ngs"
-  }
-}
-
-task ref_guided_consensus_aligned_with_dups {
-  String sample
-
-  File realignedBam
-
-  command {
-    samtools view -b -F 4 -o "${sample}.realigned.only_aligned.bam" "${realignedBam}"
-  }
-
-  output {
-    File onlyAlignedReadsBam = "${sample}.realigned.only_aligned.bam"
-  }
-  runtime {
-    memory: "8 GB"
-    docker: "quay.io/broadinstitute/viral-ngs"
-  }
-}
-
-#task ref_guided_diversity {
-#  String sample
-#
-#  File assembly # .fasta
-#  File variantCalls # .vcf.gz
-#  File referenceGenome # .fasta, but also .fasta.fai and .dict
-#
-#  command {
-#    GenomeAnalysisTK.jar
-#      -T CombineVariants -R "${referenceGenome}" {inFilesString} -o {outFile}
-#      --genotypemergeoption UNIQUIFY
-#  }
-#
-#  output {
-#
-#  }
-#  runtime {
-#    docker: "quay.io/broadinstitute/viral-ngs"
-#  }
-#}
 
 task multi_align_mafft_ref {
   File           reference_fasta
   Array[File]+   assemblies_fasta # fasta files, one per sample, multiple chrs per file okay
-  String?        out_prefix = basename(reference_fasta, '.fasta')
+  String         out_prefix = basename(reference_fasta, '.fasta')
   Int?           mafft_maxIters
-  Int?           mafft_ep
+  Float?         mafft_ep
+  Float?         mafft_gapOpeningPenalty
 
   command {
     interhost.py multichr_mafft \
       ${reference_fasta} ${sep=' ' assemblies_fasta} \
       . \
       ${'--ep' + mafft_ep} \
+      ${'--gapOpeningPenalty' + mafft_gapOpeningPenalty} \
       ${'--maxiters' + mafft_maxIters} \
       --outFilePrefix ${out_prefix} \
       --preservecase \
@@ -96,29 +22,31 @@ task multi_align_mafft_ref {
   }
 
   output {
-    File        sampleNamesFile   = "${out_prefix}-sample_names.txt"
-    Array[File] alignments_by_chr = glob("${out_prefix}*.fasta")
+    File         sampleNamesFile   = "${out_prefix}-sample_names.txt"
+    Array[File]+ alignments_by_chr = glob("${out_prefix}*.fasta")
   }
 
   runtime {
     docker: "quay.io/broadinstitute/viral-ngs"
     memory: "7 GB"
-    cpu: 4
-    dx_instance_type: "mem1_ssd1_x4"
+    cpu: 8
+    dx_instance_type: "mem1_ssd1_x8"
   }
 }
 
 task multi_align_mafft {
   Array[File]+   assemblies_fasta # fasta files, one per sample, multiple chrs per file okay
   String         out_prefix
   Int?           mafft_maxIters
-  Int?           mafft_ep
+  Float?         mafft_ep
+  Float?         mafft_gapOpeningPenalty
 
   command {
     interhost.py multichr_mafft \
       ${sep=' ' assemblies_fasta} \
       . \
       ${'--ep' + mafft_ep} \
+      ${'--gapOpeningPenalty' + mafft_gapOpeningPenalty} \
       ${'--maxiters' + mafft_maxIters} \
       --outFilePrefix ${out_prefix} \
       --preservecase \
@@ -135,8 +63,8 @@ task multi_align_mafft {
   runtime {
     docker: "quay.io/broadinstitute/viral-ngs"
     memory: "7 GB"
-    cpu: 4
-    dx_instance_type: "mem1_ssd1_x4"
+    cpu: 8
+    dx_instance_type: "mem1_ssd1_x8"
   }
 }
 

pipes/WDL/workflows/tasks/ncbi.wdl

@@ -58,22 +58,31 @@ task download_annotations {
 }
 
 task annot_transfer {
-  File chr_mutli_aln_fasta # fasta; multiple alignments of sample sequences for a single chr
-  File reference_fasta # fasta (may contain multiple chrs, only one with the same name as reference_feature_table will be used)
-  File reference_feature_table # feature table corresponding to the chr in the alignment
+  Array[File]+ multi_aln_fasta         # fasta; multiple alignments of sample sequences for each chromosome
+  File         reference_fasta         # fasta; all chromosomes in one file
+  Array[File]+ reference_feature_table # tbl; feature table corresponding to each chromosome in the alignment
+
+  Array[Int]   chr_nums=range(length(multi_aln_fasta))
 
   command {
-    ncbi.py tbl_transfer_prealigned \
-        ${chr_mutli_aln_fasta} \
-        ${reference_fasta} \
-        ${reference_feature_table} \
-        . \
-        --oob_clip \
-        --loglevel DEBUG
+    set -ex -o pipefail
+    echo ${sep=' ' multi_aln_fasta} > alignments.txt
+    echo ${sep=' ' reference_feature_table} > tbls.txt
+    for i in ${sep=' ' chr_nums}; do
+      _alignment_fasta=`cat alignments.txt | cut -f $(($i+1)) -d ' '`
+      _feature_tbl=`cat tbls.txt | cut -f $(($i+1)) -d ' '`
+      ncbi.py tbl_transfer_prealigned \
+          $_alignment_fasta \
+          ${reference_fasta} \
+          $_feature_tbl \
+          . \
+          --oob_clip \
+          --loglevel DEBUG
+    done
   }
 
   output {
-    Array[File] transferred_feature_tables = glob("*.tbl")
+    Array[File]+ transferred_feature_tables = glob("*.tbl")
   }
   runtime {
     docker: "quay.io/broadinstitute/viral-ngs"
@@ -87,12 +96,13 @@ task prepare_genbank {
   Array[File]+ assemblies_fasta
   Array[File]+ annotations_tbl
   File         authors_sbt
-  File?        coverage_table # summary.assembly.txt (from Snakemake)
-  File?        genbankSourceTable
-  File?        biosampleMap
-  String?      sequencingTech
-  String?      comment
-  String       out_prefix = "ncbi_package"
+  File         biosampleMap
+  File         genbankSourceTable
+  File?        coverage_table # summary.assembly.txt (from Snakemake) -- change this to accept a list of mapped bam files and we can create this table ourselves
+  String       sequencingTech
+  String       comment # TO DO: make this optional
+  String       organism
+  String       molType = "cRNA"
 
   command {
     set -ex -o pipefail
@@ -101,19 +111,25 @@ task prepare_genbank {
         ${authors_sbt} \
         ${sep=' ' assemblies_fasta} \
         . \
-        ${'--master_source_table=' + genbankSourceTable} \
-        ${'--sequencing_tech=' + sequencingTech} \
-        ${'--biosample_map=' + biosampleMap} \
-        ${'--coverage_table=' + coverage_table} \
-        ${'--comment=' + comment} \
+        --mol_type ${molType} \
+        --organism "${organism}" \
+        --biosample_map ${biosampleMap} \
+        --master_source_table ${genbankSourceTable} \
+        ${'--coverage_table ' + coverage_table} \
+        --comment "${comment}" \
+        --sequencing_tech "${sequencingTech}" \
         --loglevel DEBUG
-    tar -czpvf ${out_prefix}.tar.gz *.val *.cmt *.fsa *.gbf *.sqn *.src *.tbl
+    mv errorsummary.val errorsummary.val.txt # to keep it separate from the glob
   }
 
   output {
     Array[File] sequin_files = glob("*.sqn")
-    File        ncbi_package = "${out_prefix}.tar.gz"
-    File        errorSummary = "errorsummary.val"
+    Array[File] structured_comment_files = glob("*.cmt")
+    Array[File] genbank_preview_files = glob("*.gbf")
+    Array[File] source_table_files = glob("*.src")
+    Array[File] fasta_per_chr_files = glob("*.fsa")
+    Array[File] validation_files = glob("*.val")
+    File        errorSummary = "errorsummary.val.txt"
   }
 
   runtime {
@@ -122,4 +138,4 @@ task prepare_genbank {
     cpu: 2
     dx_instance_type: "mem1_ssd1_x2"
   }
-}
+}