ATAC-seq经典分析流程（上）
之前挖过一篇NC的文章，我们利用ATAC-seq和RNA-seq联合分析的手段筛选了与所关注性状的密切相关的基因。后来，通过还做了一系列的分子实验，结果发现这两种手段联合分析筛选出来的基因具有很强的功能。所以，这里我想和大家分享一下我们之前的分析过程，以供初学者在分析中有所参考！
下面我们步入正题：
1）下载sra数据
2）sra格式转fastq格式
mkdir 2_FASTQ
fastq-dump --split-3 --gzip -O ./2_FASTQ ./1_sra/*.1--split-3 把双端测序数据拆分成两个文件,对单端测序数据不起作用.fastq-dump默认会把双端测序结果保存到一个文件里
--gzip 输出gz格式, 能够节省空间的同时也不会给后续比对软件造成压力
-o|--outdir : 输出到指定文件夹
3) FastQC检测数据质量
mkdir 3_raw_fastqc
fastqc -o ./3_raw_fastqc -t 4 ./2_FASTQ/*.fastq.gz -o --outdir FastQC生成的报告文件的储存路径，生成的报告的文件名是根据输入来定的
-t --threads 选择程序运行的线程数，
-q --quiet 安静运行模式，一般不选这个选项的时候，程序会实时报告运行的状况
4) Trimmomatic过滤低质量片段和接头
当前，由于ATAC-seq普遍使用Illumina的Nextera文库，Nextera测序接头比例过高，进行准确的read比对前需将接头删除。目前，去除接头的工具大多采用不同的动态编程，包括 cutadapt，AdapterRemoval v2 ，Skewer和trimmomatic，它们都需要输入已知的接头序列。Trim-galore，fastp可以自动识别序头序列。对于Nextera和Truseq文库使用trimmomatic和内置接头序列是一种直接简单的办法，不仅可以去除接头，同时可以去除低质量的碱基。各种read过滤工具在有效去除低质量和污染接头序列的性能方面通常表现差不多。
其实我的数据本身已经进行过质控了(QC的为一个范围且短于提供给NCBI上的数据)：
为了整个步骤的完整性，这里再进行一次数据过滤，以Trimmomatic为例：
# 使用简单命令安装软件
Sudoapt install trimmomatic
# 查看软件安装目录
dpkg -Ltrimmomatic
# 单端过滤代码参考
java -jar /usr/share/java/trimmomatic-0.36.jar SE -threads 4 -phred33 /home/wang/桌面/SRX091.fastq.gz /home/wang/桌面/resultout.fq.gz ILLUMINACLIP:TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:25
# 双端过滤代码参考
java -jar/usr/share/java/trimmomatic-0.36.jar PE -threads 4 -phred33 s_1_1_sequence.txt.gz s_1_2_sequence.txt.gz lane1_forward_paired.fq.gz lane1_forward_unpaired.fq.gz lane1_reverse_paired.fq.gz
lane1_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3
TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:25 结果发现，使用以上参数有三个文件中的数据100%被过滤掉了，但是其他数据draped率不超2%，个人推测trimmomatic参数设置过于严格不适合我的短片段数据。
因此，换用了目前比较新的fastp过滤软件，如需了解详细信息找找度娘咯！
#软件下载
wget https://opengene.org/fastp/fastp
#增加该文件的可执行权限
chmod a+x ./fastp
ls
./fastp -h
#将fastp软件放在Fastp文件夹里，以便后续配置路径
mkdir Fastp
mv fastp Fastp
ls
pwd #查看软件的安装路径
cd
vim .bashrc #将fastp的路径配置在.bashrc
source .bashrc
#以后在使用fastp的时候就不需要在输入他的完整路径了
fastp#会同时生成包含.json和.html的质控报告，便于结果的可视化查看
fastp -i ./GSE15/fastq/SRR126${i}_1.fastq.gz -o ./GSE158/Fastp_trim_fq/SRR126${i}_1.fq.gz -I ./GSE158/fastq/SRR126${i}_2.fastq.gz -O ./GSE158/Fastp_trim_fq/SRR126${i}_2.fq.gz -h fastp_out_./GSE158/Fastp_trim_fq/SRR126${i}.html -j fastp_out_./GSE154/Fastp_trim_fq/SRR126${i}.json#--------------------------重要参数解析--------------------------------#
#默认会进行低质过滤，如果不希望进行，可使用参数-Q或者“disable_quality_filtering”目前支持过滤通过N碱基的数量和不合格碱基的百分比。
#-q,--disable_quality_filtering 设置碱基的质量值。默认平均质量值>=Q15是合格的
#-u,--unqualified_percent_limit 允许多少百分比的碱基不合格(0~100)。默认40%
#或者通过平均质量值过滤数据
-e, --average_qual 如果一条reads的平均质量值小于平均质量值，即将这条reads或paire reads丢弃。默认0
# 长度过滤
默认进行长度过滤，可以通过参数-L或--disable_length_filtering禁用。最小长度要求可以用-l或--length_required指定
# 低复杂过滤
默认不进行低复杂过滤，可以通过参数-y或--low_complexity_filter进行低复杂序列过滤（ase[i] != base[i+1]）。
低复杂序列阈值可以使用参数-Y或者--complexity_threshold进行指定，阈值的范围需在0~100之间，默认是30
#接头过滤
默认进行接头过滤，但可以使用参数-A或者--disable_adapter_trimming禁用，接头序列可自动地在PE/SE数据中发现。
5) BWA-MEM进行序列的参考基因组比对
BWA-MEM和Bowtie2 对于短的双端read存储效率高且快速。两个比对工具的软限位策略允许在read的两端有突出碱基，这可以进一步提高unique mapping rate。我们建议，unique mapping rate达到80％以上时认为ATAC-seq实验成功。对于哺乳动物物种，基于经验和计算估计，建议染色质开放区域检测和差异分析至少需要5000万mapped read，TF footprinting至少需要2亿。
之前我在分析Chip-seq数据时一直用的Bowtie2，后面换了物种Bowtie2的index在建立时一直报错，可能原因是我用自己的笔记本跑数据，内存不够用，加之我的reads长度只有三四十，所以转用BWA来进行mapping分析！
#一些推文表示Bwa-MEM和bowtie2的比对结果差不多
#建索引
bwa index ./botie2_index/Bos_taurus.ARS-UCD1.2.dna.toplevel.fa -p./BWA_index/Bos_taurus_BWA #可以不加-p genome，这样建立索引都是以ref.fa为前缀
#双端比对
bwa mem -t 8 -R "@RG\tID:${filename}\tSM:${filename}\tLB:WXS\tPL:Illumina" ./Bos_taurus.ARS-UCD1.2.dna.toplevel.fa ./GSE158/Fastp_trim_fq/SRR12${i}_1.fq.gz \
./GSE1584/Fastp_trim_fq/SRR126${i}_2.fq.gz | samtools sort -O bam -@ 10 -o ./GSE158/BWA_mapping/SRR126${i}.bam
#给bam文件建索引
samtools index -b ./BWA_mapping/SRR126${i}.bam
#输出比对结果
samtools flagstat ./BWA_mapping/SRR126${i}.bam > ./BWA_mapping/SRR126${i}.bam.stat6）比对后数据的处理&质检
序列比对后，可利用Picard和SAMtools 获取BAM文件的基本指标，包括唯一比对率，duplicated read的百分比以及片段大小分布等。
在完成基因比对后，要对ATAC-seq的数据完成两个基本的统计（ATACseqQC），1是测序片段长度分布；2是数据在转录起始位点的信号强度：
i. 成功的ATAC-seq实验应生成片段大小分布图，其具有递减的和周期性的峰，对应于无核小体区域（NFR）（<100 bp）和单核、双核和三核小体（〜200， 400，600碱基对）。
ii. NFR的片段应该在基因的转录起始位点（TSS）周围富集，而核小体结合区域的片段应该在TSS处被形成低谷，TSS周围的侧翼区域会稍微富集。
iii. 最后，对于正链和负链，read应分别偏移 +4 bp和 -5 bp，以便实现TF足迹和基序的碱基对解析相关分析，因为Tn5转座酶对缺口进行DNA修复产生了9 bp重复。大多数上述质量控制和分析报告可以使用MultiQC 汇总以进行集成的、用户友好的交互式的呈现。
7) 利用Genrich进行peak-calling
Peak-calling需要经历多个步骤，包括去除线粒体基因组和PCR重复等。Genich是新推出的软件，优点在于只需一行命令便可完成去除线粒体reads，PCR重复，多重mapping reads以及多生物学重复的处理，输出结果为ENCODE narrowPeak format。另外发现，这一软件calling速度较快.
#使用Genrish前需要Sort比对结果 by read name
samtools sort -n SRR126{i}.bam > SRR126{i}.sort.bam
Genrich -t ./BWA_bam/SRR126{i}.sort.bam -o ./SRR126{i}.narrowPeak -f ./SRR12697{i}.log -r -e MT -q 0.05
wc SRR126{i}.narrowPeak #结果是4000多个peaksGenrich -t ./BWA_bam/SRR126{i}.sort.bam -o ./SRR126{i}.narrowPeak -f ./SRR126{i}.log -j -r -y -e MT -p 0.01
wc SRR126{i}.narrowPeak #结果1w以上-t: 输入文件
-o: 输出文件
-j: 使用ATAC-seq模式
-r: remove PCR duplicates
-e MT/chrM : to exclude mitochondrial chromosome
-p 0.01: 筛选p值小于0.01.(默认)，一旦指定p值，就不再考虑q值
-q 默认是1
-y: Keep unpaired alignments (def. false)
至此，我们就拿到了所有样本的peaks，后续可以进行下游的注释分析和差异分析了。今天的分享先到此为止，希望还有机会和小伙伴们一起学习下游的分析！
腾讯云开发者
