100天生信-Day4

最近在做類ChIP-seq，看了一下使用ChIPseeker注釋peaks的方法，保存一下代碼：

install.packages('R.utils')
BiocManager::install("GenomicFeatures")
BiocManager::install("ChIPseeker")

library(R.utils)
library(GenomicFeatures)
library(ChIPseeker)
library(ggplot2)

##制作TxDb注釋文件
# 下載gff3
download.file('ftp://ftp.ensemblgenomes.org/pub/plants/release-51/gff3/triticum_aestivum/Triticum_aestivum.IWGSC.51.gff3.gz',destfile = 'Ta.gff3.gz')
# 解壓
R.utils::gunzip('Ta.gff3.gz')
# 制作
wheat2 <- makeTxDbFromGFF("Ta.gff3")

## 讀入BED文件
setwd('path')
Ta_peaks_file <- readPeakFile('Ta_peaks.bed')

## 可視化peaks
# 染色體上分布
covplot(Ta_peaks_file)
# TSS上下游折線圖
plotAvgProf2(Ta_peaks_file, TxDb=wheat2, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency",
             conf = 0.95, resample = 1000)
# TSS上下游熱圖
peakHeatmap(Ta_peaks_file, TxDb=wheat2, 
            upstream=3000, downstream=3000, 
            color=rainbow(length(Ta_peaks_file)))

## 注釋peaks
peakAnnoList <- annotatePeak(Ta_peaks_file, tssRegion=c(-2500,2500), TxDb=wheat2, flankDistance=5000)
# 轉(zhuǎn)換為dataframe方便查看
peakAnnoList_frame <- as.data.frame(peakAnnoList)
# 保存dataframe
write.csv(peakAnnoList_frame, file = 'Ta_peakAnnoList.csv')

## 可視化注釋情況
# 條形圖
plotDistToTSS(peakAnnoList)
# 餅圖
plotAnnoPie(peakAnnoList)

參考教程：https://www.jieandze1314.com/post/cnposts/190/
http://www.itdecent.cn/p/a7b6ce208f98