-
Notifications
You must be signed in to change notification settings - Fork 0
/
smkalign.py
82 lines (70 loc) · 2 KB
/
smkalign.py
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
import os
# Load the configuration file
configfile: "config.yaml"
# Extract variables from the config file
reference = config["reference"]
fastq_dir = config["fastq_dir"]
# Detect all .fastq.gz files in the specified directory
input_files = glob_wildcards(f"{fastq_dir}/{{sample}}.fastq.gz").sample
rule all:
input:
expand("results/{sample}.bam.bai", sample=input_files),
expand("results/{sample}.fasta", sample=input_files),
expand("results/{sample}_depth.txt", sample=input_files) # Include depth files in all rule
rule cutadapt:
input:
fastq = f"{fastq_dir}/{{sample}}.fastq.gz"
output:
cleaned_fastq = "results/cleaned_{sample}.fastq.gz"
params:
adapter = "CTGTCTCTTATACACATCT"
shell:
"""
cutadapt -q 20 -a {params.adapter} --minimum-length 50 -o {output.cleaned_fastq} {input.fastq}
"""
rule bwa_mem2:
input:
reference = reference,
cleaned_fastq = "results/cleaned_{sample}.fastq.gz"
output:
sam = "results/{sample}.sam"
shell:
"""
bwa-mem2 mem -t 16 {input.reference} {input.cleaned_fastq} > {output.sam}
"""
rule samtools_sort:
input:
sam = "results/{sample}.sam"
output:
bam = "results/{sample}.bam"
shell:
"""
samtools view -bS {input.sam} | samtools sort -o {output.bam}
"""
rule samtools_index:
input:
bam = "results/{sample}.bam"
output:
index = "results/{sample}.bam.bai"
shell:
"""
samtools index {input.bam}
"""
rule samtools_consensus:
input:
bam = "results/{sample}.bam"
output:
fasta = "results/{sample}.fasta"
shell:
"""
samtools consensus -f fasta {input.bam} -o {output.fasta}
"""
rule samtools_depth:
input:
bam = "results/{sample}.bam"
output:
depth_txt = "results/{sample}_depth.txt"
shell:
"""
samtools depth {input.bam} > {output.depth_txt}
"""