2020-06-17 RIPSeeker: 用于从RIP-seq
RIPSeeker
clp
17 June, 2020
RIPSeeker: 用于从RIP-seq实验中识别蛋白质相关转录本的统计R包
RIPSeeker使用具有负二项分布概率的双态HMM(two-state HMM with negative binomial emission probability)从RIP-Seq比对中推断和区分RIP峰值。虽然RIPSeeker是专门为RIP-seq数据分析量身定做的，但它也提供了一套集成在这个独立软件包中的生物信息学工具，全面解决了从比对后处理到可视化和注释的各种问题。此外，还提供了一种基于规则的方法，作为一个名为[rulebaseRIPSeek](http://127.0.0.1:18829/help/library/RIPSeeker/help/rulebaseRIPSeek)的附加函数，用户可以根据给定单端或双端比对的自动检索在线Ensembl注释，获得RIP(和对照)中基因/转录本表达的RPKM/FPKM(和fold-change)。
[ripSeek](http://127.0.0.1:18829/help/library/RIPSeeker/help/ripSeek)的前端主功能对于大多数应用程序来说已经足够了。该函数将比对获得文件(BAM/BED/SAM)的路径作为唯一必需的参数，并输出预测的RIP区域。可选参数，用户可以通过'cNAME'指示控制第一文件参数列表中的哪些文件以实现经验错误发现率(eFDR)计算。如果设置了参数'biomaRt_dataset'和/或'goAnno'，则ripSeek将返回与RIP预测的基因组上下游相对应的带注释的RIP预测和GO富集通路。用户也可以通过logOddCutoff、pvalCutoff、pvalAdjCutoff、eFDRCutoff指定统计显著性分数的阈值。
下载并安装包
rm(list = ls())
options(stringsAsFactors = F)
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="http://mirrors.cloud.tencent.com/CRAN/"))
options(download.file.method = 'libcurl')
options(url.method='libcurl')
#BiocManager::install("RIPSeeker")
#BiocManager::install("RIPSeekerData")
library(RIPSeeker)
#> Warning: package 'Rsamtools' was built under R version 3.6.2
#?RIPSeeker
ls("package:RIPSeeker")
#> [1] "addDummyProb" "addPseudoAlignment" "annotateRIP"
#> [4] "binCount" "combineAlignGals" "combineRIP"
#> [7] "computeLogOdd" "computeRPKM" "disambiguateMultihits"
#> [10] "empiricalFDR" "evalBinSize" "exportGRanges"
#> [13] "galp2gal" "getAlignGal" "logScoreWithControl"
#> [16] "logScoreWithoutControl" "mainSeek" "mainSeekSingleChrom"
#> [19] "nbh" "nbh_chk" "nbh_em"
#> [22] "nbh_gen" "nbh_init" "nbh_vit"
#> [25] "nbh.GRanges" "nbh.integer" "nbm_chk"
#> [28] "nbm_em" "plotCoverage" "plotStrandedCoverage"
#> [31] "randindx" "ripSeek" "rulebaseRIPSeek"
#> [34] "scoreMergedBins" "seekRIP" "selectBinSize"
#> [37] "statdis" "viewRIP"
示例
library(RIPSeekerData)
extdata.dir <- system.file("extdata", package="RIPSeekerData")
bamFiles <- list.files(extdata.dir, "\\.bam$",
recursive=TRUE, full.names=TRUE)
bamFiles <- grep("PRC2/", bamFiles, value=TRUE)
查看bam文件
# Retrieve system files
# OR change it to the extdata.dir from the code chunk above
# to get RIP predictions on the full alignment data
extdata.dir <- system.file("extdata", package="RIPSeeker")
bamFiles <- list.files(extdata.dir, "\\.bam$",
recursive=TRUE, full.names=TRUE)
bamFiles <- grep("PRC2", bamFiles, value=TRUE)
# RIP alignment for Ezh2 in mESC
ripGal <- combineAlignGals(grep(pattern="SRR039214", bamFiles, value=TRUE, invert=TRUE),
reverseComplement=TRUE, genomeBuild="mm9")
# Control RIP alignments for mutant Ezh2 -/- mESC
ctlGal <- combineAlignGals(grep(pattern="SRR039214", bamFiles, value=TRUE, invert=FALSE),
reverseComplement=TRUE, genomeBuild="mm9")
ripGal
#> GAlignments object with 31054 alignments and 1 metadata column:
#> seqnames strand cigar qwidth start end
#> <Rle> <Rle> <character> <integer> <integer> <integer>
#> SRR039210.4322179 chrX + 36M 36 3000326 3000361
#> SRR039210.5524106 chrX - 36M 36 3001293 3001328
#> SRR039210.5069294 chrX - 36M 36 3002790 3002825
#> SRR039210.1476279 chrX - 36M 36 3016328 3016363
#> SRR039210.711491 chrX + 36M 36 3021161 3021196
#> ... ... ... ... ... ... ...
#> SRR039211.1427139 chrX - 36M 36 166443604 166443639
#> SRR039211.447965 chrX - 36M 36 166445241 166445276
#> SRR039211.1352670 chrX - 36M 36 166446255 166446290
#> SRR039211.918096 chrX - 36M 36 166446315 166446350
#> SRR039211.67451 chrX - 36M 36 166446621 166446656
#> width njunc | uniqueHit
#> <integer> <integer> | <logical>
#> SRR039210.4322179 36 0 | FALSE
#> SRR039210.5524106 36 0 | FALSE
#> SRR039210.5069294 36 0 | FALSE
#> SRR039210.1476279 36 0 | TRUE
#> SRR039210.711491 36 0 | FALSE
#> ... ... ... . ...
#> SRR039211.1427139 36 0 | FALSE
#> SRR039211.447965 36 0 | FALSE
#> SRR039211.1352670 36 0 | FALSE
#> SRR039211.918096 36 0 | FALSE
#> SRR039211.67451 36 0 | TRUE
#> -------
#> seqinfo: 22 sequences from mm9 genome
ctlGal
#> GAlignments object with 8276 alignments and 1 metadata column:
#> seqnames strand cigar qwidth start end
#> <Rle> <Rle> <character> <integer> <integer> <integer>
#> SRR039214.4949641 chrX + 36M 36 3028539 3028574
#> SRR039214.5310910 chrX - 35M 35 3039791 3039825
#> SRR039214.4625677 chrX - 36M 36 3084567 3084602
#> SRR039214.1854227 chrX + 36M 36 3139109 3139144
#> SRR039214.5753635 chrX + 36M 36 3139131 3139166
#> ... ... ... ... ... ... ...
#> SRR039214.904524 chrX - 32M 32 166445788 166445819
#> SRR039214.2473565 chrX - 21M 21 166445960 166445980
#> SRR039214.576581 chrX - 36M 36 166446074 166446109
#> SRR039214.3756853 chrX - 36M 36 166446781 166446816
#> SRR039214.2347055 chrX + 36M 36 166650104 166650139
#> width njunc | uniqueHit
#> <integer> <integer> | <logical>
#> SRR039214.4949641 36 0 | FALSE
#> SRR039214.5310910 35 0 | FALSE
#> SRR039214.4625677 36 0 | FALSE
#> SRR039214.1854227 36 0 | FALSE
#> SRR039214.5753635 36 0 | FALSE
#> ... ... ... . ...
#> SRR039214.904524 32 0 | FALSE
#> SRR039214.2473565 21 0 | FALSE
#> SRR039214.576581 36 0 | FALSE
#> SRR039214.3756853 36 0 | FALSE
#> SRR039214.2347055 36 0 | TRUE
#> -------
#> seqinfo: 22 sequences from mm9 genome
我们都知道RIP（RNA Binding Protein Immunoprecipitation）是RNA结合蛋白免疫共沉淀。顾名思义，是进行RNA结合蛋白（RBP）与RNA之间的研究。利用免疫共沉淀的方法将蛋白-RNA复合物拉下来，将RNA结合蛋白所结合的RNA提取出来，并对其进行高通量测序，鉴定出基因和蛋白结合位点，那么将会对RNA与蛋白质相互作用在调节mRNA和非编码RNA功能方面提供研究的途径。理想情况下，通过RIP对目的蛋白结合的RNA进行富集后，目的蛋白结合的RNA或者目的蛋白在RNA上的结合区域，在对应的参考基因组位置上，测序reads的覆盖度会显著升高，相对其他非结合区域形成明显的peak。所以，通过测序数据检测peak，可以获知目的蛋白结合的RNA，以及可能的结合的区域等信息。
RIPSeeker是一个基于隐马尔可夫模型进行从头分析RIP peaks预测的免费开源Bioconductor R包，有着较高的敏感性和特异性。RIPSeeker区分正负链，可以鉴定链特异性的peak区域，有利于链特异性非编码RNA的鉴定。并且，RIPSeeker不局限于鉴定狭义上的peak区域，其适用于检测更大范围的峰值分布，以鉴定不同长度范围分布的整条结合转录本。所以，这也是为什么进行RIP-seq数据分析的时候，我们选择RIPSeeker的主要原因。不同于其他的peak calling方式，可以说RIPSeeker是为RIP-seq量身定做的。
image
虽然RIP-seq实验和ChIP-seq以及RNA-seq有着相似之处，但是RIP-seq有一个最根本不同的目标，就是发现目的结合蛋白相关的转录本。如上图所示，展示了3种不同模式下的比较。在ChIP-seq中，目的蛋白结合的双链DNA被抗体拉下来，然后进行高通量测序。由于reads通常短于双链DNA片段，测序数据会显示出一个特性，即真正的结合位点，会在与片段长度接近的一个距离d上，产生正负链的一个对称峰，如上图（a）所示。而由于RNA转录本是单链的，在ChIP-seq中观察到的双链DNA的双峰性质就不适用于RNA-seq和RIP-seq。因此，也就不适合用通过查找这种双峰的模式来进行RIP-seq的分析。并且，针对RNA来说，我们还需要考虑到剪切比对的形式，因为剪切事件的存在，也不能跟双链DNA的模式相同，比对reads不能直接沿着基因组延伸。
Peak calling 方法比较
RIPSeeker的作者，将RIPSeeker与各种高通量测序分析中流行的其他算法进行了比较。作者选择了三种ChIP-seq算法，包括MACS、QuEST和HPeak；两种RNA-seq算法Cufflinks+Cuffdiff和Rulebased算法和一种PAR-CLIP算法PARalyzer。
image
通过各方法的比较，发现对于相同的数据来源，所得到的峰的个数差别很大。这可能是由于不同方法之间使用的打分规则和peak长度不同所致。MACS在ENCODE数据上表现出来是检测的峰数量更多，这是由于别的方法在很大程度上会将接近的peak峰区段连接成一段较长的连续区间。针对RIPSeeker和个别其他方法来说，针对两个生物学重复，它们鉴定得到的peak重复度一般高于50%(如下图柱状图灰色部分)。
image
用同一目的蛋白富集的数据，使用两种不同的对照进行分析，来研究RIPSeeker的鲁棒性。作者针对同一细胞系的同一RIP处理组，分别将非特异性抗体富集的结果（T7Tag）和RIP input作为阴性对照，得到约60-80%的重合率。这也意味着，两种不同的对照条件可以交换使用。
image
为了比较不同方法之间结果异质性，作者使用同一数据集，对任意两种方法的结果重叠情况进行了比较。结果显示，RIPSeeker和其他方法之间有比较好的结果重叠（一般>50%）。
image
image
基于敏感度和特异性的ROC评价显示，RIPSeeker在大多数测试中占主导地位。RIPSeeker在识别信号峰方面的敏感度和特异性优于其他方法。
image
image
image
image
为了证明RIPSeeker软件包的实用性，作者将RIPSeeker和其他发表的6个工具，作用于3个RIP-seq数据集和2个PAR-CLIP数据集。基于受试者曲线，RIPSeeker表现出优越的敏感性和特异性。来自RIPSeeker鉴定所得的peaks，在后续的生物学研究中，可通过特定基因的富集情况、已发表的一些motif结果、与目的蛋白相关的典型转录本等，被进一步证实。
参考文献：
[1]. Li Yue,Zhao Dorothy Yanling,Greenblatt Jack F et al. RIPSeeker: a statistical package for identifying protein-associated transcripts from RIP-seq experiments.[J] .Nucleic Acids Res., 2013, 41: e94.
[2]. Zhang,Y., Liu,T., Meyer,C.A., Eeckhoute,J., Johnson,D.S.,Bernstein,B.E., Nusbaum,C., Myers,R.M., Brown,M., Li,W. et al.(2008) Model-based analysis of ChIP-Seq (MACS). Genome Biol.,9, R137
[3]. Valouev,A., Johnson,D.S., Sundquist,A., Medina,C., Anton,E.,Batzoglou,S., Myers,R.M. and Sidow,A. (2008) Genome-wide analysis of transcription factor binding sites based on ChIP-Seq data. Nat. Methods, 5, 829–834.
[4]. Qin,Z.S., Yu,J., Shen,J., Maher,C.A., Hu,M., KalyanaSundaram,S., Yu,J. and Chinnaiyan,A.M. (2010) HPeak: an HMM-based algorithm for defining read-enriched regions in ChIP-Seq data. BMC Bioinformatics, 11, 369.
[5]. Trapnell,C., Williams,B.A., Pertea,G., Mortazavi,A., Kwan,G., van Baren,M.J., Salzberg,S.L., Wold,B.J. and Pachter,L. (2010) Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol., 28, 516–520
[6]. Roberts,A., Goff,L., Pertea,G., Kim,D., Kelley,D.R., Pimentel,H., Salzberg,S.L., Rinn,J.L., Pachter,L. and Trapnell,C. (2012) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc., 7, 562–578
[7]. Corcoran,D.L., Georgiev,S., Mukherjee,N., Gottwein,E., Skalsky,R.L., Keene,J.D. and Ohler,U. (2011) PARalyzer: definition of RNA binding sites from PAR-CLIP short-read sequence data. Genome Biol., 12, R79
[8]. Zhao,J., Ohsumi,T.K., Kung,J.T., Ogawa,Y., Grau,D.J., Sarma,K., Song,J.J., Kingston,R.E., Borowsky,M. and Lee,J.T. (2010) Genome-wide identification of polycomb-Associated RNAs by RIP-seq. Mol. Cell, 40, 939–953
