小王的Chipseeker注释傻瓜学习教程

参考：https://blog.csdn.net/qq_52813185/article/details/121977773

激活R环境，library一切

library(ChIPseeker)

library(ChIPpeakAnno)

library(org.Mm.eg.db)

library(GenomeInfoDb)

library(GenomicFeatures)

library(diffloop)

library(ggplot2)

library(rio)

1、明确自己使用的注释基因组参考文件

师兄说TxDb.Mmusculus.UCSC.mm10.knownGene的注释太老了，让我自己做一下。哎哎哎！新技能get

makeTxDbFromGFF：通过解析GFF文件制作TxDb

去ensmble下载release——102的gff文件

>anno_ens <- makeTxDbFromGFF("Mus_musculus.GRCm38.102.chr.gff3",organism="Mus musculus")

###注意，ensmble的注释染色体是没有chr的

> seqlevels(anno_ens)

[1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"

[15] "15" "16" "17" "18" "19" "X"  "Y"  "MT"

2、读入自己需要注释的peak文件

批量注释，把需要注释的文件都放在peaks里

>h0<-readPeakFile("H3K4me3_hCG_0h_peaks.narrowPeak")

>h4<-readPeakFile("H3K4me3_hCG_4h_peaks.narrowPeak")

>peaks<-list(h0=h0,h4=h4)

peaks <- lapply(peaks, rmchr)###去掉chr

3、注释吧！

#定义TSS

promoter_ens <- getPromoters(TxDb=anno_ens, upstream=3000, downstream=3000)

peakAnnoList <- lapply(peaks, annotatePeak, TxDb=anno_ens,tssRegion=c(-3000, 3000), verbose=FALSE,addFlankGeneInfo=TRUE, flankDistance=5000,annoDb = "org.Mm.eg.db")

write.table(peakAnnoList[1], file = 'h0.txt',sep = '\t', quote = FALSE, row.names = FALSE)

##有个问题，一个peak所在的位置，可能是一个基因的外显子，同时又是另一个基因的内含子。为了解决这个问题，使用参数genomicAnnotationPriority指定优先顺序，默认顺序是：Promoter => 5’ UTR => 3’ UTR => Exon => Intron => Downstream => Distal Intergenic