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