Shared variants tasks and QC improvements for kSNP3 and Snippy

[michellescribner]

Closes https://github.com/theiagen/public_health_bioinformatics/issues/258

:hammer_and_wrench: Changes Being Made

This PR makes a number of changes to kSNP3 and Snippy-related workflows with the broad goal of creating summary files that show the SNPs shared among a set of samples and adding relevant QC assessments.

Adds a cat_files task which will concatenate variant files

• in cat_files file /tasks/utilities/file_handling/task_cat_files.wdl, adds a task cat_variants, which is similar to the cat_files task except adds the sample names as an additional column in the concatenated file

samplename	CHROM	POS	TYPE	REF	ALT	EVIDENCE	FTYPE	STRAND	NT_POS	AA_POS	EFFECT	LOCUS_TAG	GENE	PRODUCT
sample1	PEKT02000007	5224	snp	C	G	G:21 C:0
sample2	PEKT02000007	34112	snp	C	G	G:32 C:0	CDS	+	153/1620	51/539	missense_variant c.153C>G p.His51Gln	B9J08_002604	hypothetical protein
sample3	PEKT02000007	34487	snp	T	A	A:41 T:0	CDS	+	528/1620	176/539	missense_variant c.528T>A p.Asn176Lys	B9J08_002604	hypothetical protein

Adds a new task tasks/phylogenetic_inference/utilities/task_shared_variants.wdl

• reshapes concatenated variants file into shared variants table such that samples are columns and presence and absence of each mutation is indicated by 1s and 0s shared_variants_table.

CHROM	POS	TYPE	REF	ALT	FTYPE	STRAND	NT_POS	AA_POS	EFFECT	LOCUS_TAG	GENE	PRODUCT	sample1	sample2	sample3
PEKT02000007	2693938	snp	T	C	CDS	-	1008/3000	336/999	synonymous_variant c.1008A>G p.Lys336Lys	B9J08_003879	NA	chitin synthase 1	1	1	0
PEKT02000007	2529234	snp	G	C	CDS	+	282/336	94/111	missense_variant c.282G>C p.Lys94Asn	B9J08_003804	NA	cytochrome c	1	1	1
PEKT02000002	1043926	snp	A	G	CDS	-	542/1464	181/487	missense_variant c.542T>C p.Ile181Thr	B9J08_000976	NA	dihydrolipoyl dehydrogenase	1	1	0

Adds Find_Shared_Variants_PHB standalone workflow

• Adds a workflow that will concatenate the csv outputs of the snippy variants task into a single csv file using the cat_variants task above. It then reshapes this output using the shared_variants task above.

Adds cat_variants task and shared_variants task to workflows/phylogenetics/wf_snippy_tree.wdl as an optional modules

• adds input boolean call_shared_variants = true
• modifies the snippy core task tasks/phylogenetic_inference/task_snippy_core.wdl to glob the snippy results csvs after untar-ing snippy results tarballs

Adds snippy variants task QC improvements

• computes and outputs number or reads aligned to reference genome in snippy variants task using samtools view -c
• creates and outputs a coverage stats file for each sample using samtools coverage
• parses and outputs the number of variants detected in a sample from snippy summary file
• creates and outputs a file listing read depth for every position in reference genome using samtools depth
• computes and outputs proportion of genome with depth >= min_cov

Adds kSNP3 task QC improvements

• exposes number of SNPs in the sample set ksnp3_number_snps
• exposes number of core SNPs in the sample set ksnp3_number_core_snps

Adds kSNP3 Shared SNP Task

• adds new task tasks/phylogenetic_inference/utilities/task_ksnp3_shared_snps.wdl which takes in ksnp3_vcf_ref_genome file as input, filters for only core SNPs, sorts SNPs by number of occurrences in the sample set, and outputs the table as a csv file ksnp3_core_snp_table
• adds task to kSNP3 workflow

Impacted Workflows/Tasks

The following tasks and workflows are modified based on the changes described above. All TheiaProk workflows are also impacted due to the changes to Merlin Magic.

• Modifies file tasks/utilities/file_handling/task_cat_files.wdl to add new task cat_variants

• Adds new task tasks/phylogenetic_inference/utilities/task_shared_variants.wdl

• Adds new workflow workflows/utilities/file_handling/wf_find_shared_variants.wdl

• Adds Find_Shared_Variants_PHB workflow to dockstore.yml

• Adds shared variants task to Snippy Tree wf workflows/phylogenetics/wf_snippy_tree.wdl

• Globs results csvs from snippy core task for downstream use by concatenate variants tasks/phylogenetic_inference/task_snippy_core.wdl

