ATAC-seq数据分析流程

1.数据质量评估

1.fastqc质控

# 进入到到自己的文件夹，使用FastQC软件对fastq文件进行质量评估，结果输出到qc/文件夹下
vim qc.sh
i
fastqc -t 6 -o ./ *.fq.gz
# 按esc
:wq
nohup bash qc.sh >qc.log &
# 报告整合
multiqc *zip

cat qc.sh
cat qc.log
# 质控的结果给每个fq文件生成一个zip文件和一个html网页报告
# zip文件里是一些报告图片等，主要需要看的是html网页报告，可以将其下载到本地通过浏览器打开

2.使用Trim Galore修剪低质量reads

这里选择trim_galore软件，自动批量运行：

# 到文件夹下先生成一个变量,为样本ID
cd /Data/wu_lab/zyy/project/ATAC/
ls /Data/wu_lab/zyy/project/ATAC/*_1.fq.gz |awk -F '/' '{print $NF}' | cut -d '_' -f 1,2 >ID
cat ID

vim trim_galore.sh

rawdata=/Data/wu_lab/zyy/project/ATAC/
cleandata=/Data/wu_lab/zyy/project/ATAC/clean/
ID=/Data/wu_lab/zyy/project/ATAC/ID
cat $ID | while read id
do
  trim_galore -q 20 --length 33 --max_n 3 --stringency 3 --fastqc --gzip --paired -o ${cleandata}    ${rawdata}/${id}_1.fq.gz     ${rawdata}/${id}_2.fq.gz
done

nohup bash trim_galore.sh >trim_galore.log &

查看过滤后的qc

3.数据比对到参考基因组

使用比对软件（如BWA、Bowtie2等）将读段比对到参考基因组上。这个过程中，读段会根据其序列相似性被映射到基因组的特定位置。

#使用Bowtie2进行序列比对
#构建参考基因组的索引或者已经有现成的索引。
#如果没有，需要先使用 Bowtie2 构建索引。
#1.构建索引
# 进入参考基因组目录
mkdir bowtie2Index
vim bowtie2Index.sh

fa=Homo_sapiens.GRCh38.dna.primary_assembly.fa
fa_baseName=GRCh38.dna
bowtie2-build  -f ${
   fa} bowtie2Index/${
   fa_baseName}

#运行
nohup bash bowtie2Index.sh > bowtie2Index.log &

## ----2.比对
#2.1 bowtie2批量处理多个样本：
cd /Data/wu_lab/zyy/project/ATAC/align/

vim bowtie2.sh

# 设置输入数据目录
inputdir="/Data/wu_lab/zyy/project/ATAC/clean/"
# 设置ID文件路径
ID=/Data/wu_lab/zyy/project/ATAC/ID
# 设置输出结果目录
outdir="/Data/wu_lab/zyy/project/ATAC/align/"
# 设置bowtie2索引目录
index=/Data/wu_lab/zyy/database/GRCh38.p14/bowtie2Index/GRCh38.dna

# 确保输出目录存在
mkdir -p "$outdir"

# 读取ID文件并进行批量比对
cat $ID | while read id; do

    # 执行bowtie2比对并将结果排序为BAM格式
    bowtie2 -p 5 --very-sensitive -X 2000 -x ${
   index} \
    -1 ${
   inputdir}/${
   id}_1_val_1.fq.gz \
    -2 ${
   inputdir}/${
   id}_2_val_2.fq.gz \
    --phred33 --end-to-end --no-mixed --no-discordant --no-unal | samtools sort -@ 16 -O bam -o ${
   outdir}/${
   id}.bam
       
# 2.2 为排序后的BAM文件创建索引
samtools index ${
   outdir}/${
   id}.bam
samtools flagstat ${
   outdir}/${
   id}.bam > ${
   outdir}/${
   id}.raw.stat

 echo "Alignment and sorting for $id completed."
done

nohup bash bowtie2.sh > bowtie2.sh.log & 

4.序列筛选

vim redup.sh
# 设置输入数据目录
inputdir="/Data/wu_lab/zyy/project/ATAC/align/"
# 设置ID文件路径
ID=/Data/wu_lab/zyy/project/ATAC/ID
# 设置输出结果目录
outdir="/Data/wu_lab/zyy/project/ATAC/align/"

# 读取ID文件并进行批量比对
cat $ID | while read id; do

# 1. samtools idxstats 命令用于从已建立索引的 BAM 文件中检索并打印统计信息。具体来说，它会输出一个制表符分隔的表格，每行包含以下四列信息：
#参考序列名称：即染色体或参考序列的名称。
#序列长度：该参考序列的长度。
#比对上的reads数：在该参考序列上比对上的reads数量。
#未比对上的reads数：在该参考序列上未比对上的reads数量。

#samtools idxstats *.rmdup.bam
#查看发现有线粒体reads

# 2. 过滤线粒体reads
 ##先看过滤前后是否有差别
#samtools view -h -f 2 -q 30  WT.bam |grep -v MT|wc
#samtools view -h -f 2 -q 30  WT.bam |wc

samtools view -@ 6 -h ${
   inputdir}/${
   id}.bam | grep -v MT | samtools sort -@ 6 -O bam -o ${
   outdir}/${
   id}.rmChrM.bam
samtools index ${
   outdir}/${
   id}.rmChrM.bam
samtools flagstat ${
   outd