zika/phylogenetic/rules/prepare_sequences.smk at 531904396345d7b8c67f58e9fc9c4ea35b02c41d · nextstrain/zika · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
"""
This part of the workflow prepares sequences for constructing the phylogenetic tree.

REQUIRED INPUTS:

    metadata_url    = url to metadata.tsv.zst
    sequences_url   = url to sequences.fasta.zst
    reference   = path to reference sequence or genbank

OUTPUTS:

    prepared_sequences = results/aligned.fasta

This part of the workflow usually includes the following steps:

    - augur index
    - augur filter
    - augur align
    - augur mask

See Augur's usage docs for these commands for more details.
"""

rule filter:
    """
    Filtering to
      - {params.sequences_per_group} sequence(s) per {params.group_by!s}
      - from {params.min_date} onwards
      - excluding strains in {input.exclude}
      - minimum genome length of {params.min_length} (50% of Zika virus genome)
    """
    input:
        sequences = input_sequences,
        metadata = input_metadata,
        exclude = resolve_config_path(config["exclude"]),
    output:
        sequences = "results/filtered.fasta"
    params:
        group_by = as_list(config["filter"]["group_by"]),
        sequences_per_group = config["filter"]["sequences_per_group"],
        min_date = config["filter"]["min_date"],
        min_length = config["filter"]["min_length"],
        strain_id = config.get("strain_id_field", "strain"),
    log:
        "logs/filter.txt",
    benchmark:
        "benchmarks/filter.txt"
    shell:
        r"""
        exec &> >(tee {log:q})

        augur filter \
            --sequences {input.sequences:q} \
            --metadata {input.metadata:q} \
            --metadata-id-columns {params.strain_id:q} \
            --exclude {input.exclude:q} \
            --output {output.sequences:q} \
            --group-by {params.group_by:q} \
            --sequences-per-group {params.sequences_per_group:q} \
            --min-date {params.min_date:q} \
            --min-length {params.min_length:q}
        """


# rule upload_filter:
#     """TESTING ONLY TODO XXX REMOVE"""
#     input:
#         sequences = input_sequences,
#         metadata = input_metadata,
#         exclude = resolve_config_path(config["exclude"]),
#     output:
#         sequences = path_or_url("s3://nextstrain-scratch/zika-pr-89/filtered.fasta")
#     params:
#         group_by = as_list(config["filter"]["group_by"]),
#         sequences_per_group = config["filter"]["sequences_per_group"],
#         min_date = config["filter"]["min_date"],
#         min_length = config["filter"]["min_length"],
#         strain_id = config.get("strain_id_field", "strain"),
#     log:
#         "logs/filter.txt",
#     benchmark:
#         "benchmarks/filter.txt"
#     shell:
#         r"""
#         exec &> >(tee {log:q})

#         augur filter \
#             --sequences {input.sequences:q} \
#             --metadata {input.metadata:q} \
#             --metadata-id-columns {params.strain_id:q} \
#             --exclude {input.exclude:q} \
#             --output-sequences {output.sequences:q} \
#             --group-by {params.group_by:q} \
#             --sequences-per-group {params.sequences_per_group:q} \
#             --min-date {params.min_date:q} \
#             --min-length {params.min_length:q}
#         """


rule align:
    """
    Aligning sequences to {input.reference}
      - filling gaps with N
    """
    input:
        sequences = "results/filtered.fasta",
        reference = resolve_config_path(config["reference"]),
    output:
        alignment = "results/aligned.fasta"
    log:
        "logs/align.txt",
    benchmark:
        "benchmarks/align.txt"
    shell:
        r"""
        exec &> >(tee {log:q})

        augur align \
            --sequences {input.sequences:q} \
            --reference-sequence {input.reference:q} \
            --output {output.alignment:q} \
            --fill-gaps \
            --remove-reference
        """