• Adds new qc columns to tasks/gene_typing/task_snippy_variants.wdl

• Adds snippy variants qc columns to workflows/standalone_modules/wf_snippy_variants.wdl

• Adds snippy variants qc columns to workflows/utilities/wf_merlin_magic.wdl

• Adds snippy variants qc columns to workflows/theiaeuk/wf_theiaeuk_illumina_pe.wdl

• Adds new qc columns to tasks/phylogenetic_inference/task_ksnp3.wdl

• Adds new task for creating shared SNP table tasks/phylogenetic_inference/utilities/task_ksnp3_shared_snps.wdl

• Adds shared SNP task to kSNP3 workflow and new qc columns workflows/phylogenetics/wf_ksnp3.wdl

:brain: Context and Rationale

Broadly, all of these changes were made to improve the ability of the user to assess the quality of phylogenetic analyses using the kSNP3 and Snippy workflows.

• The find shared variants workflow was needed because the existing workflows/utilities/file_handling/wf_concatenate_column.wdl workflow does not include sample information in the final csv. This workflow will therefore enable users to concatenate snippy results across multiple samples effectively.
• The new QC columns for Snippy variants is intended to expose when a suspiciously low proportion of reads align to a reference genome, a very high number of mutations are detected, or the percent of the reference genome with adequate coverage is low. These pieces of data may reveal instances of sample contamination, poor choice of reference genome, or low genome coverage.
• The kSNP3 shared SNP task was requested by a public health partner and will show which mutations in a sample set are shared among multiple samples.

:clipboard: Workflow/Task Steps

Inputs

Outputs

New outputs for workflows that invoke snippy tree: (snippy tree, snippy streamline)

- snippy_concatenated_variants
- snippy_shared_variants

New outputs for workflows that invoke the snippy variants task: (snippy_variants standalone, theiaeuk, merlin magic)

- snippy_variants_num_reads_aligned
- snippy_variants_num_variants
- snippy_variants_coverage_tsv
- snippy_variants_percent_ref_coverage

New outputs for kSNP3:

- ksnp3_number_snps
- ksnp3_number_core_snps
- ksnp3_core_snp_table

Impacted Outputs

None

:test_tube: Testing

Locally

Terra

• Find_Shared_Variants_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/70bd8836-d821-4fb9-963a-1301ea93fd9e
• Snippy_Variants_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/941d4fea-6abf-4672-a139-4c8a52b51c49
• Snippy_Tree_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/c5b11c65-7bb8-416d-95d6-3413b4888818
• Snippy_Streamline_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/97786cf5-7fda-40bc-9f17-61f7e4077c8b
• TheiaEuk_Illumina_PE_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/4d4d77b7-ee42-4a80-9cb0-7ac8d7d869cf
• kSNP3_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/4bfb1885-4663-495d-b20e-a4178d3d2c96
• TheiaProk_Illumina_PE_PHB (a few samples failed at first launch but passed after relaunching failures): https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/9a97fb52-16c7-4732-8d35-8c59a04eede4
• TheiaProk_Illlumina_SE_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/696879bc-97a9-4612-b35f-52812b18bec0
• TheiaProk_FASTA_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/83267716-d1b0-4777-a33d-471c5af271a9
• TheiaProk_ONT_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/3bc3bf02-4b5d-4b9f-ad46-59e914a28726

Scenarios for Reviewer to Test

• [ ] Snippy_Variants, checking expected outputs as below are produced and make sense
  • snippy_variants_num_reads_aligned
  • snippy_variants_num_variants
  • snippy_variants_coverage_tsv
  • snippy_variants_percent_ref_coverage
• [ ] Snippy_Streamline, checking expected outputs as below are produced and make sense
  • snippy_concatenated_variants
  • snippy_shared_variants_table
• [ ] Snippy_Tree, checking expected outputs as below are produced and make sense
  • snippy_concatenated_variants
  • snippy_shared_variants_table
• [ ] Find_Shared_Variants, checking expected outputs as below are produced and make sense
  • concatenated_variants
  • shared_variants_table
• [ ] kSNP3, checking expected outputs as below are produced and make sense
  • ksnp3_number_snps
  • ksnp3_number_core_snps
  • ksnp3_core_snp_table
• [ ] TheiaEuk, checking expected outputs as below are produced and make sense
  • snippy_variants_num_reads_aligned
  • snippy_variants_num_variants
  • snippy_variants_coverage_tsv
  • snippy_variants_percent_ref_coverage
• [ ] workflows/phylogenetics/wf_find_shared_variants.wdl to ensure sample names column is produced in concatenated_variants file and shared variants table is produced

:microscope: Final Developer Checklist

