ATAC学习

2022-06-25 本文已影响0人小qqq

最近一周开始学习ATAC数据分析，目前只是尝试了简单分析，新手小白，如果有小失误欢迎留言~

1.准备环境

创建conda环境（如果已有成熟环境不需再次创建）

conda create -n atac -y python=2 bwa
conda info --envs

激活环境

conda activate atac

2.安装所需要的包

conda install -c bioconda -y fastqc
conda install -c bioconda -y trimmomatic
conda install -y hisat2
conda install -c bioconda -y samtools
conda install -c bioconda -y picard
conda install -c bioconda -y macs2
conda install -c bioconda -y bedtools
conda install -c bioconda -y deeptools
conda install -c bioconda -y homer
conda install -y sambamba
conda install -y deeptools

3.对下机数据rowdata进行简单质控--可选

执行命令

fastqc -t 8 -o /public/home/sss/ss/ATAC/raw ATAC1-2_FKDL210215276-1a_1.fq.gz
fastqc -t 8 -o /public/home/sss/ss/ATAC/rawATAC1-2_FKDL210215276-1a_2.fq.gz

说明

-t: 线程
-o: 存放路径，不用指定前缀，默认为.fastq.gz前面的字段
*.fastq.gz:fastqc可以同时接受多个fastq.gz文件，因此采用正则表达式*表示全部

结果会产生两个文件：

image

将HTLM文件下载到本地打开

image

注： 接头含量太多需要过滤以及去接头

4.使用trim_galore去除接头，trim_galore可以自动识别接头

下载trim_galore

curl -fsSL[https://github.com/FelixKrueger/TrimGalore/archive/0.6.6.tar.gz](https://github.com/FelixKrueger/TrimGalore/archive/0.6.6.tar.gz)-o trim_galore.tar.gz
tar xvzf trim_galore.tar.gz

######################################

去除NexteraPE

trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired./ATAC1-2_FKDL210215276-1a_1.fq.gz ./ATAC1-2_FKDL210215276-1a_2.fq.gz -o ./

产生新的.gz文件

image

5.对去除接头后的数据质检

执行命令

 fastqc -t 8 -o /public/home/sss/ss/ATAC/raw ATAC1-2_FKDL210215276-1a_1_val_1.fq.gz \fastqc -t 8 -o/public/home/sss/ss/ATAC/raw ATAC1-2_FKDL210215276-1a_2_val_2.fq.gz

产生两个文件：

image

将htlm文件下载并打开:

image

目前已没有接头序列，重复序列比较多，后面比对后需要进一步去重

6.使用hisat2比对

最重要建立参考基因组的索引

mkdir reference && cd reference

参考基因组

wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.con

注释文件

wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.gff3

Hisat2 建立索引

hisat2-build -p 10 reference/all.con reference/Osativa

建立后的索引文件

image

建立align.sh文件

clean data

/public/home/sss/software/anaconda3/envs/chipseq/bin/hisat2 --dta -t-x /public/home/sss/ss/genome/NIP -1/public/home/sss/ss/ATAC/cleandata/ATAC1-2_FKDL210215276-1a_2_val_2.fq.gz -2 /public/home/sss/sss/ATAC/cleandata/NIP-1_FRAS210228372-2r_2.clean.fq.gz | /public/home/sss/anaconda2/envs/chipseq/bin/samtools sort -@ 8 -O bam -o-/public/home/sss/ss/ATAC/cleandata/ATAC_1-2.bam

注意：由于sam文件较大，因此我们直接跳过sam，使用samtools转换为排序后的bam文件，注意【|】不要漏掉，且排序是依据reads比对到基因组的位置排的

-O: 表示输出的bam文件

-o: 输出的bam文件名

-t与@: 线程数

比对后产生：

image

7.比对后的序列去重和排序

选用sambamba来去重复

sambamba markdup -r ATAC_1-2.bam ATAC_1-2.sambamba.rmdup.bam

smarttool提取可靠的比对结果

samtools view -f 2 -q 30  -o ATAC_1-2.rmdup.q30.bamATAC_1-2.sambamba.rmdup.bam

smarttool对比对结果排序

samtools sort -O bam -@ 4 -o ./ATAC_1-2.fraw.bam ATAC_1-2.rmdup.q30.bam

产生的bam 文件:

image

8.计算插入片段长度，FRiP值，IDR计算重复情况

提取片段的长度

从sam文件中求得insert size：在“read mapped in proper pair”的前提下，取第九列

samtools  view ATAC_1-2.fraw.bam |awk -F'\t' 'function abs(x){return ((x < 0.0) ? -x : x)} {print$1"\t"abs($9)}' | sort | uniq | cut -f2 > fragment.length.txt

提取出片段长度，R代码如下：

image

得到结果图:

image

FRiP_score的计算

Reads=$(bedtools intersect -a ATAC_1-2_summits.bed -bATAC_1-2_peaks.narrowPeak |wc -l|awk '{print $1}')

totalReads=$(wc -l ATAC_1-2_summits.bed|awk '{print $1}')

echo ${Reads} ${totalReads} ATAC_1-2 'FRiP_score:' $(echo"scale=2;100*${Reads}/${totalReads}"|bc)'%'

9.使用过滤好的数据call peak

samtools构建index

samtools index   ATAC_1-2.fraw.bam

bedtools转换为bed 文件

bedtools bamtobed -i ATAC_1-2.fraw.bam > ATAC_1-2.bed

使用macs2 callpeak

macs2 callpeak -t ATAC_1-2.bed  -g380699722 --nomodel --shift -100 --extsize 200 -n  ATAC_1-2 --outdir  ./peak

得到的文件：

image

10.deeptools 进行TSS可视化

将 bam 文件转化为 bw 文件

bamCoverage -p 8 --bam ATAC_1-2.sambamba.rmdup.bam --binSize 10--centerReads --smoothLength 14 --normalizeUsing RPKM -ocoverageATAC_1-2.bigwig

bamCoverage -p 8 --bam ATAC_1-2.fraw.bam --binSize 10 --centerReads--smoothLength 14 --normalizeUsing RPKM -o coverageATAC_1-2.fraw.bigwig

建立.deeptools_TSS.sh文件

计算reference-point---需要基因组的Refseq文件--下载或者自己创制

mkdir -p  /public/home/sss/ss/ATAC/tss

cd   /public/home/sss/ss/ATAC/tss

 source activate atac

computeMatrix reference-point --referencePoint TSS  -p 15  \

-b 2000 -a 2000    \

-S /public/home/sss/ss/ATAC/cleandata/coverageATAC_1-2.bw  \

-R  /public/home/sss/ss/ATAC/NIP.deptools.bed \

--skipZeros  \

-o matrix1_test_TSS.gz  \

--outFileSortedRegions regions1_test_genes.bed

##     both plotHeatmap andplotProfile will use the output from  computeMatrix

plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.png

plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.pdf --plotFileFormat pdf  --dpi 720

plotProfile -m matrix1_test_TSS.gz -out test_Profile.png

plotProfile -m matrix1_test_TSS.gz -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720

##########################################################################################

结果文件：

image

11.peaks注释

结构注释----会将peak所落在基因组上的区域结构注释出来，比如说启动子区域，UTR区域，内含子区域等等。同时，也会将与peak最临近的基因注释出来，非常好用

使用gff/gtf构建一个TxDb

R代码：

if (!requireNamespace("BiocManager", quietly = TRUE))

 install.packages("BiocManager")

BiocManager::install()

BiocManager::install(c("GenomicFeatures","AnnotationDbi"))

#install.packages("AnnotationDbi")

library(clusterProfiler)

library(ChIPseeker)

library(GenomicFeatures)

library(AnnotationDbi)

txdb <- makeTxDbFromGFF("D:/生信学习/ATAC/ricegenome/NIP.gtf",

                       format="gtf")    #也可以使用gtf

keytypes(txdb)    #感兴趣的话，可以用以下方法探索txdb都包含了什么内容

keys(txdb)

读入单个summits文件

peaks <- readPeakFile("D:/生信学习/ATAC/ATAC_1-2_summits.bed")

结构注释

peakAnno <- annotatePeak(peaks,

                         TxDb=txdb,

                        tssRegion=c(-1000, 1000))

#注释完，进行可视化，多种图可供选择

plotAnnoBar(peakAnno)

plotDistToTSS(peakAnno)

vennpie(peakAnno)

plotAnnoPie(peakAnno)

install.packages("ggupset")

library(ggupset)

upsetplot(peakAnno)

install.packages("ggimage")

library(ggimage)

upsetplot(peakAnno, vennpie=TRUE)

library(ggplot2)

最后将我们的注释结果转为数据框，便于查看

df <- as.data.frame(peakAnno)

将注释到的基因提取出来（第14列），用于后续功能分析

gene <- df[,14]

head(gene)

获得结果文件;

image

功能注释-----需要物种的OrgDb注释库

R代码：

#GO 富集分析

library(AnnotationHub)

hub <- AnnotationHub()

query(hub, "Oryza sativa")

nip_org <- hub[['AH96213']] ##下载目标物种到org数据

keytypes(nip_org)#查看可使用的key

length(keys(nip_org))

##我们可以使用riceidconverter来转换ID，甚至于直接检索GO注释：

require(clusterProfiler)

head(gene)

gene_SYMBOL <- riceidconverter::RiceIDConvert(gene,fromType='MSU',toType = 'SYMBOL')

head(gene_SYMBOL)

gene_SYMBOL1 <- as.character(gene_SYMBOL$SYMBOL)

head(gene_SYMBOL1)

gene_SYMBOL2 <- gene_SYMBOL1[grep("LOC*" ,gene_SYMBOL1)]

C4_3DAP_go = enrichGO(gene_SYMBOL2,keyType = "SYMBOL",OrgDb=nip_org, pvalueCutoff=1, qvalueCutoff=1)

dotplot(C4_3DAP_go,showCategory = 20)

cnetplot(C4_3DAP_go,circular=T,  ###画为圈图

         colorEdge=T,showCategory =15)

cnetplot(C4_3DAP_go, showCategory = 5)

write.table(C4_3DAP_go, 'C4_3DAP_go.txt', sep = '\t', quote = FALSE,row.names = FALSE)

kegg

install.packages("R.utils")

R.utils::setOption("clusterProfiler.download.method",'auto')

head(gene)

gene_list <- bitr(geneID = gene_SYMBOL1,

                  fromType ="SYMBOL",

                  toType =c("ENTREZID"),

                  OrgDb = nip_org)

head(gene_list)

genelist2<-(gene_list$ENTREZID)

head(genelist2)

compKEGG <- enrichKEGG(genelist2,organism="osa",keyType="kegg",pvalueCutoff=0.05,pAdjustMethod="BH",qvalueCutoff=0.2)

dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway EnrichmentAnalysis")

################################################################

获得结果文件：

image

12.homer寻找motif

使用homer 进行 motif分析

conda 环境下安装 homor环境

conda install -c bioconda -y homer

awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}'\

ATAC_1-2_peaks.narrowPeak > ATAC_1-2_peaks.tmp

findMotifsGenome.pl ATAC_1-2_peaks.tmp \

/public/home/sss/ss/genome/NIP.fa \

./ATAC_1-2.motif \

-size 200 -len 8,10,12

#########################################

image

产生的文件：

image

     将该目录下所有文件下载到本地，打开homerResults.html 文件，即可得到如下结果展示：

image

输出结果按照p-value排序，最后一列是一个链接到motif文件的超链接，可以从这个文件中找到包含此motif的其他序列。

在Best Match/Details 列中，HOMER将会展示与denovo motif最匹配的已知motif（该列还包含有一个名为More information 的超链接，

点击后会出现以下页面：该页面包含motif的一些基本信息，如链接到motfi文件的超链接，且可以生成motif logo的pdf格式。

参考文章链接：
原文链接：https://blog.csdn.net/u012110870/article/details/102804623
原文链接： https://www.jianshu.com/p/9aa719faa4b5
原文链接：https://cloud.tencent.com/developer/article/1360799
原文链接：[https://www.jianshu.com/p/9a31f5f01e7b

ATAC学习

1.准备环境

创建conda环境（如果已有成熟环境不需再次创建）

激活环境

2.安装所需要的包

3.对下机数据rowdata进行简单质控--可选

4.使用trim_galore去除接头，trim_galore可以自动识别接头

下载trim_galore

去除NexteraPE

产生新的.gz文件

5.对去除接头后的数据质检

6.使用hisat2比对

最重要建立参考基因组的索引

参考基因组

注释文件

Hisat2 建立索引

建立后的索引文件

建立align.sh文件

clean data

7.比对后的序列去重和排序

选用sambamba来去重复

smarttool提取可靠的比对结果

smarttool对比对结果排序

产生的bam 文件:

8.计算插入片段长度，FRiP值，IDR计算重复情况

提取片段的长度

FRiP_score的计算

9.使用过滤好的数据call peak

samtools构建index

bedtools转换为bed 文件

使用macs2 callpeak

得到的文件：

10.deeptools 进行TSS可视化

将 bam 文件转化为 bw 文件

建立.deeptools_TSS.sh文件

11.peaks注释

使用gff/gtf构建一个TxDb

读入单个summits文件

结构注释

最后将我们的注释结果转为数据框，便于查看

将注释到的基因提取出来（第14列），用于后续功能分析

功能注释-----需要物种的OrgDb注释库

kegg

12.homer寻找motif

使用homer 进行 motif分析

conda 环境下安装 homor环境

猜你喜欢

热点阅读