diff --git a/Snakefile.genome b/Snakefile.genome new file mode 100644 index 0000000..1ec918f --- /dev/null +++ b/Snakefile.genome @@ -0,0 +1,260 @@ +import json +import glob + +SEGMENTS = ["ha", "np", "pb2", "pb1", "pa", "np", "mp", "ns"] # order important! +SUBTYPES = ['h5n1'] + + +rule all: + input: expand("results/genome/aligned_{subtype}_{segment}.fasta", subtype=SUBTYPES, segment=SEGMENTS) + + +rule files: + params: + # dropped_strains = "config/dropped_strains_{subtype}.txt", + # include_strains = "config/include_strains_{subtype}_{time}.txt", + reference = "config/reference_{subtype}_{segment}.gb", + include = "config/include_genome-strains_{subtype}.txt", + colors = "config/colors_{subtype}.tsv", + lat_longs = "config/lat_longs_{subtype}.tsv", + auspice_config = "config/auspice_config_{subtype}.json", + # clades_file = "clade-labeling/{subtype}-clades.tsv", + description = "config/description.md" + +files = rules.files.params + +rule filter: + input: + sequences = "results/sequences_{subtype}_{segment}.fasta", + metadata = "results/metadata_{subtype}_{segment}.tsv", + include = files.include, + output: + sequences = "results/genome/sequences_{subtype}_{segment}.fasta", + metadata = "results/genome/metadata_{subtype}_{segment}.tsv", + shell: + """ + augur filter \ + --sequences {input.sequences} \ + --metadata {input.metadata} \ + --exclude-all --include {input.include} \ + --output-sequences {output.sequences} \ + --output-metadata {output.metadata} + """ + +rule align: + input: + sequences ="results/genome/sequences_{subtype}_{segment}.fasta", + reference = files.reference, + output: + alignment = "results/genome/aligned_{subtype}_{segment}.fasta", + shell: + """ + augur align \ + --sequences {input.sequences} \ + --reference-sequence {input.reference} \ + --output {output.alignment} \ + --remove-reference \ + --nthreads 1 + """ + +rule join_sequences: + input: + alignment = expand("results/genome/aligned_{{subtype}}_{segment}.fasta", segment=SEGMENTS), + output: + alignment = "results/genome/genome_{subtype}.fasta", + shell: + """ + python scripts/join-segments.py \ + --segments {input.alignment} \ + --output {output.alignment} + """ + +# NOTE - I'm not joining metadata, I'm relying on everything being there in the starting HA metadata... + +# NOTE - from now on the rules mirror the main snakefile... +rule tree: + message: "Building tree" + input: + alignment = "results/genome/genome_{subtype}.fasta", + output: + tree = "results/genome/tree-raw_{subtype}.nwk", + params: + method = "iqtree" + shell: + """ + augur tree \ + --alignment {input.alignment} \ + --output {output.tree} \ + --method {params.method} \ + --nthreads 1 + """ + +def clock_rate(w): + """ + MVP!!!!!!!!! + """ + clock_rates_h5n1 = { + 'pb2': 0.00287, + 'pb1': 0.00264, + 'pa': 0.00248, + 'ha': 0.00455, + 'np': 0.00252, + 'na': 0.00349, + 'mp': 0.00191, + 'ns': 0.00249 + } + lengths = { + 'pb2': 2341, + 'pb1': 2341, + 'pa': 2233, + 'ha': 1760, + 'np': 1565, + 'na': 1458, + 'mp': 1027, + 'ns': 865 + } + mean = sum([cr * lengths[seg] for seg,cr in clock_rates_h5n1.items()])/sum(lengths.values()) + stdev = mean/2 + return f"--clock-rate {mean} --clock-std-dev {stdev}" + +rule refine: + message: + """ + Refining tree + - estimate timetree + - use {params.coalescent} coalescent timescale + - estimate {params.date_inference} node dates + """ + input: + tree = "results/genome/tree-raw_{subtype}.nwk", + alignment = "results/genome/genome_{subtype}.fasta", + metadata = "results/metadata_{subtype}_ha.tsv", # assuming all salient info is in HA + output: + tree = "results/genome/tree_{subtype}.nwk", + node_data = "results/branch-lengths_{subtype}.json" + params: + coalescent = "const", + date_inference = "marginal", + clock_rate = clock_rate, + shell: + """ + augur refine \ + --tree {input.tree} \ + --alignment {input.alignment} \ + --metadata {input.metadata} \ + --output-tree {output.tree} \ + --output-node-data {output.node_data} \ + --timetree \ + --coalescent {params.coalescent} \ + --date-confidence \ + --date-inference {params.date_inference} \ + {params.clock_rate} \ + """ + +rule ancestral: + input: + tree = "results/genome/tree_{subtype}.nwk", + alignment = "results/genome/genome_{subtype}.fasta", + output: + node_data = "results/genome/nt-muts_{subtype}.json" + params: + inference = "joint" + shell: + """ + augur ancestral \ + --tree {input.tree} \ + --alignment {input.alignment} \ + --output-node-data {output.node_data} \ + --inference {params.inference}\ + --keep-ambiguous + """ + +# TODO XXX translate needs joined genbank + +# rule translate: +# message: "Translating amino acid sequences" +# input: +# tree = "results/genome/tree_{subtype}.nwk", +# node_data = "results/genome/nt-muts_{subtype}.json" +# reference = files.reference +# output: +# node_data = "results/aa-muts_{subtype}_{segment}_{time}.json" +# shell: +# """ +# augur translate \ +# --tree {input.tree} \ +# --ancestral-sequences {input.node_data} \ +# --reference-sequence {input.reference} \ +# --output {output.node_data} +# """ + + +def traits_columns(w): + traits = {'h5nx':'region','h5n1': 'region country', 'h7n9': 'country division', 'h9n2': 'region country'} + return traits[w.subtype] + +rule traits: + message: "Inferring ancestral traits for {params.columns!s}" + input: + tree = "results/genome/tree_{subtype}.nwk", + metadata = "results/metadata_{subtype}_ha.tsv", + output: + node_data = "results/genome/traits_{subtype}.json", + params: + columns = traits_columns, + shell: + """ + augur traits \ + --tree {input.tree} \ + --metadata {input.metadata} \ + --output {output.node_data} \ + --columns {params.columns} \ + --confidence + """ + +rule cleavage_site: + input: + alignment = "results/genome/aligned_{subtype}_ha.fasta", + output: + cleavage_site_annotations = "results/genome/cleavage-site_{subtype}_ha.json", + cleavage_site_sequences = "results/genome/cleavage-site-sequences_{subtype}_ha.json" + shell: + """ + python scripts/annotate-ha-cleavage-site.py \ + --alignment {input.alignment} \ + --furin_site_motif {output.cleavage_site_annotations} \ + --cleavage_site_sequence {output.cleavage_site_sequences} + """ + +rule export: + message: "Exporting data files for for auspice" + input: + tree = "results/genome/tree_{subtype}.nwk", + metadata = "results/metadata_{subtype}_ha.tsv", + node_data = [ + rules.refine.output.node_data, + rules.traits.output.node_data, + rules.ancestral.output.node_data, + # rules.translate.output.node_data, + rules.cleavage_site.output.cleavage_site_annotations, + rules.cleavage_site.output.cleavage_site_sequences, + ], + colors = files.colors, + lat_longs = files.lat_longs, + auspice_config = files.auspice_config, + description = files.description, + output: + auspice_json = "auspice/avian-flu_{subtype}_genome.json" + shell: + """ + augur export v2 \ + --tree {input.tree} \ + --metadata {input.metadata} \ + --node-data {input.node_data} \ + --colors {input.colors} \ + --lat-longs {input.lat_longs} \ + --auspice-config {input.auspice_config} \ + --description {input.description} \ + --include-root-sequence-inline \ + --output {output.auspice_json} + """ diff --git a/config/include_genome-strains_h5n1.txt b/config/include_genome-strains_h5n1.txt new file mode 100644 index 0000000..fe64a78 --- /dev/null +++ b/config/include_genome-strains_h5n1.txt @@ -0,0 +1,257 @@ +A/Texas/37/2024 +A/Goose/USA/23-038138-001-original/2024 +A/skunk/NewMexico/24006483001original/2024 +A/Skunk/USA/24-006483-001-original/2024 +A/Cattle/USA/24-009110-003-original-repeat/2024 +A/Cattle/USA/24-010193-006-original/2024 +A/Cattle/USA/24-009110-023-original/2024 +A/Chicken/USA/24-009654-001-original/2024 +A/Cattle/USA/24-010193-007-original/2024 +A/Cattle/USA/24-009110-019-original/2024 +A/Cattle/USA/24-010192-006-original/2024 +A/Cattle/USA/24-010193-003-original/2024 +A/Cattle/USA/24-010330-012-original-300/2024 +A/Grackle/USA/24-008765-002-original/2024 +A/Cattle/USA/24-010330-013-original-300/2024 +A/Cattle/USA/24-010512-001-original/2024 +A/Cattle/USA/24-010330-002-original-300/2024 +A/Cattle/USA/24-010330-005-original/2024 +A/Cattle/USA/24-010330-007-original/2024 +A/Cattle/USA/24-010330-003-original/2024 +A/Cattle/USA/24-010330-009-original-300/2024 +A/Cattle/USA/24-010330-008-original/2024 +A/Cattle/USA/24-010330-006-original-300/2024 +A/Cattle/USA/24-010330-004-original-300/2024 +A/Cattle/USA/24-010514-002-original/2024 +A/Cattle/USA/24-010327-003-original/2024 +A/Cattle/USA/24-010071-001-original/2024 +A/Cattle/USA/24-009776-002-original/2024 +A/Cattle/USA/24-009110-020-original/2024 +A/Cattle/USA/24-009306-004-original/2024 +A/Cattle/USA/24-010513-002-original/2024 +A/Chicken/USA/24-010482-001-original/2024 +A/Cattle/USA/24-009495-001-original/2024 +A/Cattle/USA/24-009109-002-original/2024 +A/Cattle/USA/24-010195-001-original/2024 +A/Cattle/USA/24-010071-008-original/2024 +A/Cattle/USA/24-009109-003-original/2024 +A/Cattle/USA/24-010513-006-original/2024 +A/Cattle/USA/24-009497-012-original/2024 +A/Cattle/USA/24-009497-004-original/2024 +A/Chicken/USA/24-007264-001-original/2024 +A/Chicken/USA/24-007264-003-original/2024 +A/Cattle/USA/24-009309-001-original/2024 +A/Cattle/USA/24-009309-002-original/2024 +A/Cattle/USA/24-009027-005-original-MTM/2024 +A/Cattle/USA/24-009109-001-original/2024 +A/Cattle/USA/24-010303-001-original-300/2024 +A/Cattle/USA/24-009027-003-original-repeat/2024 +A/Cattle/USA/24-008749-006-v/2024 +A/Cattle/USA/24-008660-002-original/2024 +A/dairycattle/Texas/24008749006original/2024 +A/Cattle/USA/24-008749-006-original/2024 +A/Cattle/USA/24-009586-003-original/2024 +A/Cattle/USA/24-009586-009-original/2024 +A/Cattle/USA/24-009586-007-original/2024 +A/Cattle/USA/24-009586-006-original/2024 +A/Cattle/USA/24-009491-005-original/2024 +A/Cattle/USA/24-009491-001-original/2024 +A/Cattle/USA/24-009491-003-original/2024 +A/Cattle/USA/24-009491-006-original/2024 +A/feline/USA/24009116005original/2024 +A/feline/USA/24009116002original/2024 +A/Cat/USA/24-009116-005-original/2024 +A/Cat/USA/24-009116-004-original/2024 +A/Cat/USA/24-009116-002-original/2024 +A/feline/USA/24009116004original/2024 +A/Cattle/USA/24-009775-001-original/2024 +A/Cattle/USA/24-009367-009-original/2024 +A/Cattle/USA/24-009367-010-original/2024 +A/Cattle/USA/24-009367-014-original/2024 +A/Cattle/USA/24-009367-003-original/2024 +A/Cattle/USA/24-009367-005-original/2024 +A/Cattle/USA/24-009367-013-original/2024 +A/Cattle/USA/24-009028-019-original/2024 +A/Cattle/USA/24-009028-001-original/2024 +A/Cattle/USA/24-009499-001-original/2024 +A/Cattle/USA/24-009028-005-original-repeat/2024 +A/Cattle/USA/24-009367-011-original/2024 +A/Cattle/USA/24-009028-003-original-MTM/2024 +A/Cattle/USA/24-009028-007-original/2024 +A/Cattle/USA/24-009028-008-original/2024 +A/dairycattle/Texas/24008749002v/2024 +A/Cattle/USA/24-008749-002-v/2024 +A/Cattle/USA/24-008749-004-original/2024 +A/Cattle/USA/24-008749-005-original/2024 +A/Cattle/USA/24-008749-004-v/2024 +A/Cattle/USA/24-008749-007-original/2024 +A/Cattle/USA/24-008749-008-original-repeat/2024 +A/Cattle/USA/24-008749-003-original/2024 +A/Cattle/USA/24-008660-001-original/2024 +A/Cattle/USA/24-008749-001-original/2024 +A/dairycattle/Texas/24008749003original/2024 +A/dairycattle/Texas/24008749007original/2024 +A/dairycattle/Texas/24008749005original/2024 +A/dairycattle/Texas/24008749001original/2024 +A/dairycattle/Texas/24008749004original/2024 +A/Cattle/USA/24-010354-009-original-300/2024 +A/Cattle/USA/24-010354-011-original/2024 +A/Cattle/USA/24-010354-001-original/2024 +A/Cattle/USA/24-009497-007-original/2024 +A/Cattle/USA/24-010354-018-original-300/2024 +A/Cattle/USA/24-010354-005-original/2024 +A/Cattle/USA/24-010354-016-original-300/2024 +A/Cattle/USA/24-010354-013-original-300/2024 +A/Cattle/USA/24-010354-002-original/2024 +A/Cattle/USA/24-010354-004-original-300/2024 +A/Cattle/USA/24-010354-008-original/2024 +A/Cattle/USA/24-010354-020-original/2024 +A/Cattle/USA/24-010354-014-original-300/2024 +A/Cattle/USA/24-010354-017-original-300/2024 +A/Cattle/USA/24-010354-006-original-300/2024 +A/Cattle/USA/24-010354-012-original/2024 +A/Cattle/USA/24-010354-015-original-300/2024 +A/Cattle/USA/24-010354-010-original-300/2024 +A/Cattle/USA/24-010354-019-original-300/2024 +A/Cattle/USA/24-010354-003-original/2024 +A/Cattle/USA/24-008767-001-v/2024 +A/Cattle/USA/24-009110-018-original/2024 +A/Cattle/USA/24-009110-026-original/2024 +A/Cattle/USA/24-009110-015-original/2024 +A/Cattle/USA/24-009110-021-original/2024 +A/Cattle/USA/24-009110-027-original/2024 +A/Cattle/USA/24-009110-002-original/2024 +A/Cattle/USA/24-009088-002-original/2024 +A/Cattle/USA/24-009110-001-original/2024 +A/Cattle/USA/24-009290-005-original/2024 +A/Cattle/USA/24-009110-025-original/2024 +A/Cattle/USA/24-009110-005-original/2024 +A/Cattle/USA/24-009110-016-original/2024 +A/Cattle/USA/24-009088-001-original/2024 +A/Cattle/USA/24-009110-006-original/2024 +A/Cattle/USA/24-009110-012-original/2024 +A/Cattle/USA/24-009110-004-original/2024 +A/Cattle/USA/24-009110-024-original/2024 +A/Cattle/USA/24-009290-006-original/2024 +A/Cattle/USA/24-009110-022-original/2024 +A/Cattle/USA/24-009110-010-original/2024 +A/Cattle/USA/24-009110-009-original/2024 +A/Cattle/USA/24-009110-007-original/2024 +A/Cattle/USA/24-009110-011-original/2024 +A/Cattle/USA/24-009110-014-original/2024 +A/Cat/USA/24-008850-002-original/2024 +A/Raccoon/USA/24-009496-002-original/2024 +A/Chicken/USA/24-010308-007-original/2024 +A/Cattle/USA/24-009497-010-original/2024 +A/Cattle/USA/24-009586-002-original/2024 +A/feline/USA/24008764002original/2024 +A/Cattle/USA/24-009309-004-original/2024 +A/feline/USA/24008850002original/2024 +A/Cattle/USA/24-009110-017-original/2024 +A/Cattle/USA/24-009109-008-original/2024 +A/Chicken/USA/24-008355-002-original/2024 +A/Cattle/USA/24-009367-018-original/2024 +A/Cattle/USA/24-009586-001-original/2024 +A/Cattle/USA/24-009367-012-original/2024 +A/Cattle/USA/24-009108-004-original/2024 +A/Cattle/USA/24-010071-002-original/2024 +A/Cattle/USA/24-009586-008-original/2024 +A/Cattle/USA/24-009495-008-original/2024 +A/Cattle/USA/24-010071-003-original/2024 +A/Cattle/USA/24-009367-001-original/2024 +A/Cattle/USA/24-009497-014-original/2024 +A/Cattle/USA/24-009110-008-original-repeat/2024 +A/Cattle/USA/24-010071-009-original/2024 +A/Cattle/USA/24-009491-007-original/2024 +A/Cattle/USA/24-008766-001-original/2024 +A/Cattle/USA/24-010303-002-original/2024 +A/Cattle/USA/24-009028-009-original/2024 +A/Cattle/USA/24-009497-005-original/2024 +A/Cattle/USA/24-009776-004-original/2024 +A/Cattle/USA/24-009308-004-original/2024 +A/Cattle/USA/24-010513-007-original/2024 +A/Cattle/USA/24-009109-006-original/2024 +A/Cattle/USA/24-009367-004-original/2024 +A/Chicken/USA/24-010308-006-original-300/2024 +A/Cattle/USA/24-009776-015-original/2024 +A/Cattle/USA/24-009310-010-original/2024 +A/feline/USA/24009311004original/2024 +A/Cattle/USA/24-009497-003-original/2024 +A/Cattle/USA/24-009586-005-original/2024 +A/Grackle/USA/24-008356-003-original/2024 +A/Chicken/USA/24-010308-012-original/2024 +A/Cattle/USA/24-009310-007-original-repeat/2024 +A/Cattle/USA/24-009367-002-original/2024 +A/Cattle/USA/24-009029-001-original/2024 +A/Cattle/USA/24-010331-001-original/2024 +A/Cattle/USA/24-009495-010-original/2024 +A/Cattle/USA/24-009367-023-original/2024 +A/Cattle/USA/24-009109-007-original/2024 +A/Cattle/USA/24-009027-010-original-repeat/2024 +A/Chicken/USA/24-009582-001-original/2024 +A/Cattle/USA/24-009367-016-original/2024 +A/Cattle/USA/24-009310-006-original/2024 +A/Cattle/USA/24-009308-006-original/2024 +A/Cattle/USA/24-009108-007-original/2024 +A/Cattle/USA/24-009586-010-original/2024 +A/feline/USA/24008850001original/2024 +A/Cat/USA/24-009311-004-original/2024 +A/Cattle/USA/24-010513-004-original/2024 +A/Cat/USA/24-009311-006-original/2024 +A/Cattle/USA/24-009087-001-original/2024 +A/Cattle/USA/24-009367-007-original/2024 +A/Cattle/USA/24-009027-002-original/2024 +A/Cattle/USA/24-009108-006-original/2024 +A/Cattle/USA/24-010071-006-original/2024 +A/Blackbird/USA/24-008357-001-original/2024 +A/Chicken/USA/24-008355-003-original/2024 +A/Cattle/USA/24-009491-002-original/2024 +A/Chicken/USA/24-009582-002-original/2024 +A/Cattle/USA/24-010513-003-original/2024 +A/Cattle/USA/24-009108-001-original/2024 +A/Chicken/USA/24-008355-004-original/2024 +A/Cattle/USA/24-009108-003-original/2024 +A/Chicken/USA/24-008355-005-original/2024 +A/Cattle/USA/24-010194-001-original/2024 +A/Cattle/USA/24-009108-005-original/2024 +A/Cattle/USA/24-008766-001-v/2024 +A/Cattle/USA/24-009497-009-original/2024 +A/Cattle/USA/24-010071-005-original/2024 +A/Cat/USA/24-008764-002-original/2024 +A/feline/USA/24009311006original/2024 +A/Chicken/USA/24-009582-003-original/2024 +A/Cattle/USA/24-009581-002-original/2024 +A/Cattle/USA/24-009027-001-v/2024 +A/Cattle/USA/24-009497-008-original/2024 +A/Cattle/USA/24-010071-007-original/2024 +A/Chicken/USA/24-008355-006-original/2024 +A/Cattle/USA/24-009306-003-original/2024 +A/Cattle/USA/24-009108-002-original/2024 +A/feline/USA/24008764001original/2024 +A/Cattle/USA/24-010330-011-original-300/2024 +A/Cattle/USA/24-009027-002-v/2024 +A/Cattle/USA/24-009308-008-original/2024 +A/Cattle/USA/24-009109-004-original/2024 +A/Cattle/USA/24-009497-006-original/2024 +A/Cattle/USA/24-010071-010-original/2024 +A/Cattle/USA/24-009290-001-original/2024 +A/Cattle/USA/24-009776-001-original/2024 +A/Grackle/USA/24-008356-001-original/2024 +A/Cattle/USA/24-009027-004-original-repeat/2024 +A/Blackbird/USA/24-008354-001-original/2024 +A/Chicken/USA/24-010308-001-original-300/2024 +A/Cattle/USA/24-009495-009-original/2024 +A/Cattle/USA/24-009495-012-original/2024 +A/Cattle/USA/24-009109-005-original/2024 +A/Cattle/USA/24-009497-011-original/2024 +A/Cattle/USA/24-009495-011-original/2024 +A/Cat/USA/24-008764-001-original/2024 +A/Cattle/USA/24-009495-002-original/2024 +A/Cat/USA/24-008850-001-original/2024 +A/Cattle/USA/24-009310-005-original/2024 +A/Chicken/USA/24-008355-001-original/2024 +A/Cattle/USA/24-009495-003-original/2024 +A/blackbird/Texas/24008354001original/2024 +A/blackbird/Texas/24008357001original/2024 +A/commongrackle/Texas/24008356003original/2024 +A/commongrackle/Texas/24008356001original/2024 \ No newline at end of file diff --git a/scripts/join-segments.py b/scripts/join-segments.py new file mode 100644 index 0000000..ba0a7b4 --- /dev/null +++ b/scripts/join-segments.py @@ -0,0 +1,26 @@ + +from Bio import SeqIO +from collections import defaultdict +import argparse + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument('--segments', type = str, required = True, nargs='+', help = "per-segment alignments") + parser.add_argument('--output', type = str, required = True, help = "output whole genome alignment") + args = parser.parse_args() + + records = {} + strain_counts = defaultdict(int) + for fname in args.segments: + records[fname] = {record.name:str(record.seq) for record in SeqIO.parse(fname, 'fasta')} + for key in records[fname]: strain_counts[key]+=1 + print(f"{fname}: parsed {len(records[fname].keys())} sequences") + + with open(args.output, 'w') as fh: + print("writing genome to ", args.output) + for name,count in strain_counts.items(): + if count!=len(args.segments): + print(f"Excluding {name} as it only appears in {count} segments") + continue + genome = "".join([records[seg][name] for seg in args.segments]) + print(f">{name}\n{genome}", file=fh)