ATAC-Seq分析教程：用ChIPseeker对peaks进行注释和可视化

来源：生信拾光评论5,919

ATAC-Seq分析教程系列

ATAC-Seq分析教程：ATAC-seq的背景介绍以及与ChIP-Seq的异同

ATAC-Seq分析教程：原始数据的质控、比对和过滤

ATAC-Seq分析教程：用MACS2软件call peaks

ATAC-Seq分析教程：对ATAC-Seq/ChIP-seq的质量评估（一）phantompeakqualtools

ATAC-Seq分析教程：对ATAC-Seq/ChIP-seq的质量评估（二）ChIPQC

ATAC-Seq分析教程：重复样本的处理-IDR

ATAC-Seq分析教程：用ChIPseeker对peaks进行注释和可视化

ATAC-Seq分析教程：用网页版工具做功能分析和motif分析

ATAC-Seq分析教程：差异peaks分析——DiffBind

ATAC-Seq分析教程：ATAC-Seq、ChIP-Seq、RNA-Seq整合分析

 

上一步骤用IDR对重复样本peaks的一致性进行了评估，同时得到了merge后的一致性的peaks——sample-idr，接下来就是对peaks的注释。这篇主要介绍用Y叔的R包ChIPseeker对peaks的位置（如peaks位置落在启动子、UTR、内含子等）以及peaks临近基因的注释。

ChIPseeker

ChIPseeker虽然最初是为了ChIP-seq注释而写的一个R包，但它不只局限于ChIP-seq，也可用于ATAC-Seq等其他富集peaks注释，也可用于lincRNA注释、DNA breakpoints的断点位置注释等所有genomic coordination的注释，另外提供了丰富的可视化功能。

使用方法

使用ChIPseeker需要准备两个文件：一个就是要注释的peaks的文件，需满足BED格式。另一个就是注释参考文件，即需要一个包含注释信息的TxDb对象。Bioconductor提供了30个TxDb包，如果其中有研究的物种就可以直接下载安装此物种的TxDb信息。如果研究的物种没有已知的TxDb，可以使用GenomicFeatures包的函数（makeTxDbFromUCSC，makeTxDbFromBiomart）制作TxDb对象：

makeTxDbFromUCSC： 通过UCSC在线制作TxDb
makeTxDbFromBiomart: 通过ensembl在线制作TxDb
makeTxDbFromGRanges：通过GRanges对象制作TxDb
makeTxDbFromGFF：通过解析GFF文件制作TxDb

制作TxDb方法示例：

  • 如用人的参考基因信息来做注释，从UCSC生成TxDb:
    1. library(GenomicFeatures)
    2. hg19.refseq.db <- makeTxDbFromUCSC(genome="hg19", table="refGene")
  • 用GFF文件做裂殖酵母的注释
    1. download.file("ftp://ftp.ebi.ac.uk/pub/databases/pombase/pombe/Chromosome_Dumps/gff3/schizosaccharomyces_pombe.chr.gff3", "schizosaccharomyces_pombe.chr.gff3")require(GenomicFeatures)
    2. spombe <- makeTxDbFromGFF("schizosaccharomyces_pombe.chr.gff3")

具体步骤如下：

第1步：下载安装ChIPseeker注释相关的包

从Bioconductor直接下载，或从github安装最新版本

  1. source ("https://bioconductor.org/biocLite.R")
  2. biocLite("ChIPseeker")
  3. # 下载人的基因和lincRNA的TxDb对象
  4. biocLite("org.Hs.eg.db")
  5. biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
  6. biocLite("TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts")
  7. biocLite("clusterProfiler")
  8. #载入各种包
  9. library("ChIPseeker")
  10. library(clusterProfiler)
  11. library("org.Hs.eg.db")
  12. library(TxDb.Hsapiens.UCSC.hg19.knownGene)
  13. txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
  14. library("TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts")
  15. lincRNA_txdb=TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts

第2步：读入peaks文件

函数readPeakFile读入peaks文件

  1. nanog <- readPeakFile("./idr_out.bed/nanog_idr-bed")
  2. pou5f1 <- readPeakFile("./idr_out.bed/pou5f1_idr-bed")

第3步：注释peaks

peaks的注释是用的annotatePeak函数，可以单独对每个peaks文件进行注释，也可以将多个peaks制作成一个list，进行比较分析和可视化。

  1. # 制作多个样本比较的list
  2. peaks <- list(Nanog=nanog,Pou5f1=pou5f1)
  3. # promotor区间范围可以自己设定
  4. promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
  5. tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)
  6. #annotatePeak传入annoDb参数,可进行基因ID转换（Entrez，ENSEMBL，SYMBOL，GENENAME）
  7. peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb,tssRegion=c(-3000, 3000), verbose=FALSE,addFlankGeneInfo=TRUE, flankDistance=5000,annoDb="org.Hs.eg.db")

annotatePeak传入annoDb参数,即可进行基因ID转换，将Entrez ID转化为ENSEMBL，SYMBOL，GENENAME，peakAnnoList的结果如下：

  1. seqnames start end width strand annotation geneChr geneStart geneEnd geneLength geneStrand geneId transcriptId distanceToTSS ENSEMBL SYMBOL GENENAME flank_txIds flank_geneIds flank_gene_distances
  2. 5 chr3 196625522 196625873 352 * Intron (uc003fwz.4/205564, intron 2 of 9) 3 196594727 196661584 66858 1 205564 uc011bty.2 30795 ENSG00000119231 SENP5 SUMO specific peptidase 5 uc003fwz.4;uc011bty.2 205564;205564 0;0

第4步：可视化

提供了多种可视化方法，如plotAnnoBar(),vennpie(),plotAnnoPie(),plotDistToTSS()等，下面展示了两个样本在基因组特征区域的分布以及转录因子在TSS区域的结合。

  1. plotAnnoBar(peakAnnoList)
  2. plotDistToTSS(peakAnnoList,title="Distribution of transcription factor-binding loci n relative to TSS")

ATAC-Seq分析教程：用ChIPseeker对peaks进行注释和可视化-图片1

第5步：功能富集分析

提取peakAnnolist中的基因，结合clusterProfiler包对peaks内的邻近基因进行富集注释。

  1. # Create a list with genes from each sample
  2. gene = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
  3. # Run GO enrichment analysis 
  4. ego <- enrichGO(gene = entrez, 
  5.                     keytype = "ENTREZID", 
  6.                     OrgDb = org.Hs.eg.db, 
  7.                     ont = "BP", 
  8.                     pAdjustMethod = "BH", 
  9.                     qvalueCutoff = 0.05, 
  10.                     readable = TRUE)
  11. # Dotplot visualization
  12. dotplot(ego, showCategory=50)
  13. # Multiple samples KEGG analysis
  14. compKEGG <- compareCluster(geneCluster = gene, 
  15.                          fun = "enrichKEGG",
  16.                          organism = "human",
  17.                          pvalueCutoff  = 0.05, 
  18.                          pAdjustMethod = "BH")
  19. dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")

ATAC-Seq分析教程：用ChIPseeker对peaks进行注释和可视化-图片2

第6步：保存文件

  1. # Output peakAnnolist file
  2. save(peakAnnoList,file="peakAnnolist.rda")
  3. write.table(as.data.frame(peakAnnoList$Nanog),file="Nanog.PeakAnno",sep='t',quote = F)
  4. # Output results from GO analysis to a table
  5. cluster_summary <- data.frame(ego)
  6. write.csv(cluster_summary, "results/clusterProfiler_Nanog.csv")

发表评论

匿名网友

拖动滑块以完成验证
加载失败