Commit 1579db1f authored by Slaheddine Kastalli's avatar Slaheddine Kastalli
Browse files

modifications

parent 7cc7db32
...@@ -4,7 +4,7 @@ derepF <- derepFastq(snakemake@input[[1]], verbose = TRUE) ...@@ -4,7 +4,7 @@ derepF <- derepFastq(snakemake@input[[1]], verbose = TRUE)
derepR <- derepFastq(snakemake@input[[2]], verbose = TRUE) derepR <- derepFastq(snakemake@input[[2]], verbose = TRUE)
errF <- learnErrors(snakemake@input[[1]], multithread = snakemake@threads, verbose = TRUE) errF <- learnErrors(snakemake@input[[1]], multithread = snakemake@threads, verbose = TRUE)
errR <- learnErrors(snakemake@input[[2]], multithread = snakemake@threads, verbose = TRUE) errR <- learnErrors(snakemake@input[[2]], multithread = snakemake@threads, verbose = TRUE)
dadaFs <- dada(derep, err = errF, multithread = snakemake@threads, verbose = TRUE) dadaFs <- dada(derepF, err = errF, multithread = snakemake@threads, verbose = TRUE)
dadaRs <- dada(derep, err = errR, multithread = snakemake@threads, verbose = TRUE) dadaRs <- dada(derepR, err = errR, multithread = snakemake@threads, verbose = TRUE)
mergers <- mergePairs(dadaFs, dadaRs, verbose=TRUE) mergers <- mergePairs(dadaFs, derepF, dadaRs, derepR, verbose=TRUE)
saveRDS(dada, snakemake@output[[1]]) saveRDS(mergers, snakemake@output[[1]])
rule dada2: rule dada2:
input: input:
filt = "work/filter/{sample}.sra_1.fastq.gz", filt = "work/filter/{sample}_1.fastq.gz",
filtrev = "work/filter/{sample}.sra_2.fastq.gz" filtrev = "work/filter/{sample}_2.fastq.gz"
output: output:
"work/dada/{sample}.rds" "work/dada/{sample}.rds"
threads: threads:
...@@ -25,7 +25,7 @@ rule makeSequenceTable: ...@@ -25,7 +25,7 @@ rule makeSequenceTable:
rule metrics: rule metrics:
input: input:
rds = "work/dada/seqtab.rds", rds = "work/dada/seqtab.rds",
fastqc = expand("work/fastqc/{sample}_fastqc.zip", sample=SAMPLES) fastqc = expand("work/fastqc/{sample}_{unit}_fastqc.zip", sample=SAMPLES, unit=unit)
output: output:
tsv = "work/dada/metrics.tsv" tsv = "work/dada/metrics.tsv"
script: script:
......
library(dada2) library(dada2)
filterAndTrim(snakemake@input[[1]],snakemake@input[[2]] snakemake@output[[1]], maxN=0, rm.phix=TRUE, compress=TRUE, verbose = TRUE) filterAndTrim(snakemake@input[[1]],snakemake@output[[1]], snakemake@input[[2]], snakemake@output[[2]], maxN=0, rm.phix=TRUE, compress=TRUE, verbose = TRUE)
\ No newline at end of file
...@@ -4,23 +4,26 @@ shell.prefix("source /usr/local/genome/Anaconda2-5.1.0/etc/profile.d/conda.sh;") ...@@ -4,23 +4,26 @@ shell.prefix("source /usr/local/genome/Anaconda2-5.1.0/etc/profile.d/conda.sh;")
configfile: "./config.json" configfile: "./config.json"
SAMPLES=config["SAMPLES"] SAMPLES=config["SAMPLES"]
unit=config["unit"]
rule all: rule all:
input:expand("DATA/{sample}.sra.fastq.gz", sample=SAMPLES), input:#expand("DATA/{sample}.fastq.gz", sample=SAMPLES),
"report/multiqc_report.html", "report/multiqc_report.html",
expand("work/filter/{sample}_1.fastq.gz", sample=SAMPLES),
expand("work/filter/{sample}_2.fastq.gz", sample=SAMPLES),
expand("work/dada/{sample}.rds", sample=SAMPLES), expand("work/dada/{sample}.rds", sample=SAMPLES),
"work/dada/seqtab.fasta", "work/dada/seqtab.fasta",
"work/dada/seqtab.tsv", "work/dada/seqtab.tsv",
"work/dada/metrics.tsv", "work/dada/metrics.tsv",
"work/frogs/affiliation.biom", #"work/frogs/affiliation.biom",
"report/clusters_metrics.html", #"report/clusters_metrics.html",
"report/affiliations_metrics.html", #"report/affiliations_metrics.html",
"report/abundance.biom", #"report/abundance.biom",
"report/abundance.tsv", #"report/abundance.tsv",
"report/tree.nwk" #"report/tree.nwk"
include: "prefetech.smk" #include: "prefetech.smk"
include: "quality.smk" include: "quality.smk"
include: "preprocess.smk" include: "preprocess.smk"
include: "dada2.smk" include: "dada2.smk"
include: "affiliation.smk" #include: "affiliation.smk"
rule cutadapt: rule cutadapt:
input: input:
fwd = "DATA/{sample}.sra_1.fastq.gz", fwd = "DATA/{sample}_1.fastq.gz",
rev = "DATA/{sample}.sra_2.fastq.gz" rev = "DATA/{sample}_2.fastq.gz"
output: output:
cut = "work/cutadapt/{sample}.sra_1.fastq.gz", cut = "work/cutadapt/{sample}_1.fastq.gz",
cutrev = "work/cutadapt/{sample}.sra_2.fastq.gz" cutrev = "work/cutadapt/{sample}_2.fastq.gz"
params: params:
five = lambda wildcards: config["FIVE_PRIMER"][wildcards.sample], five = lambda wildcards: config["FIVE_PRIMER"][wildcards.sample],
three = lambda wildcards: config["THREE_PRIMER"][wildcards.sample] three = lambda wildcards: config["THREE_PRIMER"][wildcards.sample]
...@@ -14,20 +14,23 @@ rule cutadapt: ...@@ -14,20 +14,23 @@ rule cutadapt:
"cutadapt " "cutadapt "
"-g {params.five} " "-g {params.five} "
"-a {params.three} " "-a {params.three} "
"-G {params.five} "
"-A {params.three} "
"--error-rate 0.1 " "--error-rate 0.1 "
"--discard-untrimmed " "--discard-untrimmed "
"--match-read-wildcards " "--match-read-wildcards "
"-o {output} " "-o {output.cut} "
"-p {output.cutrev} "
"{input} " "{input} "
"&& " "&& "
"conda deactivate" "conda deactivate"
rule filter: rule filter:
input: input:
cut = "work/cutadapt/{sample}.sra_1.fastq.gz", cut = "work/cutadapt/{sample}_1.fastq.gz",
cutrev = "work/cutadapt/{sample}.sra_2.fastq.gz" cutrev = "work/cutadapt/{sample}_2.fastq.gz"
output: output:
filt = "work/filter/{sample}.sra_1.fastq.gz", filt = "work/filter/{sample}_1.fastq.gz",
filtrev = "work/filter/{sample}.sra_2.fastq.gz" filtrev = "work/filter/{sample}_2.fastq.gz"
script: script:
"filterAndTrim.R" "filterAndTrim.R"
rule fastqc: rule fastqc:
input: input:
fwd = "DATA/{sample}.sra_1.fastq.gz", fwd = "DATA/{sample}_1.fastq.gz",
rev = "DATA/{sample}.sra_2.fastq.gz" rev = "DATA/{sample}_2.fastq.gz"
output: output:
zip = "work/fastqc/{sample}_fastqc.zip", zip = "work/fastqc/{sample}_{unit}_fastqc.zip",
html = "work/fastqc/{sample}_fastqc.html" html = "work/fastqc/{sample}_{unit}_fastqc.html"
params: params:
output = "work/fastqc/" output = "work/fastqc/"
shell: shell:
...@@ -19,7 +19,7 @@ rule fastqc: ...@@ -19,7 +19,7 @@ rule fastqc:
rule multiqc: rule multiqc:
input: input:
expand("work/fastqc/{sample}_fastqc.zip", sample=SAMPLES) expand("work/fastqc/{sample}_{unit}_fastqc.zip", sample=SAMPLES, unit=unit)
output: output:
html = "report/multiqc_report.html", html = "report/multiqc_report.html",
params: params:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment