Mapler/Snakefile at main · Nimauric/Mapler · 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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
###### Utility function ######

# Return a list containing the attribute "attribute" of each sample
# attribute = {name, read_path, short_reads_1, short_reads_2}
def get_samples(attribute) :
    return [sample[attribute] for sample in config["samples"]]

# Return the attribute "attribute" of an assembly with the name "sample"
def get_sample(attribute, wildcards):
    index = get_samples("name").index(wildcards.sample)
    return get_samples(attribute)[index]

# Return short reads either from the sample section or from global config
def get_short_read(attribute, wildcards):
    index = get_samples("name").index(wildcards.sample)
    sample = config["samples"][index]

    if attribute in sample:
        return sample[attribute]
    elif attribute in config:
        return config[attribute]
    else:
        raise ValueError(
            f"{attribute} not found neither in sample '{wildcards.sample}' nor in global config"
        )
def sample_has_long_reads(wildcards):
    return get_sample("read_path", wildcards) != "none"

def sample_has_short_reads(wildcards):
    index = get_samples("name").index(wildcards.sample)
    sample = config["samples"][index]

    if "short_reads_1" in sample and "short_reads_2" in sample:
        return True
    elif "short_reads_1" in config and "short_reads_2" in config:
        return True
    return False

def get_samples_with_long_reads():
    out = []
    for sample in config["samples"]:
        if "read_path" in sample and sample["read_path"] != "none":
            out.append(sample["name"])
    return out

def get_samples_with_short_reads():
    out = []
    for sample in config["samples"]:
        if "short_reads_1" in sample and "short_reads_2" in sample:
            out.append(sample["name"])
        elif "short_reads_1" in config and "short_reads_2" in config:
            out.append(sample["name"])
    return out

# Return the path to the reads of a fraction of an assembly
def get_read_path(wildcards) :
    if(wildcards.fraction == "full") :
        return get_sample("read_path", wildcards)
    return "outputs/" + wildcards.sample + "/" + wildcards.assembler + "/" + wildcards.fraction + "_reads.fastq"

# Return the path to the reads of all fractions of an assembly
def get_all_read_path(wildcards) :
    out = []
    for f in config["fractions"] :
        wildcards.fraction = f
        out.append(get_read_path(wildcards))
    return out

def get_reference_names() :
    return [os.path.splitext(os.path.basename(r))[0] for r in config["reference_genomes"]]

def get_reference(reference_name) :
    for ref in config["reference_genomes"] :
        if reference_name in ref :
            return ref
    return None


##### Sequencing technology (global) #####
TECH = config.get("technology", "hifi")

allowed_tech = ["hifi", "ont"]
if TECH not in allowed_tech:
    raise ValueError(f"technology must be one of {allowed_tech} (got: {TECH})")

# minimap2 preset for long reads mapping
LONGREAD_PRESET = "lr:hq" if TECH == "ont" else "map-hifi"


##### Additional rules #####
include : "rules/assembly.smk"
include : "rules/mapping.smk"
include : "rules/binning.smk"
include : "rules/bin_quality_analysis.smk"
include : "rules/read_quality_analysis.smk"
include : "rules/contig_quality_analysis.smk"

##### Improvements #####
"""
replace minimap2 by mapquick for hifi long reads mapping ?
replace kraken2 by sourmash for taxonomic assignation ?
"""


binnings = []
if(config["binning"] == True) :
    binnings += expand("{binner}_bins_reads_alignement", binner = config["binners"])
if(config["short_read_binning"] == True) :
    binnings += expand("{binner}_bins_short_reads_alignement", binner = config["binners"])
if(config["short_read_cobinning"] == True) :
    binnings += expand("{binner}_bins_cobinning_alignement", binner = config["binners"])
if(config["additional_reads_cobinning"] == True) :
    binnings += expand("{binner}_bins_additional_reads_cobinning_alignement", binner = config["binners"])


rule all :
    input :
        expand("outputs/{sample}/{assembler}/assembly.fasta", sample=get_samples("name"), assembler = config["assemblers"]),

        # Read quality analysis (fastqc, kraken2, kat)
        expand("outputs/{sample}/{assembler}/fastqc/{fraction}/fastqc_report.html", sample=get_samples("name"), assembler = config["assemblers"], fraction=config["fractions"])
            if(config["fastqc"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/kraken2/{fraction}/krona.html", sample=get_samples("name"), assembler = config["assemblers"], fraction=config["fractions"])
            if(config["kraken2"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/kat/{fraction}-stats.tsv", sample=get_samples("name"), assembler = config["assemblers"], fraction=config["fractions"])
            if(config["kat"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/kat/kat-plot.pdf", sample=get_samples("name"), assembler = config["assemblers"])
            if(config["kat"] == True and "mapped" in config["fractions"] and "unmapped" in config["fractions"]) else "Snakefile",

        # Contig quality analysis (read mapping, short read mapping, metaquast, reference mapping)
        expand("outputs/{sample}/{assembler}/reads_on_contigs_mapping_evaluation/report.txt", sample=get_samples("name"), assembler = config["assemblers"])
            if(config["read_mapping_evaluation"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/metaquast/report.txt", sample=get_samples("name"), assembler = config["assemblers"])
            if(config["metaquast"] == True and ("abundance_information" in config)) else "Snakefile",
        expand("outputs/{sample}/{assembler}/metaquast/results/summary/TSV/", sample=get_samples("name"), assembler = config["assemblers"])
            if(config["metaquast"] == True) else "Snakefile",
        expand("outputs/{sample}/long_reads_on_reference.{reference}.bam", sample=get_samples_with_long_reads(), reference=get_reference_names())
            if(config["reference_mapping_evaluation"] == True) else "Snakefile",
        expand("outputs/{sample}/short_reads_on_reference.{reference}.bam", sample=get_samples_with_short_reads(), reference=get_reference_names())
            if(config["reference_mapping_evaluation"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/short_reads_on_contigs.bam", sample=get_samples("name"), assembler=config["assemblers"])
            if(config["short_read_mapping_evaluation"] == True or config["short_read_binning"] == True or config["short_read_cobinning"] == True) else "Snakefile",

        # Bins quality analysis (checkm, separate read and contig quality analysis by bin quality)
        expand("outputs/{sample}/{assembler}/{binning}/checkm/checkm_report.txt", sample=get_samples("name"), assembler = config["assemblers"], binning=binnings)
            if(config["checkm"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/{binning}/checkm/checkm-plot.pdf", sample=get_samples("name"), assembler = config["assemblers"], binning=binnings)
            if(config["checkm"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/{binning}/gtdbtk/results/gtdbtk.bac120.summary.tsv", sample=get_samples("name"), assembler = config["assemblers"], binning=binnings)
            if(config["gtdbtk"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/{binning}/kraken2/bin.{target_bin}/krona.html", sample=get_samples("name"), assembler = config["assemblers"], binning=binnings, target_bin=config["kraken2_target_bins"])
            if(config["kraken2_on_bins"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/{binning}/read_contig_mapping_plot.pdf", sample=get_samples("name"), assembler = config["assemblers"], binning=binnings)
            if(config["checkm"] == True and config["read_mapping_evaluation"] == True) else "Snakefile",
        expand("outputs/{sample}/{assembler}/{binning}/read_contig_mapping.txt", sample=get_samples("name"), assembler = config["assemblers"], binning=binnings)
            if(config["checkm"] == True and config["read_mapping_evaluation"] == True) else "Snakefile",