RNA-seq是目前最常用的高通量测序技术之一。其中,有一种很常见的用途,就是测量在不同条件的实验组之间哪些基因的表达是有显著差异的,或许还会对这些差异表达的基因相关的通路或者基因集进行进一步的分析。对于学习生物信息的同学来说,RNA-seq数据分析几乎是必备技能了;而对于做湿实验的朋友们,了解一些RNA-seq相关的知识也有助于更好地理解合作者或者公司给出的数据分析结果。 很多人可能以为RNA-seq数据分析很难学,然而实际上并不是这样的。在limma-voom的作者Charity Law和她的同事们写的RNA-seq数据差异表达分析流程的教程中,RNA-seq数据分析被形容为“如同数1-2-3一样简单”。这篇教程中的所有分析都是在R语言中进行的,而且用到的所有包都是完全免费的。不仅如此,这篇教程最近由作者实验室的中国学生(也就是笔者本人2333)翻译成了中文 ,更为广大同胞免去了看英文教程时难以理解的痛苦!上网搜索“RNAseq123”并进入其Bioconductor包页面,即可查看教程中英文原文和其中所有步骤的代码。 这篇推文里,我们缩减了一下教程原文,简单地介绍了RNA-seq差异表达分析的数据整合、数据预处理、差异表达分析、热图和MD图的绘制以及camera基因集检验。如果想要进一步学习,更推荐大家查看教程的原文(戳“阅读原文”即可访问)。如果在学习的过程中遇到任何问题,都可以在Bioconductor support论坛提问。 1、背景介绍这篇教程中用到的数据(GEO序列号: GSE63310)由雌性处女小鼠的三种乳腺细胞类型组成,即基底(basal)、管腔祖细胞(liminal progenitor, LP)和成熟管腔(mature luminal, ML),每组各设三个生物学重复。 文中所描述的数据分析是从已经经过比对和基因水平计数的表达计数(counts)数据起步的。如果你手头想要分析的数据是原始的序列也不要慌,作者推荐使用Rsubread包的 首先载入分析需要用到的软件包。 library(limma) library(Glimma) library(edgeR) library(Mus.musculus) 2、数据整合使用edgeR的 files <- c('gsm1545535_10_6_5_11.txt',="" 'gsm1545536_9_6_5_11.txt',="">-> 组织样品信息,将样品分组、处理、批次等与实验设计有关的样品信息与计数矩阵的列关联。 samplenames <- substring(colnames(x),="" 12,="">-> colnames(x) <->-> group <- as.factor(c('lp',="" 'ml',="" 'basal',="" 'basal',="" 'ml',="" 'lp',="">-> 'Basal', 'ML', 'LP')) x$samples$group <->-> lane <- as.factor(rep(c('l004','l006','l008'),="">-> x$samples$lane <->-> 使用Mus.musculus包,利用数据集中的Entrez基因ID来检索相对应的基因符号和染色体信息,将这些基因注释信息存储到DGEList对象中的第二个数据框 geneid <->-> 为了处理重复的基因ID,选取其中第一次出现的染色体来代表具有重复注释的基因。 genes <->-> 将基因注释的数据框加入数据对象,数据即被整洁地整理入一个DGEList对象中,它包含原始计数数据和相关的样品信息和基因注释。 x$genes <->-> 3、数据预处理更深的测序总会产生更多的序列。使用edgeR中的 cpm <->-> lcpm <- cpm(x,="" log="TRUE," prior.count="">-> 在任何样本中都没有足够高的表达量的基因应该从下游分析中过滤掉。edgeR包中的 keep.exprs <- filterbyexpr(x,="" group="">-> 下游分析假设所有样品表达值的范围和分布都相似,因此我们需要进行归一化来确保整个实验中每个样本的表达分布都相似。使用edgeR中的 x <- calcnormfactors(x,="" method='TMM'>-> 用limma中的 lcpm <- cpm(x,="" log="">-> 在这个数据集中,可以看出样本在维度1和2能很好地按照实验分组聚类,随后在维度3按照测序道(样品批次)分离(如下图所示)。 此外,也可以使用Glimma包的 glMDSPlot(lcpm, labels=paste(group, lane, sep='_'), groups=x$samples[,c(2,5)], launch=FALSE) 4、差异表达分析在此研究中,我们想知道哪些基因在我们研究的三组细胞之间以不同水平表达。在我们的分析中,假设基础数据是正态分布的,为其拟合一个线性模型。在此之前,需要创建一个包含细胞类型以及测序泳道(批次)信息的设计矩阵。 design <->-> 对于一个给定的实验,通常有几种等价的方法可以创建一个合适的设计矩阵。比如说, 用于细胞群之间成对比较的对比可以在limma中用 contr.matrix <->-> BasalvsLP = Basal-LP, BasalvsML = Basal - ML, LPvsML = LP - ML, levels = colnames(design)) limma的线性模型方法的核心优势之一便是其可以适应不同的实验复杂程度。无论这篇教程中关于细胞类型和批次的简单实验设计,还是更复杂的因子设计和含有交互作用项的模型,都能够得到相对简单的处理。 对于RNA-seq原始计数或其log-CPM值,方差与均值并不独立,而使用负二项分布来模拟计数的方法假设均值与方差间具有二次的关系。在limma中,假设log-CPM值符合正态分布,并使用由 下图左侧展示了这个数据集log-CPM值的均值-方差关系。通常而言,voom图中均值随着方差增加而递减。生物学差异高的实验通常会有更平坦的趋势,其方差值在高表达处稳定。生物学差异低的实验更倾向于急剧下降的趋势。 此外,如果对于低表达基因的过滤不够充分,在图上表达低的一端可以观察到方差水平的下降。如果观察到了这种情况,应当回到最初的过滤步骤并提高用于该数据集的表达阈值。 par(mfrow=c(1,2)) limma的线性建模使用 为快速查看差异表达水平,以表格列出各对比中显著上调和下调的基因数量。 summary(decideTests(efit)) ## BasalvsLP BasalvsML LPvsML 如果需要检测差异表达基因之间的表达量差异倍数是否大于某一阈值(比如本文中的 tfit <- treat(vfit,="" lfc="">-> dt <->-> summary(dt) ## BasalvsLP BasalvsML LPvsML 在多种比较中皆差异表达的基因可以从 de.common <- which(dt[,1]!="0" &="" dt[,2]!="">-> vennDiagram(dt[,1:2], circle.col=c('turquoise', 'salmon')) write.fit(tfit, dt, file='results.txt') 使用 basal.vs.lp <- toptreat(tfit,="" coef="1," n="">-> ## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val ## 12759 12759 Clu chr14 -5.46 8.86 -33.6 1.72e-10 1.71e-06 ## 53624 53624 Cldn7 chr11 -5.53 6.30 -32.0 2.58e-10 1.71e-06 ## 242505 242505 Rasef chr4 -5.94 5.12 -31.3 3.08e-10 1.71e-06 ## 67451 67451 Pkp2 chr16 -5.74 4.42 -29.9 4.58e-10 1.74e-06 ## 228543 228543 Rhov chr2 -6.26 5.49 -29.1 5.78e-10 1.74e-06 ## 70350 70350 Basp1 chr15 -6.08 5.25 -28.3 7.27e-10 1.74e-06 basal.vs.ml <- toptreat(tfit,="" coef="2," n="">-> ## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val ## 242505 242505 Rasef chr4 -6.53 5.12 -35.1 1.23e-10 1.24e-06 ## 53624 53624 Cldn7 chr11 -5.50 6.30 -31.7 2.77e-10 1.24e-06 ## 12521 12521 Cd82 chr2 -4.69 7.07 -31.4 2.91e-10 1.24e-06 ## 20661 20661 Sort1 chr3 -4.93 6.70 -30.7 3.56e-10 1.24e-06 ## 71740 71740 Nectin4 chr1 -5.58 5.16 -30.6 3.72e-10 1.24e-06 ## 12759 12759 Clu chr14 -4.69 8.86 -28.0 7.69e-10 1.48e-0 5、为结果绘图
plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1], Glimma的 glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1], side.main='ENTREZID', counts=lcpm, groups=group, launch=FALSE) 上方指令生成的均值-差异图可以在这里查看:https:///2MGbzKF 使用gplots包的 library(gplots) 7、使用camera的基因集检验从http://bioinf./software/MSigDB/以RData对象格式下载Broad Institute的MSigDB c2中的适应小鼠的c2基因表达特征,并运行camera。此外,对于人类和小鼠,来自MSigDB的其他有用的基因集也可从此网站获取。 load(system.file('extdata', 'mouse_c2_v5p1.rda', package = 'RNAseq123')) idx <- ids2indices(mm.c2,id="">-> cam.BasalvsLP <- camera(v,idx,design,contrast="">-> cam.BasalvsML <- camera(v,idx,design,contrast="">-> cam.LPvsML <- camera(v,idx,design,contrast="">->
barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP, limma中也有其他的基因集检验,比如mroast的自包含检验。简单来说,camera更适用于在大型数据库中搜寻具有意义的基因集,而mroast测试的是已经确定有意义的一个或少个基因集的显著性。 此外,我们也推荐使用EGSEA进行基因集检验,结合12种优秀的的基因集测试算法的结果来获得生物学结果的整体排名,并且可以生成结果报告网页。如果感兴趣的话,可以搜索EGSEA123查看其使用教程(目前只有英文版,但我们也有在今后将其翻译成中文的打算)。 8、参考文献 Bioconductor Core Team. 2016a. Homo.sapiens: Annotation package for the Homo.sapiens object. https:///packages/release/data/annotation/html/Homo.sapiens.html. ———. 2016b. Mus.musculus: Annotation package for the Mus.musculus object. https:///packages/release/data/annotation/html/Mus.musculus.html. Durinck, S., Y. Moreau, A. Kasprzyk, S. Davis, B. De Moor, A. Brazma, and W. Huber. 2005. “BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis.” Bioinformatics 21:3439–40. Durinck, S., P. Spellman, E. Birney, and W. Huber. 2009. “Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt.” Nature Protocols 4:1184–91. Huber, W., V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2):115–21. http://www./nmeth/journal/v12/n2/full/nmeth.3252.html. Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. “Voom: Precision Weights Unlock Linear Model Analysis Tools for RNA-seq Read Counts.” Genome Biology 15:R29. Liao, Y., G. K. Smyth, and W. Shi. 2013. “The Subread Aligner: Fast, Accurate and Scalable Read Mapping by Seed-and-Vote.” Nucleic Acids Res 41 (10):e108. ———. 2014. “featureCounts: an Efficient General-Purpose Program for Assigning Sequence Reads to Genomic Features.” Bioinformatics 30 (7):923–30. Lim, E., D. Wu, B. Pal, T. Bouras, M. L. Asselin-Labat, F. Vaillant, H. Yagita, G. J. Lindeman, G. K. Smyth, and J. E. Visvader. 2010. “Transcriptome analyses of mouse and human mammary cell subpopulations reveal multiple conserved genes and pathways.” Breast Cancer Research 12 (2):R21. Liu, R., K. Chen, N. Jansz, M. E. Blewitt, and M. E. Ritchie. 2016. “Transcriptional Profiling of the Epigenetic Regulator Smchd1.” Genomics Data 7:144–7. Liu, R., A. Z. Holik, S. Su, N. Jansz, K. Chen, H. S. Leong, M. E. Blewitt, M. L. Asselin-Labat, G. K. Smyth, and M. E. Ritchie. 2015. “Why weight? Combining voom with estimates of sample quality improves power in RNA-seq analyses.” Nucleic Acids Res 43:e97. McCarthy, D. J., and G. K. Smyth. 2009. “Testing significance relative to a fold-change threshold is a TREAT.” Bioinformatics 25:765–71. Ritchie, M. E., B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth. 2015. “limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Res 43 (7):e47. Robinson, M. D., D. J. McCarthy, and G. K. Smyth. 2010. “edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26:139–40. Robinson, M. D., and A. Oshlack. 2010. “A Scaling Normalization Method for Differential Expression Analysis of RNA-seq data.” Genome Biology 11:R25. Sheridan, J. M., M. E. Ritchie, S. A. Best, K. Jiang, T. J. Beck, F. Vaillant, K. Liu, et al. 2015. “A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1.” BMC Cancer 15 (1). BioMed Central:221. Smyth, G. K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Stat Appl Genet Mol Biol 3 (1):Article 3. Subramanian, A., P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proc Natl Acad Sci U S A 102 (43):15545–50. Wu, D., E. Lim, F. Vaillant, M. L. Asselin-Labat, J. E. Visvader, and G. K. Smyth. 2010. “ROAST: rotation gene set tests for complex microarray experiments.” Bioinformatics 26 (17):2176–82. Wu, D., and G. K. Smyth. 2012. “Camera: a competitive gene set test accounting for inter-gene correlation.” Nucleic Acids Res 40 (17):e133. |
|