ATAC-seq Snakemake 流程

参考1
参考2

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