• [ ] The workflow/task has been tested locally and results, including file contents, are as anticipated
• [ ] The workflow/task has been tested on Terra and results, including file contents, are as anticipated
• [ ] The CI/CD has been adjusted and tests are passing (to be completed by Theiagen developer)
• [x] Code changes follow the style guide

🎯 Reviewer Checklist

• [ ] All impacted workflows/tasks have been tested on Terra with a different dataset than used for development
• [ ] All reviewer-suggested scenarios have been tested and any additional
• [ ] All changed results have been confirmed to be accurate
• [ ] All workflows/tasks impacted by change/s have been tested using a standard validation dataset to ensure no unintended change of functionality
• [ ] All code adheres to the style guide
• [ ] MD5 sums have been updated
• [ ] The PR author has addressed all comments

[emmadoughty]

This is a beast! Initiating function testing for the following scenarios and datasets:

• [ ] Snippy_Variants, salmonella dataset, checking expected outputs as below are produced and make sense
  • snippy_variants_num_reads_aligned
  • snippy_variants_num_variants
  • snippy_variants_coverage_tsv
  • snippy_variants_percent_ref_coverage
• [x] Snippy_Streamline, HAV dataset, checking expected outputs as below are produced and make sense (HAV specifically requested for review by Emily on behalf of a PHL)
  • snippy_concatenated_snps
  • snippy_shared_snps
• [ ] Snippy_Streamline, Salmonella dataset, checking expected outputs as below are produced and make sense
  • snippy_concatenated_snps
  • snippy_shared_snps
• [ ] Snippy_Tree, salmonella dataset, checking expected outputs as below are produced and make sense
  • snippy_concatenated_snps
  • snippy_shared_snps
• [ ] kSNP3, salmonella dataset, checking expected outputs as below are produced and make sense
  • ksnp3_number_snps
  • ksnp3_number_core_snps
  • ksnp3_core_snp_table
• [ ] TheiaEuk, fungal dataset TBD, checking expected outputs as below are produced and make sense
  • snippy_variants_num_reads_aligned
  • snippy_variants_num_variants
  • snippy_variants_coverage_tsv
  • snippy_variants_percent_ref_coverage
• [ ] workflows/utilities/file_handling/wf_concatenate_column.wdl to ensure sample names column is produced when optional input used

Do not intend to test:

• Merlin magic stand-alone
• workflows/utilities/file_handling/wf_find_shared_variants.wdl (As it's a sub-wf that's not intended to be run independently?)

@michellescribner, does this cover all the scenarios to test? Is the Find_Shared_Variants_PHB workflow intended to be run independently by users?

[emmadoughty]

@michellescribner This PR has a lot to test, so I am providing some initial feedback for discussion before testing all scenarios:

• When setting up the workflow, an optional input, include_gbff being set to true doesn't seem to lead to a gbff (gbk) file being used as the reference, nor being output to the Terra table. I'm not sure what this input is doing, but if we could use the gbk file as the input reference genome (this is feasible for snippy), we would get the gene annotations in the - snippy_concatenated_snps and snippy_shared_snps output files, which I think would be very helpful for users (rather than simply genome positions).
• Is there a reason that call_shared_variants = false? This seem like a very useful output so it would be good to have it by default
• The file name for the snippy_concatenated_variants output is just the tree_name input. Could we make this file name more descriptive and give it a .csv file extension so that it opens up in Excel properly?
• Note, the outputs on Terra are snippy_concatenated_variants and snippy_shared_variants_table, using the word "variants" rather than "snps" as mentioned in the PR.
• The output for snippy_variants_num_reads_aligned and other outputs coming from Snippy_Variants are simply listed as below. Whilst helpful for checking if there are any samples with insufficient quality to include in the tree, it doesn't help to identify which samples to exclude. Could each result perhaps be appended with the sample name?
• Similarly, snippy_variants_coverage_tsv outputs don't seem to be particularly helpful as a Snippy_Streamline output as the set-level output is just giving a list of GS URIs that aren't easy to access. Would it be worth removing this output for Snippy_Streamline?

[michellescribner]

@emmadoughty Thank you so much for all of your comments!! Answers to come below... When setting up the workflow, an optional input, include_gbff being set to true doesn't seem to lead to a gbff (gbk) file being used as the reference, nor being output to the Terra table. I'm not sure what this input is doing, but if we could use the gbk file as the input reference genome (this is feasible for snippy), we would get the gene annotations in the - snippy_concatenated_snps and snippy_shared_snps output files, which I think would be very helpful for users (rather than simply genome positions). -To my limited knowledge, include_gbff being set to true causes the ncbi_datasets task to download both the fasta file and the gbk file within Snippy Streamline, but the fasta is still the file passed on to the snippy variants task. I already provide gbk files as the input reference genome for when I directly provide a reference genome to the workflow, and can confirm that you wind up receiving the gene annotations in the concatenated snps file, so this is definitely already possible! I like your idea to make it possible for Snippy Streamline to use the gbk but not sure if we want to do that in this branch or create a new one.

