R-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")