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软件，注意物种。