[michellescribner]

Is there a reason that call_shared_variants = false? This seem like a very useful output so it would be good to have it by default

• good point, happy to make this change

The file name for the snippy_concatenated_variants output is just the tree_name input. Could we make this file name more descriptive and give it a .csv file extension so that it opens up in Excel properly?

• There is a reason for this particular madness. I used the existing cat_files task to concatenate the variants from multiple samples into one file, and this is the existing name of the file that is outputted by the task. I couldn't figure out how to add a ".csv" extension or other name to the file without altering the name of all files that would be concatenated with this task (some of which are not csvs). There might be a very easy solution to this, I will ask other folks!

[michellescribner]

Note, the outputs on Terra are snippy_concatenated_variants and snippy_shared_variants_table, using the word "variants" rather than "snps" as mentioned in the PR.

• This was intentional- snippy reports non-SNP mutations like indels as well as SNPs so I moved away from the SNP language. I should change the name of the PR.

[michellescribner]

The output for snippy_variants_num_reads_aligned and other outputs coming from Snippy_Variants are simply listed as below. Whilst helpful for checking if there are any samples with insufficient quality to include in the tree, it doesn't help to identify which samples to exclude. Could each result perhaps be appended with the sample name?

• The image that you showed with this question isn't the intended usage of the snippy_variants_num_reads_aligned column - it's meant to be examined after running snippy variants to determine which samples to include in a subsequent tree as shown below. I'm curious what you were looking at above?

[michellescribner]

Similarly, snippy_variants_coverage_tsv outputs don't seem to be particularly helpful as a Snippy_Streamline output as the set-level output is just giving a list of GS URIs that aren't easy to access. Would it be worth removing this output for Snippy_Streamline?

• I see your point, and I now understand your last question! These outputs should definitely be removed from Snippy_Streamline. They are helpful as outputs for snippy variants but the arrays are not useful.

[emmadoughty]

@michellescribner I have updated my comment above to reflect the scenarios to test. I hadn't realized that the QC outputs were not intended to be in Snippy_Streamline. That said, I think they would be incredibly useful in this workflow so I'm going to open a feature request issue to integrate this. I'm also going to open an enhancement issue to use the gbk file from reference seeker

[michellescribner]

Relaunched tests post modifications above

Tests launched on 23 C. auris specimens: TheiaEuk_Illumina_PE_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/93b97e30-18d6-407d-8c44-d33f45be4ad9 Snippy_Variants_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/e468f95d-dd3b-486c-992e-b82c8b4c2468 Snippy_Tree_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/7f6daae1-85bd-425a-be4e-cb2efb548dde kSNP3_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/f352cf52-79ff-44e3-ac86-8af9f3ec2824 -Verified that 3 new output columns are present and core tree remains identical to PHB v1.3 Snippy_Streamline_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/b366d926-27f4-4935-94e4-cbacb0156ab7

Other function tests to make sure nothing broke: TheiaProk_Illumina_SE_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/4fd029dc-d1a0-4333-8aa4-c06fe52991fc TheiaProk_FASTA_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/b42dd53d-e750-40ee-a184-e387925c6956 TheiaProk_ONT_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/4b30117b-c952-4545-871e-c101e44030cb Concatenate_Column_Content_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/43082c2a-b4a4-47a6-8fca-fdc089427b25 TheiaCoV_Illumina_PE_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/3e619181-cd6f-41c0-88fe-c66842c7280d Snippy_Streamline_PHB on some bacterial samples: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/5f4f0b53-5f24-4cdc-afb4-640dddb297f3

[michellescribner]

Retesting after silly mistake where I didn't update variable name:

TheiaEuk_Illumina_PE_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/9ace08ba-035a-4d32-bb1b-ce4421e32192 Snippy_Variants_PHB: https://app.terra.bio/#workspaces/theiagen-validations/Theiagen_Scribner_Sandbox/job_history/9aa4f317-0943-4609-aa46-d181ba8c8988

[sage-wright]

Code changes approved.

• [x] kSNP3 testing here
• [x] Snippy_Streamline testing here
  • tests both Snippy_Variants and Snippy_Tree
  • [x] TheiaEuk testing here
  • [x] Find_Shared_Variants testing here