Merge pull request #200 from broadinstitute/dp-consensus

sarscov2 large workflow improvements

Author: dpark01

pipes/WDL/tasks/tasks_assembly.wdl

@@ -418,7 +418,7 @@ task refine_assembly_with_aligned_reads {
 
     runtime {
         docker: "${docker}"
-        memory: select_first([machine_mem_gb, 7]) + " GB"
+        memory: select_first([machine_mem_gb, 15]) + " GB"
         cpu: 8
         disks: "local-disk 375 LOCAL"
         dx_instance_type: "mem1_ssd1_v2_x8"

pipes/WDL/tasks/tasks_ncbi.wdl

@@ -252,8 +252,8 @@ task prefix_fasta_header {
   input {
     File    genome_fasta
     String  prefix
+    String  out_basename = basename(genome_fasta, ".fasta")
   }
-  String  out_basename = basename(genome_fasta, ".fasta")
   command <<<
     set -e
     python3 <<CODE
@@ -308,6 +308,12 @@ task gisaid_meta_prep {
     String out_name
     String continent = "North America"
     Boolean strict = true
+    String? username
+    String? submitting_lab_name
+    String? submitting_lab_addr
+    String? originating_lab_addr
+    String? authors
+    String? fasta_filename
   }
   command <<<
     python3 << CODE
@@ -347,12 +353,12 @@ task gisaid_meta_prep {
             'covv_seq_technology': sample_to_cmt[row['Sequence_ID']]['Sequencing Technology'],
 
             'covv_orig_lab': row['collected_by'],
-            'covv_subm_lab': 'REQUIRED',
-            'covv_authors': 'REQUIRED',
-            'covv_orig_lab_addr': 'REQUIRED',
-            'covv_subm_lab_addr': 'REQUIRED',
-            'submitter': 'REQUIRED',
-            'fn': 'REQUIRED',
+            'covv_subm_lab': "~{default='REQUIRED' submitting_lab_name}",
+            'covv_authors': "~{default='REQUIRED' authors}",
+            'covv_orig_lab_addr': "~{default='REQUIRED' originating_lab_addr}",
+            'covv_subm_lab_addr': "~{default='REQUIRED' submitting_lab_addr}",
+            'submitter': "~{default='REQUIRED' username}",
+            'fn': "~{default='REQUIRED' fasta_filename}",
           })
 
           #covv_specimen
@@ -551,10 +557,12 @@ task biosample_to_genbank {
         --biosample_in_smt \
         --iso_dates \
         --loglevel DEBUG
+    cut -f 1 "${base}.genbank.src" | tail +2 > "${base}.sample_ids.txt"
   }
   output {
     File genbank_source_modifier_table = "${base}.genbank.src"
     File biosample_map                 = "${base}.biosample.map.txt"
+    File sample_ids                    = "${base}.sample_ids.txt"
   }
   runtime {
     docker: docker

pipes/WDL/tasks/tasks_ncbi_tools.wdl

@@ -146,6 +146,7 @@ task group_sra_bams_by_biosample {
     Array[String] biosamples
     Array[File]   biosample_attributes_jsons
     Array[String] library_strategies
+    Array[String] seq_platforms
   }
   parameter_meta {
     bam_filepaths: {
@@ -165,14 +166,16 @@ task group_sra_bams_by_biosample {
     biosample_accs = '~{sep="*" biosamples}'.split('*')
     attributes = '~{sep="*" biosample_attributes_jsons}'.split('*')
     libstrats = '~{sep="*" library_strategies}'.split('*')
-    assert len(bam_uris) == len(biosample_accs) == len(attributes) == len(libstrats)
+    seqplats = '~{sep="*" seq_platforms}'.split('*')
+    assert len(bam_uris) == len(biosample_accs) == len(attributes) == len(libstrats) == len(seqplats)
 
     # lookup table files to dicts
     sample_to_bams = {}
     sample_to_attributes = {}
     sample_to_libstrat = {}
+    sample_to_seqplat = {}
     attr_keys = set()
-    for samn,bam,attr_file,libstrat in zip(biosample_accs,bam_uris, attributes, libstrats):
+    for samn,bam,attr_file,libstrat,seqplat in zip(biosample_accs,bam_uris, attributes, libstrats, seqplats):
       sample_to_bams.setdefault(samn, [])
       sample_to_bams[samn].append(bam)
       with open(attr_file, 'rt') as inf:
@@ -181,24 +184,30 @@ task group_sra_bams_by_biosample {
       attr_keys.update(k for k,v in attr.items() if v)
       sample_to_libstrat.setdefault(samn, set())
       sample_to_libstrat[samn].add(libstrat)
+      sample_to_seqplat.setdefault(samn, set())
+      sample_to_seqplat[samn].add(seqplat)
 
     # write outputs
-    with open('attributes.json', 'wt') as out_attr:
-        json.dump(sample_to_attributes, out_attr)
-    with open('attributes.tsv', 'wt') as out_attr:
+    with open('attributes.json', 'wt') as outf:
+        json.dump(sample_to_attributes, outf)
+    with open('attributes.tsv', 'wt') as outf:
         headers = tuple(sorted(attr_keys))
-        out_attr.write('\t'.join(headers)+'\n')
+        outf.write('\t'.join(headers)+'\n')
         for sample in sorted(sample_to_bams.keys()):
-            out_attr.write('\t'.join(sample_to_attributes[sample].get(h,'') for h in headers)+'\n')
+            outf.write('\t'.join(sample_to_attributes[sample].get(h,'') for h in headers)+'\n')
     with open('grouped_bams', 'wt') as out_groups:
       with open('samns', 'wt') as out_samples:
         for sample in sorted(sample_to_bams.keys()):
           out_samples.write(sample+'\n')
           out_groups.write('\t'.join(sample_to_bams[sample])+'\n')
-    with open('library_strategies.json', 'wt') as out_attr:
+    with open('library_strategies.json', 'wt') as outf:
         for k,v in sample_to_libstrat.items():
             sample_to_libstrat[k] = ';'.join(sorted(v))
-        json.dump(sample_to_libstrat, out_attr)
+        json.dump(sample_to_libstrat, outf)
+    with open('sequencing_platforms.json', 'wt') as outf:
+        for k,v in sample_to_seqplat.items():
+            sample_to_seqplat[k] = ';'.join(sorted(v))
+        json.dump(sample_to_seqplat, outf)
     CODE
   >>>
   output {
@@ -207,6 +216,7 @@ task group_sra_bams_by_biosample {
     Map[String,Map[String,String]] samn_to_attributes = read_json('attributes.json')
     File                biosample_attributes_tsv = 'attributes.tsv'
     Map[String,String]  samn_to_library_strategy = read_json('library_strategies.json')
+    Map[String,String]  samn_to_sequencing_platform = read_json('sequencing_platforms.json')
   }
   runtime {
     docker: "python:slim"

pipes/WDL/tasks/tasks_read_utils.wdl

@@ -12,7 +12,7 @@ task max {
     CODE
   >>>
   output {
-    Int max = read_int(stdout())
+    Int out = read_int(stdout())
   }
   runtime {
     docker: "python:slim"

pipes/WDL/workflows/demux_deplete.wdl

@@ -152,6 +152,7 @@ workflow demux_deplete {
         Array[Int]  read_counts_depleted = deplete.depletion_read_count_post
 
         File?       sra_metadata          = sra_meta_prep.sra_metadata
+        File?       cleaned_bam_uris      = sra_meta_prep.cleaned_bam_uris
 
         Array[File] demux_metrics         = illumina_demux.metrics
         Array[File] demux_commonBarcodes  = illumina_demux.commonBarcodes

pipes/WDL/workflows/sarscov2_illumina_full.wdl

@@ -3,6 +3,7 @@ version 1.0
 import "../tasks/tasks_read_utils.wdl" as read_utils
 import "../tasks/tasks_ncbi.wdl" as ncbi
 import "../tasks/tasks_nextstrain.wdl" as nextstrain
+import "../tasks/tasks_reports.wdl" as reports
 
 import "demux_deplete.wdl"
 import "assemble_refbased.wdl"
@@ -38,22 +39,33 @@ workflow sarscov2_illumina_full {
     }
 
     input {
+        File          flowcell_tgz
         File          reference_fasta
         String        amplicon_bed_prefix
 
-        File          biosample_attributes
+        Array[File]   biosample_attributes
         String        instrument_model
         String        sra_title
 
-        Int           min_genome_bases = 20000
+        Int           min_genome_bases = 15000
     }
     Int     taxid = 2697049
     String  gisaid_prefix = 'hCoV-19/'
+    String  flowcell_id = basename(basename(basename(basename(flowcell_tgz, ".gz"), ".zst"), ".tar"), ".tgz")
+
+    # merge biosample attributes tables
+    call reports.tsv_join as biosample_merge {
+        input:
+            input_tsvs = biosample_attributes,
+            id_col = 'accession',
+            out_basename = "biosample_attributes-merged"
+    }
 
     ### demux, deplete, SRA submission prep, fastqc/multiqc
     call demux_deplete.demux_deplete {
         input:
-            biosample_map = biosample_attributes,
+            flowcell_tgz = flowcell_tgz,
+            biosample_map = biosample_merge.out_tsv,
             instrument_model = instrument_model,
             sra_title = sra_title
     }
@@ -99,7 +111,7 @@ workflow sarscov2_illumina_full {
 
             File passing_assemblies = rename_fasta_header.renamed_fasta
             String passing_assembly_ids = orig_name
-            Array[String] assembly_cmt = [orig_name, "Broad viral-ngs v. " + demux_deplete.demux_viral_core_version, assemble_refbased.assembly_mean_coverage]
+            Array[String] assembly_cmt = [orig_name, "Broad viral-ngs v. " + demux_deplete.demux_viral_core_version, assemble_refbased.assembly_mean_coverage, instrument_model]
 
             # lineage assignment
             call sarscov2_lineages.sarscov2_lineages {
@@ -124,81 +136,100 @@ workflow sarscov2_illumina_full {
             String failed_assembly_id = orig_name
         }
 
-        Map[String,String?] assembly_stats = {
-            'sample_orig': orig_name,
-            'sample': name_reads.left,
-            'amplicon_set': demux_deplete.meta_by_sample[name_reads.left]["amplicon_set"],
-            'assembly_mean_coverage': assemble_refbased.assembly_mean_coverage,
-            'nextclade_clade':   sarscov2_lineages.nextclade_clade,
-            'nextclade_aa_subs': sarscov2_lineages.nextclade_aa_subs,
-            'nextclade_aa_dels': sarscov2_lineages.nextclade_aa_dels,
-            'pango_lineage':     sarscov2_lineages.pango_lineage
-        }
-        Map[String,File?] assembly_files = {
-            'assembly_fasta':           assemble_refbased.assembly_fasta,
-            'coverage_plot':            assemble_refbased.align_to_ref_merged_coverage_plot,
-            'aligned_bam':              assemble_refbased.align_to_ref_merged_aligned_trimmed_only_bam,
-            'replicate_discordant_vcf': assemble_refbased.replicate_discordant_vcf,
-            'nextclade_tsv': sarscov2_lineages.nextclade_tsv,
-            'pangolin_csv':  sarscov2_lineages.pangolin_csv,
-            'vadr_tgz': vadr.outputs_tgz
-        }
-        Map[String,Int?] assembly_metrics = {
-            'assembly_length_unambiguous': assemble_refbased.assembly_length_unambiguous,
-            'dist_to_ref_snps':            assemble_refbased.dist_to_ref_snps,
-            'dist_to_ref_indels':          assemble_refbased.dist_to_ref_indels,
-            'replicate_concordant_sites':  assemble_refbased.replicate_concordant_sites,
-            'replicate_discordant_snps':   assemble_refbased.replicate_discordant_snps,
-            'replicate_discordant_indels': assemble_refbased.replicate_discordant_indels,
-            'num_read_groups':             assemble_refbased.num_read_groups,
-            'num_libraries':               assemble_refbased.num_libraries,
-            'vadr_num_alerts': vadr.num_alerts
-        }
-
+        Array[String] assembly_tsv_row = [
+            orig_name,
+            name_reads.left,
+            demux_deplete.meta_by_sample[name_reads.left]["amplicon_set"],
+            assemble_refbased.assembly_mean_coverage,
+            assemble_refbased.assembly_length_unambiguous,
+            select_first([sarscov2_lineages.nextclade_clade, ""]),
+            select_first([sarscov2_lineages.nextclade_aa_subs, ""]),
+            select_first([sarscov2_lineages.nextclade_aa_dels, ""]),
+            select_first([sarscov2_lineages.pango_lineage, ""]),
+            assemble_refbased.dist_to_ref_snps,
+            assemble_refbased.dist_to_ref_indels,
+            select_first([vadr.num_alerts, ""]),
+            assemble_refbased.assembly_fasta,
+            assemble_refbased.align_to_ref_merged_coverage_plot,
+            assemble_refbased.align_to_ref_merged_aligned_trimmed_only_bam,
+            assemble_refbased.replicate_discordant_vcf,
+            select_first([sarscov2_lineages.nextclade_tsv, ""]),
+            select_first([sarscov2_lineages.pangolin_csv, ""]),
+            select_first([vadr.outputs_tgz, ""]),
+            assemble_refbased.replicate_concordant_sites,
+            assemble_refbased.replicate_discordant_snps,
+            assemble_refbased.replicate_discordant_indels,
+            assemble_refbased.num_read_groups,
+            assemble_refbased.num_libraries,
+        ]
+    }
+    Array[String] assembly_tsv_header = [
+        'sample', 'sample_sanitized', 'amplicon_set', 'assembly_mean_coverage', 'assembly_length_unambiguous',
+        'nextclade_clade', 'nextclade_aa_subs', 'nextclade_aa_dels', 'pango_lineage',
+        'dist_to_ref_snps', 'dist_to_ref_indels', 'vadr_num_alerts',
+        'assembly_fasta', 'coverage_plot', 'aligned_bam', 'replicate_discordant_vcf',
+        'nextclade_tsv', 'pangolin_csv', 'vadr_tgz',
+        'replicate_concordant_sites', 'replicate_discordant_snps', 'replicate_discordant_indels', 'num_read_groups', 'num_libraries',
+        ]
+
+    call nextstrain.concatenate as assembly_meta_tsv {
+      input:
+        infiles = [write_tsv([assembly_tsv_header]), write_tsv(assembly_tsv_row)],
+        output_name = "assembly_metadata-~{flowcell_id}.tsv"
     }
 
-    # TO DO: filter out genomes from submission that are less than ntc_bases.max
-    call read_utils.max as ntc {
+
+    # TO DO: filter out genomes from submission that are less than ntc_max.out
+    call read_utils.max as ntc_max {
       input:
         list = select_all(ntc_bases)
     }
 
     ### prep genbank submission
-    call nextstrain.concatenate as submit_genomes {
-      input:
-        infiles = select_all(submittable_genomes),
-        output_name = "assemblies.fasta"
-    }
     call ncbi.biosample_to_genbank {
       input:
-        biosample_attributes = biosample_attributes,
+        biosample_attributes = biosample_merge.out_tsv,
         num_segments = 1,
         taxid = taxid,
         filter_to_ids = write_lines(select_all(submittable_id))
     }
     call ncbi.structured_comments {
       input:
-        assembly_stats_tsv = write_tsv(flatten([[['SeqID','Assembly Method','Coverage']],select_all(assembly_cmt)])),
-        filter_to_ids = write_lines(select_all(submittable_id))
+        assembly_stats_tsv = write_tsv(flatten([[['SeqID','Assembly Method','Coverage','Sequencing Technology']],select_all(assembly_cmt)])),
+        filter_to_ids = biosample_to_genbank.sample_ids
+    }
+    call nextstrain.concatenate as passing_genomes {
+      input:
+        infiles = select_all(submittable_genomes),
+        output_name = "assemblies.fasta"
+    }
+    call nextstrain.filter_sequences_to_list as submit_genomes {
+      input:
+        sequences = passing_genomes.combined,
+        keep_list = [biosample_to_genbank.sample_ids]
     }
     call ncbi.package_genbank_ftp_submission {
       input:
-        sequences_fasta = submit_genomes.combined,
+        sequences_fasta = submit_genomes.filtered_fasta,
         source_modifier_table = biosample_to_genbank.genbank_source_modifier_table,
-        structured_comment_table = structured_comments.structured_comment_table
+        structured_comment_table = structured_comments.structured_comment_table,
+        submission_name = flowcell_id,
+        submission_uid = flowcell_id
     }
 
     ### prep gisaid submission
     call ncbi.prefix_fasta_header as prefix_gisaid {
       input:
-        genome_fasta = submit_genomes.combined,
-        prefix = gisaid_prefix
+        genome_fasta = submit_genomes.filtered_fasta,
+        prefix = gisaid_prefix,
+        out_basename = "gisaid-sequences-~{flowcell_id}"
     }
     call ncbi.gisaid_meta_prep {
       input:
         source_modifier_table = biosample_to_genbank.genbank_source_modifier_table,
         structured_comments = structured_comments.structured_comment_table,
-        out_name = "gisaid_meta.tsv"
+        fasta_filename = "gisaid-sequences-~{flowcell_id}.fasta",
+        out_name = "gisaid-meta-~{flowcell_id}.tsv"
     }
 
     output {
@@ -212,12 +243,13 @@ workflow sarscov2_illumina_full {
         Array[Int]  read_counts_depleted  = demux_deplete.read_counts_depleted
 
         File        sra_metadata          = select_first([demux_deplete.sra_metadata])
+        File        cleaned_bam_uris      = select_first([demux_deplete.cleaned_bam_uris])
 
         Array[File] assemblies_fasta = assemble_refbased.assembly_fasta
         Array[File] passing_assemblies_fasta = select_all(passing_assemblies)
         Array[File] submittable_assemblies_fasta = select_all(submittable_genomes)
 
-        Int         max_ntc_bases = ntc.max
+        Int         max_ntc_bases = ntc_max.out
 
         Array[File] demux_metrics            = demux_deplete.demux_metrics
         Array[File] demux_commonBarcodes     = demux_deplete.demux_commonBarcodes
@@ -227,9 +259,7 @@ workflow sarscov2_illumina_full {
         File        multiqc_report_cleaned = demux_deplete.multiqc_report_cleaned
         File        spikein_counts         = demux_deplete.spikein_counts
 
-        Array[Map[String,String?]] per_assembly_stats = assembly_stats
-        Array[Map[String,File?]]   per_assembly_files = assembly_files
-        Array[Map[String,Int?]]    per_assembly_metrics = assembly_metrics
+        File assembly_stats_tsv = assembly_meta_tsv.combined
 
         File submission_zip = package_genbank_ftp_submission.submission_zip
         File submission_xml = package_genbank_ftp_submission.submission_xml