Understand the method
snakemake file
shell.executable("/bin/bash")
configfile: "./config.yaml"
(SAMPLES,READS,) = glob_wildcards("01fastq/{sample}_{read}.fastq.gz")
READS=["1","2"]
rule all:
input:
expand("06peak_macs2/{sample}_macs2_peaks.broadPeak", sample=SAMPLES)
rule fastqc:
input: "01fastq/{sample}_1.fastq.gz", "01fastq/{sample}_2.fastq.gz"
output: "02fqc/{sample}_1_fastqc.zip", "02fqc/{sample}_2_fastqc.zip"
log: "00log/{sample}_fastqc"
threads: 2
params : jobname = "{sample}"
message: "fastqc {input}: {threads}"
shell:
"""
# fastqc works fine on .gz file as well
module load fastqc
fastqc -o 02fqc {input[0]} 2> {log}
fastqc -o 02fqc {input[1]} 2> {log}
"""
rule fastp:
input:
fwd="01fastq/{sample}_1.fastq.gz",
rev="01fastq/{sample}_2.fastq.gz"
output:
fwd="03trim/{sample}_1_good.fq.gz",
rev="03trim/{sample}_2_good.fq.gz",
html="03trim/{sample}.html",
json="03trim/{sample}.json"
log:
"00log/{sample}_fastp"
message: "trim_adaptor {input}"
shell:
"""
ml fastp
fastp -w 10 -h {output[html]} -j {output[json]} -i {input[fwd]} -I {input[rev]} -o {output[fwd]} -O {output[rev]} 2> {log}
"""
## the later step will remove chrM from the bam file and coordinate sort the bam
## so I did not cooridnate sort the bam at this step to save some time.
rule align:
input: "03trim/{sample}_1_good.fq.gz", "03trim/{sample}_2_good.fq.gz"
output: "04aln/{sample}.sorted.bam", "00log/{sample}.align"
threads: 10
params: jobname = "{sample}"
message: "aligning {input}: {threads} threads"
log:
bowtie2 = "00log/{sample}.align",
markdup = "00log/{sample}.markdup"
shell:
"""
## samblaster mark duplicates for read id grouped reads. I do not coordinate sort the bam
ml bowtie/2 samblaster samtools
bowtie2 --threads 10 --very-sensitive -x {config[idx_bt2]} -1 {input[0]} -2 {input[1]} 2> {log.bowtie2} \
| samblaster 2> {log.markdup} \
| samtools view -@ 10 -Sb - > {output[0]}
"""
# check number of reads mapped by samtools flagstat
rule flagstat_bam:
input: "04aln/{sample}.sorted.bam"
output: "04aln/{sample}.sorted.bam.flagstat"
log: "00log/{sample}.flagstat_bam"
threads: 2
params: jobname = "{sample}"
message: "flagstat_bam {input}: {threads} threads"
shell:
"""
ml samtools
samtools flagstat {input} > {output} 2> {log}
"""
## shifting the reads are only critical for TF footprint, for peak calling and making bigwigs, it should be fine using the bams without shifting
# https://sites.google.com/site/atacseqpublic/atac-seq-analysis-methods/offsetmethods
rule remove_chrM_bam:
input: "04aln/{sample}.sorted.bam"
output: "05aln_exclude_chrM/{sample}_exclude_chrM.sorted.bam", "05aln_exclude_chrM/{sample}_exclude_chrM.sorted.bam.bai"
log: "00log/{sample}_exclude_chrM_bam.log"
threads: 10
message: "excluding chrM from bam {input} : {threads} threads"
params: jobname = "{sample}"
shell:
"""
ml samtools samblaster
# remove duplicates and reads on chrM, coordinate sort the bam
# samblaster expects name sorted bamq
samtools view -h {input} | samblaster --removeDups \
| grep -v -P '\tchrM\t' \
| samtools view -Sb -F 4 - \
| samtools sort -m 15G -@ 10 -T {input}.tmp -o {output[0]}
samtools index {output[0]}
"""
rule flagstat_chrM_exclude_bam:
input: "05aln_exclude_chrM/{sample}_exclude_chrM.sorted.bam"
output: "05aln_exclude_chrM/{sample}_exclude_chrM.sorted.bam.flagstat"
log: "00log/{sample}_exclude_chrM_flagstat_bam"
threads: 5
params: jobname = "{sample}"
message: "flagstat_bam {input}: {threads} threads"
shell:
"""
ml samtools
samtools flagstat {input} -@ 5 > {output} 2> {log}
"""
# https://github.com/taoliu/MACS/issues/145
rule call_peaks_macs2:
input: "05aln_exclude_chrM/{sample}_exclude_chrM.sorted.bam", "05aln_exclude_chrM/{sample}_exclude_chrM.sorted.bam.bai"
output: bed = "06peak_macs2/{sample}_macs2_peaks.broadPeak"
log: "00log/{sample}_call_peaks_macs2.log"
params:
name = "{sample}_macs2",
jobname = "{sample}"
message: "call_peaks macs2 {input}: {threads} threads"
shell:
"""
ml macs
## for macs2, when nomodel is set, --extsize is default to 200bp, this is the same as 2 * shift-size in macs14.
macs2 callpeak -t {input[0]} \
--keep-dup all -f BAMPE -g {config[macs2_g]} -B \
--outdir 06peak_macs2 -n {params.name} -p {config[macs2_pvalue]} \
--broad --broad-cutoff {config[macs2_pvalue_broad]} &> {log}
"""
config file
# bowtie1 index
idx_bt2: ~/index/bowtie2/humanized/ucsc
# genome size for bedtools slop
genome_size: ~/index/ori_genomic/humanized_mice/ucsc/hg38.chrom.sizes.human.mm10.chrom.sizes.mouse
## genome fasta for nucleoATAC
genome_fasta: ~/index/ori_genomic/humanized_mice/ucsc/hg38_human_mm10_mouse.fa
macs2_g: 4.57e+09#for human just "hs", check the help
macs2_pvalue: 1e-5
macs2_pvalue_broad: 1e-5
递交文件
#!/bin/bash
#SBATCH --job-name=sbatch_ATAC_seq_flow
#SBATCH --time=24:00:00
#SBATCH --cpus-per-task=40
#SBATCH --mem=60g
module load snakemake singularity
ml python/3.6
snakemake -j 40 --use-singularity
