ChIPseq-GEO数据挖掘
ChIPseeker系列文章已经介绍了很多内容,包括注释的方方面面,也包括强大的可视化功能(《CS6: ChIPseeker的可视化方法(中秋节的视觉饕餮)》)。
今天要介绍一下数据挖掘,从大量已有的数据来产生新的hypothesis。正如我在ChIPseeker的文章里写的:
There are increasing evidences shown that combinations of TFs are important for regulating gene expression (Perez-Pinera et al., 2013; Zhu et al., 2008). However, systematically identification of TF interactions by ChIP-seq is still not available. Even if a specific TF binding is essential for a particular regulation was known, we do not have prior knowledge of all its co-factors. There are no systematic strategies available to identified un-known co-factors by ChIP- seq.
并没有方法可以大规模地预测未知的共同调控因子,而数据挖掘就是要给我们这种预测的能力。
我当年在写ChIPseeker的时候,我有纠结是写篇Bioinformatics的application note呢,还是写篇长文灌水NAR,毕竟NAR影响因子高一点,最后还是发了Bioinformatics,因为我没钱,囧,Bioinformatics不要版面费啊。然后限于篇幅,ChIPseeker有大量可视化的函数,我在文章中一张图都没放!!!如果当时决定发NAR的话,这个数据挖掘这一块我就会写多点。
做注释在Windows上有个软件CisGenome,是Hongkai Ji课题组做的,他们还做了一个hmChIP的database:
这个hmChIP可以说是先行者,它只支持人和鼠,收集了一些GEO的数据,这个数据库做得比较早,当年GEO上的ChIPseq数据并不多,当由于搞得早,也发了一篇Bioinformatics,当然灌水之后,就不维护了,数据没更新,而基因组版本号还停留在hg18和mm8之上,所以即使它还能用,我想也不会有人去用了。这个数据库做一个事情,你可以上传自己的数据,然后它会和数据库中的数据进行比较,给你报道出你的数据和数据库中的数据overlap了多少,仅此而已。
可以说这是一个雏形,我一直在想overlap了多少,不是很能说明问题,多少算多?能不能算个p值出来,overlap了这么多,随机情况下是不太可能产生的。思来想去,没办法直接算啊,于是求助于permutation,没有枪没有炮,我们自己造。
写了一个shuffle函数,也就是传说中的洗牌,你给基因位置信息,随机给你返回等长的片段,但位置改变了。
p <- GRanges(seqnames=c("chr1", "chr3"),
ranges=IRanges(start=c(1, 100), end=c(50, 130)))
shuffle(p, TxDb=txdb)
## GRanges object with 2 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 [239651460, 239651509] *
## [2] chr3 [163562934, 163562964] *
## -------
## seqinfo: 2 sequences from an unspecified genome; no seqlengths
有了这个函数,我们可以对数据洗牌洗它一千次、一万次,那么就可以计算一千一万个overlap,来形成背景分布,实际数据的overlap在这一背景分布中产生的概率就可以计算出来。
enrichPeakOverlap(queryPeak = files[[5]],
targetPeak = unlist(files[1:4]),
TxDb = txdb,
pAdjustMethod = "BH",
nShuffle = 50,
chainFile = NULL,
verbose = FALSE)
## qSample
## ARmo_0M GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz
## ARmo_1nM GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz
## ARmo_100nM GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz
## CBX6_BF GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz
## tSample qLen tLen N_OL
## ARmo_0M GSM1174480_ARmo_0M_peaks.bed.gz 1663 812 0
## ARmo_1nM GSM1174481_ARmo_1nM_peaks.bed.gz 1663 2296 8
## ARmo_100nM GSM1174482_ARmo_100nM_peaks.bed.gz 1663 1359 3
## CBX6_BF GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz 1663 1331 968
## pvalue p.adjust
## ARmo_0M 0.90196078 0.9019608
## ARmo_1nM 0.25490196 0.4444444
## ARmo_100nM 0.33333333 0.4444444
## CBX6_BF 0.05882353 0.2352941
函数里有一个chainFile
的参数,如果你提供了,数据的比较可以跨基因组版本,甚至于可以跨物种比较。
有了这个函数,那么我再打通GEO,让用户可以检索下载GEO的数据,如此一来,你就可以用自己的数据去和GEO的数据做比较,或者就可能挖出overlap非常significant但没被报道到的co-factor呢!
ChIPseeker已经收集了超过90个物种,超过2万个BED文件的信息。你可以用getGEOgenomeVersion列出所有支持的基因组及相应的BED文件数目。
getGEOgenomeVersion()
## organism genomeVersion Freq
## 1 Anolis carolinensis anoCar2 5
## 2 Bos taurus bosTau4 2
## 3 Bos taurus bosTau6 24
## 4 Bos taurus bosTau7 2
## 5 Caenorhabditis elegans ce10 4
## 6 Caenorhabditis elegans ce6 64
## 7 Danio rerio danRer6 6
## 8 Danio rerio danRer7 61
## 9 Drosophila melanogaster dm3 502
## 10 Gallus gallus galGal3 32
## 11 Gallus gallus galGal4 15
## 12 Homo sapiens hg18 2512
## 13 Homo sapiens hg19 6876
## 14 Homo sapiens hg38 43
## 15 Mus musculus mm10 214
## 16 Mus musculus mm8 507
## 17 Mus musculus mm9 6289
## 18 Monodelphis domestica monDom5 8
## 19 Pan troglodytes panTro3 48
## 20 Pan troglodytes panTro4 42
## 21 Macaca mulatta rheMac2 81
## 22 Macaca mulatta rheMac3 31
## 23 Rattus norvegicus rn5 3
## 24 Saccharomyces cerevisiae sacCer2 141
## 25 Saccharomyces cerevisiae sacCer3 227
## 26 Sus scrofa susScr2 17
## 27 Xenopus (Silurana) tropicalis xenTro3 3
比如说斑马鱼,文件也就几十个,你大可全部下载下来。
downloadGEObedFiles(genome="danRer7", destDir="danRer7")
而比如人鼠这些明显特种,实在太多,一般来说全部下载也不太现实。ChIPseeker可以给你列出信息,连pubmed ID都给出来了,也方便翻阅文献,如果simplify=FALSE
的话,还会给出protocal和data processing等信息哦。
hg19 <- getGEOInfo(genome="hg19", simplify=TRUE)
head(hg19)
## series_id gsm organism
## 111 GSE16256 GSM521889 Homo sapiens
## 112 GSE16256 GSM521887 Homo sapiens
## 113 GSE16256 GSM521883 Homo sapiens
## 114 GSE16256 GSM1010966 Homo sapiens
## 115 GSE16256 GSM896166 Homo sapiens
## 116 GSE16256 GSM910577 Homo sapiens
## title
## 111 Reference Epigenome: ChIP-Seq Analysis of H3K27me3 in IMR90 Cells; renlab.H3K27me3.IMR90-02.01
## 112 Reference Epigenome: ChIP-Seq Analysis of H3K27ac in IMR90 Cells; renlab.H3K27ac.IMR90-03.01
## 113 Reference Epigenome: ChIP-Seq Analysis of H3K14ac in IMR90 Cells; renlab.H3K14ac.IMR90-02.01
## 114 polyA RNA sequencing of STL003 Pancreas Cultured Cells; polyA-RNA-seq_STL003PA_r1a
## 115 Reference Epigenome: ChIP-Seq Analysis of H4K8ac in hESC H1 Cells; renlab.H4K8ac.hESC.H1.01.01
## 116 Reference Epigenome: ChIP-Seq Analysis of H3K4me1 in Human Spleen Tissue; renlab.H3K4me1.STL003SX.01.01
## supplementary_file
## 111 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM521nnn/GSM521889/suppl/GSM521889_UCSD.IMR90.H3K27me3.SK05.bed.gz
## 112 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM521nnn/GSM521887/suppl/GSM521887_UCSD.IMR90.H3K27ac.YL58.bed.gz
## 113 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM521nnn/GSM521883/suppl/GSM521883_UCSD.IMR90.H3K14ac.SK17.bed.gz
## 114 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1010nnn/GSM1010966/suppl/GSM1010966_UCSD.Pancreas.mRNA-Seq.STL003.bed.gz
## 115 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM896nnn/GSM896166/suppl/GSM896166_UCSD.H1.H4K8ac.AY17.bed.gz
## 116 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM910nnn/GSM910577/suppl/GSM910577_UCSD.Spleen.H3K4me1.STL003.bed.gz
## genomeVersion pubmed_id
## 111 hg19 19829295
## 112 hg19 19829295
## 113 hg19 19829295
## 114 hg19 19829295
## 115 hg19 19829295
## 116 hg19 19829295
然后经过你的过滤,你就可以选择性地下载部分文件来分析了。
gsm <- hg19$gsm[sample(nrow(hg19), 10)]
downloadGSMbedFiles(gsm, destDir="hg19")