00.背景

基于常规转录组（bulk rna）的数据进行染色质可及性（ATAC）的常规分析流程。代码仅供参考。

01.质控

和别的组学一样，拿到原始fastq数据的第一步时进行质控处理，这里可以选择的软件是trim_galore或者trimmomatic。首先需要去相关软件的官网或者是其他教程介绍中下载相应的软件。对于trim_galore，参考代码如下：

nohup  trim_galore -q 25 --phred33 --length $minlen -e 0.1 --stringency 4 --paired -o  ./02.clean  $fq1   $fq2  &
#'$minlen质控时候的最小长度
#'$fq1 左端序列
#'$fq2 右端序列

02.比对

将序列与对应的参考基因组进行比对，为序列注释信息。这里可以选择软件bowtie2和samtools，参考代码如下：

if [ ! -d "03.align" ];then
mkdir 03.align
samtools=软件路径
cut  -f1 sample.list|rush -j 8 \
"bowtie2软件路径 \
-1 ./02.clean/{}_1.clean.fq.gz -2 ./02.clean/{}_2.clean.fq.gz \
-x $bt2_genome|$samtools sort -@ 5 -O bam -o ./03.align/{}.raw.bam "
#比对完马上生成索引
for i in `cut  -f1 sample.list`;do
        $samtools index -@ 5 -b ./03.align/$i.raw.bam
        $samtools flagstat ./03.align/$i.raw.bam >./03.align/$i.raw.stat
done
else
        echo "03.align already exists."
fi
#'sample.list是一个信息列表，主要为了方便传参并行分析
#'相应的参数可以去官网或者相关教程查看

03.去重

ATAC数据的特点是需要去除PCR重复，可以使用的软件是picard。同时，线粒体和低质量的序列也需要过滤。
去重后统计每个样本的序列信息：

##mapping quality>30;注意线粒体
for i in `cut -f1 sample.list`;do
        samtools view -h -f 2 -q 30 ./03.align/$i.rmdup.bam \
        | grep -v -e "chrM" \
        | samtools sort -O bam -@ 10 -o - > ./03.align/$i.last.bam
        samtools index ./03.align/$i.last.bam
        samtools flagstat ./03.align/$i.last.bam > ./03.align/$i.last.stat
done

使用python脚本处理质控信息表格：

import sys
import os
import pandas as pd
sample_list = sys.argv[1]
def get_FQinfo(sample_list):
    with open(f"{sample_list}", 'r') as  f:
        f = f.read().split('\n')
    while '' in f:
        f.remove('')
    sample_id=[]
    final_reads=[]
    duplicate=[]
    read1_mapped=[]
    read2_mapped=[]
    for line in f:
        line = line.split('\t')
        sample_name = line[0]
        sample_id.append(sample_name)
        with open(f"./03.align/{sample_name}.last.stat", 'r') as  f1:
            f1 = f1.read()
        f1 = f1.split('\n')
        final_reads.append(str.split(f1[0]," +")[0])
        read1_mapped.append(str.split(f1[6]," +")[0])
        read2_mapped.append(str.split(f1[7]," +")[0])
        duplicate.append(str.split(f1[3]," +")[0])
    qc_table=pd.DataFrame()
    qc_table["Sample ID"]=sample_id
    qc_table["Read1 Mapped"]=read1_mapped
    qc_table["Read2 Mapped"]=read2_mapped
    qc_table["Final Reads"]=final_reads
    qc_table["Duplicates"]=duplicate
    qc_table.to_csv("./03.align/align_qc_stat.csv")
get_FQinfo(sample_list)

04.peaks分析

主要目的是找到序列中的peaks，用到的软件是macs2。参考代码如下：

cut -f1 sample.list|rush -j 8 \
        "macs2 callpeak -t ./04.peaks/{}.last.bed  -g mm --nomodel --shift  -100 --extsize 200  -n {} --outdir ./04.peaks/"

05.TSS/TES可视化

需要使用deeptools将比对的bam文件转换为bw文件。
然后使用computeMatrix进行分析并使用plotHeatmap和plotProfile进行可视化。
参考代码如下：

cut  -f1 sample.list|rush -j 8 \
"
mkdir -p ./04.peaks/TSS/{}
computeMatrix reference-point --referencePoint TSS -p 15 \
        -b 2000 -a 2000 \
        -R refseq.bed \
        -S ./04.peaks/{}.last.bw \
        -o ./04.peaks/TSS/{}/TSS.gz \
        --skipZeros \
        --outFileSortedRegions ./04.peaks/TSS/{}/Heatmap1sortedRegions.bed
"

for i in `cut -f1 sample.list`;do
plotHeatmap -m ./04.peaks/TSS/$i/TSS.gz -out ./04.peaks/TSS/$i/$i_TSS.Heatmap.pdf --plotFileFormat pdf
plotProfile -m ./04.peaks/TSS/$i/TSS.gz -out ./04.peaks/TSS/$i/$i_TSS.Profile.pdf --plotFileFormat pdf --perGroup
done

06.样本相关性

主要使用deeptools进行样本相关性分析和可视化。

07.peaks结构注释

ChIPseeker注释peaks，仅仅是结构注释：启动子等。这里需要用到R脚本。

07.peaks motif注释

注释motif，使用homer软件。注意物种。

perl ~/miniconda3/envs/atac/share/homer-4.9.1-6/configureHomer.pl  -install mm10

08.peaks差异分析

差异peaks，使用DiffBind软件。

09.footprint分析

使用hintatac软件，注意物种。