-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathSnakefile
140 lines (119 loc) · 5.52 KB
/
Snakefile
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
"""
Juno_assembly pipeline
Authors: Alejandra Hernandez-Segura, Robert Verhagen, Ernst Hamer, Dennis Schmitz, Diogo Borst, Tom van Wijk, Maaike van der Beld
Organization: Rijksinstituut voor Volksgezondheid en Milieu (RIVM)
Department: Infektieziekteonderzoek, Diagnostiek en Laboratorium Surveillance (IDS), Bacteriologie (BPD)
Date: 26-08-2021
Documentation: https://rivm-bioinformatics.github.io/ids_bacteriology_man/juno-assembly.html
"""
#################################################################################
##### Import config file, sample_sheet and set output folder names #####
#################################################################################
import os
import yaml
#################################################################################
##### Load samplesheet, load genus dict and define output directory #####
#################################################################################
# SAMPLES is a dict with sample in the form:
# sample > read number > file. E.g.: SAMPLES["sample_1"]["R1"] = "x_R1.gz"
sample_sheet = config["sample_sheet"]
SAMPLES = {}
with open(sample_sheet) as sample_sheet_file:
SAMPLES = yaml.safe_load(sample_sheet_file)
OUT = config["out"]
IN = config["input_dir"]
for param in ["threads", "mem_gb"]:
for k in config[param]:
config[param][k] = int(config[param][k])
# @################################################################################
# @#### Processes #####
# @################################################################################
#############################################################################
##### Data quality control and cleaning #####
#############################################################################
include: "bin/rules/fastqc_raw_data.smk"
include: "bin/rules/clean_fastq.smk"
include: "bin/rules/fastqc_clean_data.smk"
include: "bin/rules/subsample_fastq.smk"
#############################################################################
##### De novo assembly #####
#############################################################################
include: "bin/rules/de_novo_assembly.smk"
#############################################################################
##### Species identification #####
#############################################################################
include: "bin/rules/identify_species.smk"
#############################################################################
##### Scaffold analyses: QUAST, CheckM, picard, bbmap and QC-metrics #####
#############################################################################
include: "bin/rules/run_quast.smk"
include: "bin/rules/run_checkm.smk"
include: "bin/rules/parse_checkm.smk"
include: "bin/rules/picard_insert_size.smk"
include: "bin/rules/generate_contig_metrics.smk"
include: "bin/rules/parse_bbtools.smk"
include: "bin/rules/parse_bbtools_summary.smk"
include: "bin/rules/multiqc.smk"
include: "bin/rules/create_juno_qc.smk"
# @################################################################################
# @#### The `onstart` checker codeblock #####
# @################################################################################
onstart:
print("Checking if all specified files are accessible...")
important_files = [
config["sample_sheet"],
"files/trimmomatic_0.36_adapters_lists/NexteraPE-PE.fa",
"files/accepted_genera_checkm.txt",
"files/multiqc_config.yaml",
]
for filename in important_files:
if not os.path.exists(filename):
raise FileNotFoundError(filename)
#################################################################################
##### Specify final output: #####
#################################################################################
localrules:
all,
select_genus_checkm,
rule all:
input:
expand(
OUT + "/qc_raw_fastq/{sample}_{read}_fastqc.zip",
sample=SAMPLES,
read=["R1", "R2"],
),
expand(
OUT + "/clean_fastq/{sample}_{read}.fastq.gz",
sample=SAMPLES,
read=["pR1", "pR2"],
),
expand(
OUT + "/qc_clean_fastq/{sample}_{read}_fastqc.zip",
sample=SAMPLES,
read=["pR1", "pR2"],
),
expand(OUT + "/de_novo_assembly/{sample}/scaffolds.fasta", sample=SAMPLES),
expand(
OUT
+ "/qc_de_novo_assembly/checkm/per_sample/{sample}/checkm_{sample}.tsv",
sample=SAMPLES,
),
OUT + "/qc_de_novo_assembly/quast/report.tsv",
expand(
OUT
+ "/qc_de_novo_assembly/bbtools_scaffolds/per_sample/{sample}_MinLenFiltSummary.tsv",
sample=SAMPLES,
),
OUT + "/qc_de_novo_assembly/bbtools_scaffolds/bbtools_scaffolds.tsv",
OUT + "/qc_de_novo_assembly/bbtools_scaffolds/bbtools_summary_report.tsv",
expand(
OUT + "/identify_species/reads/{sample}/{sample}_species_content.txt",
sample=SAMPLES,
),
expand(
OUT + "/identify_species/reads/{sample}/{sample}_bracken_species.kreport2",
sample=SAMPLES,
),
OUT + "/identify_species/top1_species_multireport.csv",
OUT + "/multiqc/multiqc.html",
OUT + "/Juno_assembly_QC_report/QC_report.xlsx",