ATAC-seq分析
所用到的软件和ChIP-seq一样，所以直接在ChIP环境下进行即可
MD5值校验
md5sum C1_R1.fq.gz
用trim-galore软件质控过滤，然后直接fastqc
vim trim_galore_batch.sh
输入以下内容
#!/bin/bash
# This is for trimming a batch data
for i in C1 C2 BIX1 BIX2
do
trim_galore -q 25 --phred33 --length 25 --stringency 4 --gzip --paired ~/Data/H3K9me2-ATAC/fastq/${i}_R1.fq.gz ~/Data/H3K9me2-ATAC/fastq/${i}_R2.fq.gz -o ~/Data/H3K9me2-ATAC/fastp --fastqc
done
bash trim_galore_batch.sh
image.png
建索引
如果之前已经有索引的话就不用重新建了（之前建的索引不能用了，所以这次重新建的）
nohup >1.log & 后台运行，运行文件追加到1.log
nohup bwa index -a bwtsw ~/genome/mouse/ensembl/index_new/Mus_musculus.GRCm39.dna.toplevel.fa >1.log &
image.png
索引和gtf&gfa文件建立在一个文件夹下面，会生成5个文件，分别是.amb .ann .bwt .pac .sa
比对
比对一般用BWA或者Bowtie2
|也可以写成文件，运行文件，当时一直报错，就一个一个比对的
加上nohup ***** &，让代码在后台运行
nohup bwa mem -v 3 -t 4 ~/genome/mouse/ensembl/index_new/Mus_musculus.GRCm39.dna.toplevel.fa ~/Data/H3K9me2-ATAC/fastp/BIX2_R1_val_1.fq.gz ~/Data/H3K9me2-ATAC/fastp/BIX2_R2_val_2.fq.gz -o ~/Data/H3K9me2-ATAC/bwa/BIX2_bwa.sam >4.log &
查看比对成功率
samtools flagstat -@ 8 *****.sam
转换为bam文件及排序，或者有时候需要合并bam文件（将多个重复合并到一个）
因为bwa生成的是sam文件，所以要转换成bam文件以及进行排序（如果不排序后面用macs2 peak calling的时候可能会报错）
vim bam.sh
for i in C1 C2 BIX1 BIX2
do
samtools view -bS -h ~/Data/H3K9me2-ATAC/bwa/${i}_bwa.sam -o ~/Data/H3K9me2-ATAC/bam/${i}_bwa.bam
done
对bam文件进行排序
vim bam.sort.sh
for i in C1 C2 BIX1 BIX2
do
samtools sort -@ 5 -o ~/Data/H3K9me2-ATAC/bam/${i}_bwa.sort.bam ~/Data/H3K9me2-ATAC/bam/${i}_bwa.bam
done
bash bam.sort.sh
MACS2
vim macs2
输入以下内容
for i in C1 C2 BIX1 BIX2
do
macs2 callpeak -t ${i}_bwa.sort.bam -n ${i}_bwa.sort.bam --shift -100 --extsize 200 --nomodel -B --SPMR -g mm --outdir ~/Data/H3K9me2-ATAC/macs2_2
done
diffbind生成差异peak文件处理
使用bedtools
安装bedtools
onda install -c bioconda bedtools
第一步：文件处理成如下格式（不带表头）将你的染色体位置坐标文件整理成bed格式。
bed格式文件至少包括前3列，分别是：染色体的名字、染色体上的起始位置、染色体上的终止位置。这一步无论用写字板、excel、R等进行处理都可以，文件的后缀名也不重要，因为强行将文件后缀改为bed时，在后面的Linux系统中进行bedtools处理时也会报错。所需的bed格式文件参见下图。
第二步：获得基因组的注释文件（gtf）。
第三步：在Linux系统中处理下载的基因组注释文件，得到蛋白编码基因的位置坐标。
在Linux系统中输入下面的代码，得到蛋白编码基因的位置坐标：
cat ~/genome/mouse/ensembl/index_new/Mus_musculus.GRCm39.107.gtf | grep protein_coding |perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' >Mus_musculus.position
第四步：
mv diff_bed.txt diff.bed
perl -p -i -e 's/ /\t/g' diff.bed
此时如果直接进行下一步的话可能会报错，因为用windows生成的tab分隔文件和linux生成的tab分隔文件不同，会有隐形符号，可用如下代码查看
cat -e diff.bed
用如下代码去掉特殊符号
cat diff.bed | tr -d '\r' > clean_diff.bed
第五步：在Linux系统中利用bedtools得到包含染色体位置坐标的蛋白编码基因。
首先需要启动自己安装了bedtools软件的conda小环境，然后输入下面的代码：
bedtools intersect -a clean_diff.bed -b Mus_musculus.position -wa -wb > symbol
