精华内容
下载资源
问答
  • 此次案例我们是要比较ARID1A这个基因突变与不突变两亚类肝癌病人哪些基因的表达存在差异(寻找差异基因)。 在之前的三篇TCGA的数据采集文章中,我们已经采集到了所需要的所有数据,具体见下: ...
  • limma对基因进行差异分析全代码,包括自动判断是否需要进行log,对0和负值的处理等。
  • limma_tutorial-源码

    2021-03-19 01:58:56
    limma_tutorial
  • limma 算法总结

    多人点赞 2021-04-29 03:16:35
    我这个人完全没办法在不了解算法的情况下用一个包,那样我会疯,可是网上居然完全找不到limma算法的讲解资料,全是介绍怎么用的。 可能我太初级了,这么简单的算法大概不值得一提…… 于是搞了一天原文,和东拼西凑...

    我这个人完全没办法在不了解算法的情况下用一个包,那样我会疯,可是网上居然完全找不到limma算法的讲解资料,全是介绍怎么用的。

    可能我太初级了,这么简单的算法大概不值得一提……
    于是搞了一天原文,和东拼西凑的统计学知识,把现在了解的先总结一下:
    limma是做差异表达的包,其算法核心在两个function上:
    lmfit()以及eBays()

    lmfit()就是multiple linear regression

    假设我们有基因表达矩阵 Y = [ y 1 , y 2 , . . . y 100 ] Y=[y_1,y_2,...y_{100}] Y=[y1,y2,...y100],其中包括50个样本和100个基因,样本包括30个病人和20个健康的
    这样病人vs健康可以写成一个长度为50的0-1向量X=[11111111….0000000], 1表示病人,0表示healthy control

    那么算法的大体思路是:

    1. 假设数据质量良好,我们只想知道病人和健康人群的差异表达,对于每一个gene i,构建 一个线性模型
      y i = a x + b y_i=ax+b yi=ax+b,然后利用最小二乘,求出这条线的a和b;这就是lmfit的作用。最小二乘法就是最小化实际的 y i y_i yi和拟合直线上的 y ^ i \hat{y}_i y^i的绝对距离。

    2. 由于只有一个变量,那么这就是个简单的线性回归。(这个时候,机器学习和统计学的区别就出来了:如果是机器学习,找到一条线,就要开始交叉验证,在测试集上验证模型的鲁棒性,然后选出最好的模型来。但这些都是基于我们有大量的样本允许我们这么折腾。而很多时候我们就没几个样本,所以生信啊,大多还是基于统计模型。另外,'交叉验证,测试集’什么的主要是为了说明模型的鲁棒性,目的是为了拿这个稳健的做预测,你说它这有什么理论依据,只能拿测试集结果的‘准确率’来说明模型有效。) 扯远了,由于这是个统计模型,模型的有效性需要各种假设检验来证明。我们的目的是找差异表达基因,拟合了线性模型,怎么看我们的x和yi有多相关呢?就需要用到 R 2 R^2 R2
      R 2 = ( v a r ( y i ) − v a r ( y f i t ) ) / v a r ( y i ) R^2=(var(y_i)-var(y_{fit}))/var(y_i) R2=(var(yi)var(yfit))/var(yi),
      yi的方差就是说,不考虑任何x的因素,单看yi上的样本离均值有多远,yfit就是考虑了x之后,我们的样本离拟合线有多远, 如果x和yi很相关,那么加上x的因素之后,yi的方差就会变小, R 2 R^2 R2就会很大 。(p.s.这里机器学习和统计的不同在于,机器学习才不在乎什么相关不相关,它一般只在乎拟合的好不好,能不能用这个线去预测别的人有没有病

    3. 然而,统计模型的麻烦之处就在于,你不仅要找到有多相关( R 2 R^2 R2),还要验证你找到的相关性有多显著。这里用F分布来求出p-value来说明,模型的显著性。(为什么是F分布?)

    4. 注意F分布的p-value是为了测试模型的显著性,我们这里只有一个变量,所以模型的显著性也就是x的显著性了。但是如果有两个或以上样本,F test是为了测试单纯的yi的方差,和用所有样本x1,x2…拟合之后的 y f i t y_{fit} yfit差异的显著性。

      t test: This test checks the significance of individual regression coefficients.
      F test: This test can be used to simultaneously check the significance of a number of regression coefficients.

    5. 可是,事实上我们可能并不关心模型的显著性,我们只是想知道哪些基因和某一个疾病(特征/自变量)相关。尤其是,对于limma,如果我们的数据存在‘batch effect’,我们还需要把batch也作为自变量输入进去,对于batch这一列自变量,我们根本就不关心它拟合的好不好,甚至还有点不想让它存在。因此我们真正需要的是对我们感兴趣的单个自变量进行的t-test检验,

    6. 对某一个自变量x做t-test,就是测试,如果我们拿掉x(线性回归参数a=0),和不拿掉x(原来的拟合结果)差异的显著性,这些结果对应了imma里的t-value和p.value

    7. 多重假设检验: 然而还没有结束,因为我们不止对一个gene做了回归,我们对全部gene都做了回归,gene i的t-test的pvalue表示了结果的假阳性率,可是所有的gene都计算,假阳性总和就会变很高,结果就不可靠了。因此我们必须要进行 multiple hypothesis test。以BH方法为例,比如100个gene,ttest之后我们觉得结果里面有假阳性,所以就对这100个gene按照pvalue从小到大排序,然后计算一个新的qvalue:
      q i = p v i ∗ ( m / k ) q_i = pv_i*(m/k) qi=pvi(m/k) 其中m是总共检测数量(m=100),k是排名
      这里我用自己的数据画了一个pvalue和adjusted pvalue的关系(m=10000),可以看到这个凸函数,放大了pvalue中较小的值
      在这里插入图片描述

    eBays()就是用经验贝叶斯调整了t-test中方差的部分

    先贴原文:
    在这里插入图片描述
    对比一下经典的t-test:
    在这里插入图片描述
    可以发现这里,就是用后验 s ~ j \tilde{s}_j s~j去替代了原始的 s j s_j sj.这才是limma的主要贡献,以上都是基本的多元线性回归。可是为什么要这样做?

    没时间了,过两天补上。。

    Reference

    http://reliawiki.org/index.php/Multiple_Linear_Regression_Analysis#:~:text=Test%20on%20Individual%20Regression%20Coefficients%20(t%20Test),-The%20t%5C%2C%5C!&text=test%20is%20used%20to%20check,may%20make%20the%20model%20worse.

    https://bioconductor.statistik.tu-dortmund.de/packages/2.11/bioc/vignettes/limma/inst/doc/usersguide.pdf

    展开全文
  • limma userguide

    2015-06-03 09:55:07
    用于基因表达谱数据分析的R语言数据包,是bioconductor数据包的子数据包
  • limma包的使用技巧

    千次阅读 2019-10-08 10:52:08
    limmar package是一个功能比较全的包,既含有cDNA芯片的RAW data输入、前处理(归一化)功能,同时也有差异化基因分析的“线性”算法(limma: Linear Models for Microarray Data),特别是对于“多因素实验...
    limmar package是一个功能比较全的包,既含有cDNA芯片的RAW data输入、前处理(归一化)功能,同时也有差异化基因分析的“线性”算法(limma: Linear Models for Microarray Data),特别是对于“多因素实验(multifactor designed experiment)”。limmar包的可扩展性非常强,单通道(one channel)或者双通道(tow channel)数据都可以分析差异基因,甚至也包括了定量PCR和RNA-seq(第一次见分析microarray的包也能处理RNA-seq)。
    limmar package是一个集大成的包,对载入数据、数据前处理(背景矫正、组内归一化和组间归一化都有很多种选择方法)、差异基因分析都有很多的选择。而且,所设计的线性回归和经验贝叶斯方法找差异基因非常值得学习。

    1. 读入样本信息
    使用函数读readTargets(file="Targets.txt", path=NULL, sep="\t", row.names=NULL, quote="\"",...)。这个函数其实是一个包装了的read.table(),读入的是样本的信息,创建的对象类似于marray包的marrayInfos和Biobase包的AnnotatedDataFrame。

    2. 读入探针密度数据 
    marray包一致,Bioconductor不能读入原始的TIFF图像文件,只能输出扫描仪输入的、转换成数字信号的文本文件。使用函数read.maimages(files=NULL, source="generic", path=NULL, ext=NULL, names=NULL, columns=NULL, other.columns=NULL, annotation=NULL, green.only=FALSE, wt.fun=NULL, verbose=TRUE, sep="\t", quote=NULL, ...)
    参数说明:files需要通过函数dir(pattern = "Mypattern")配合正则表达式筛选(规范命名很重要),同时该函数可以读入符合格式的压缩过的文件,比如*.txt.gz的文件,这极大的减小的数据储存大小;source的取值分为两类,一类是“高富帅”,比如“agilent”、“spot”等等(下表),它们是特定扫描仪器的特定输出格式;如果不幸是“屌丝”,即格式是自己规定的,可以选定source="generic",这时需要规定columns;任何cDNA文件都要有R/G/Rb/Gb四列(Mean或者Median);annotation可以规定哪些是注释列;wt.fun用于对点样点进行质量评估,取值为0表示这些点将在后续的分析中被剔除,取值位1表示需要保留,对点样点的评估依赖于图像扫描软件的程序设定,比如SPOT和GenePix软件,查看QualityWeights(现成函数或者自己写函数)。
    读入单通道数据:读入单通道数据,可以设定green.only = TRUE即可,然后对应读入columns = list(G = "Col1", Gb = "Col2")。
    读入的数据,如果是单通道,则成为EListRaw class;如果是双通道,则是RGList class。
    数据操作:
    cbind():合并数据;
    [”:分割数据;
    RGList class有的names是 "R","G","Rb","Gb","weights","printer","genes","targets","notes":  R/G/Rb/Gb分别红和绿的前景和背景噪音;weight是扫描软件的质量评估;printer是点样规则(printer layout);genes是基因注释;target是样本注释;notes是一般注释。可以通过myRGList$names进行相应的取值和赋值
    limma包的使用技巧

    3. 前处理
    cDNA芯片前处理的流程是:背景去除(background correction)--> 组内归一化(with-array normalization)--> 组间归一化(between-array normalization)。最后一步组间归一化可选,注意trade-off原则。
    3.1 背景去除: 
    backgroundCorrect(RG, method="auto", offset=0, printer=RG$printer, normexp.method="saddle", verbose=TRUE)
    RG:可以是EListRaw,也可以是RGList class;
    method:取“auto”,none”,“subtract”,"half","minimum","movingmin","edwards"或者"normexp";
    normexp.method:在method取“normexp”时,可以取"saddle","mle","rma"或者"rma75"。
    作者建议使用“mle”或者“saddle”,两个运行速度也差不多。如果数据中含有“形态背景估计(morphological background estimates)”,比如SPOT和GenePix图像处理软件,那么method = "subtract"也就较好的表现。
    3.2 组内归一化
    normalizeWithinArrays(object, layout, method="printtiploess", weights=object$weights, span=0.3, iterations=4, controlspots=NULL, df=5, robust="M", bc.method="subtract", offset=0)
    method:"none","median","loess","printtiploess","composite","control"和"robustspline"。初始值设定为“printtiploess”,但是对于Agilent芯片或者小样本芯片(每个“点样块”少于150个探针)。可用选用的loess或者robustspline。“loess归一化方法假设了有相当大的一部分探针没有发生差异化表达,而不是上调或者下调的基因数差不多或者差异化程度围绕着0波动。”(参考文献1)。需要注意,组内归一化只是对单张芯片的M值进行了规整(不影响A值),但是对与组间各个通道没有进行比较。
    weight:是图像处理软件对探针权重的标记,如果使用weight,那么在归一化过程中weight为0的探针不会影响其他探针,这并不意味着这些探针会被剔除,它们同样会被归一化,也会出现在归一化的结果中。如果不想使用,那么weight设为NULL。
    iterations:设定loess的循环数,循环数越多,结果越强健(robust)。
    3.3 组间归一化
    normalizeBetweenArrays(object, method=NULL, targets=NULL, cyclic.method="fast", ...)
    method:“none","scale","quantile","Aquantile","Gquantile",“Rquantile”,“Tquantile”,“cyclicloess”。“scale”对log-ratio数值进行归一化;“quantile”对密度(intensity)进行归一化;“Aquantile”对A数值进行归一化,不调正M数值;"Gquantile"和“Rquantile”则分别对Green和Red通道进行归一化,适用于Green或者Red为“参考标准(common reference)”的样品;“Tquantile”表示样品分开归一化或者Green和Red一先一后进行归一化。
    以下是一个cDNA芯片的密度图,分别为原始、去背景(method = "normexp", normexp.method="mle")、组内归一化(method = "robustspline")和组间归一化(method = “quantile”),使用plotDensities()绘制。
    limma包的使用技巧limma包的使用技巧

    limma包的使用技巧limma包的使用技巧

    4. 线性模型设计(linear model desigh)
    “线性模型”主要分析差异化表达基因,使用这种方法需要准备两个矩阵:“实验设计矩阵(design matrix)”和“对比矩阵(contrast matrix)”。“实验设计矩阵”主要用于每张芯片上用的RNA样品种类,“对比矩阵”主要用于规定实验组和对照组。
    “实验设计矩阵”:列表示RNA样品种类,对于oligo芯片或者使用common reference的cDNA芯片(简称第一类),列数与不同RNA样品数恰好相同;行数等于芯片数。对于不同染料对应不同样品的cDNA芯片(简称第二类),列数比RNA样品数恰好少一,行数等于芯片数(如要要衡量染料效应(dye effect),则列数与第一类相同)。对于第一类芯片,使用model.matrix()函数创建“实验设计矩阵”:比如oligo芯片model.matrix(~ 0+factor(c(1,1,1,2,2,3,3,3)))表示三类RNA样品;对于cDNA芯片,modelMatrix(targets,ref="EGFP")创建以“EGFP”为common reference的矩阵。对于第二类芯片,需要手动创建。
    “对比矩阵”通过函数makeContrasts()函数创建。
    使用步骤:创建“实验设计矩阵” --> lmFit() --> 创建“对比矩阵”(此步可选,如正交实验) -->contrasts.fit()(接上步,可选)--> eBayes() --> topTable()。文档详细列出了各种各样的例子。
     比较有用的例子:
    4.1 交换染料的技术重复
    ===========================================================================
    beta7pq$targets
          FileNames SubjectID  Cy3  Cy5 Date of Blood Draw Date of Scan
    1 6Hs.195.1.gpr         1 b7 - b7 +         2002.10.11   2003.07.25
    2   6Hs.168.gpr         3 b7 + b7 -         2003.01.16   2003.08.07
    3   6Hs.166.gpr         4 b7 + b7 -         2003.01.16   2003.08.07
    4 6Hs.187.1.gpr         6 b7 - b7 +         2002.09.16   2003.07.18
    5   6Hs.194.gpr         8 b7 - b7 +         2002.09.18   2003.07.25
    6 6Hs.243.1.gpr        11 b7 + b7 -         2003.01.13   2003.08.06
             Labels
    1 6Hs.195.1.gpr
    2   6Hs.168.gpr
    3   6Hs.166.gpr
    4 6Hs.187.1.gpr
    5   6Hs.194.gpr
    6 6Hs.243.1.gpr
    # 假定这六张芯片每两两做技术重复
    design <- c(1, -1, -1, 1, 1, -1)
    corfit <- duplicateCorrelation(beta7pq, design, ndups=1, block = c(1, 1, 2, 2, 3, 3), weight = NULL)
    fit <- lmFit(beta7pq, design, block = c(1, 1, 2, 2, 3, 3), cor = corfit$consensus)
    fit <- eBayes(fit)
    topTable(fit, adjust = "dfr")
               P.Value  adj.P.Val        B
    6647  4.870266e-07 0.01129122 5.341680
    11025 1.507433e-06 0.01747417 4.663784
    4910  9.223463e-06 0.04706320 3.434271
    11470 1.054914e-05 0.04706320 3.336518
    7832  1.182200e-05 0.04706320 3.252921
    22582 1.217992e-05 0.04706320 3.230931
    20812 1.939031e-05 0.05662524 2.882772
    1755  2.042581e-05 0.05662524 2.843204
    21405 2.743542e-05 0.05662524 2.616544
    11717 2.833754e-05 0.05662524 2.591458
    # 也可以采用以下方法
    design <- modelMatrix(beta7pq$targets, ref = "b7 -")
    design <- cbind(Dye = 1, desigh)
    fit <- lmFit(beta7pq, design)
    colnames(design)[2] <- "b7"
    cont.matrix <- makeContrasts(MUvsWT = b7/2, levels = design)
    fit2 <- contrasts.fit(fit, cont.matrix)
    fit2 <- eBayes(fit2)
    topTable(fit2, adjust = "fdr")
            P.Value  adj.P.Val        B
    6647  8.351055e-07 0.01936109 4.430939
    11025 3.555781e-06 0.02976803 3.700994
    11431 3.851970e-06 0.02976803 3.656632
    6211  6.431916e-06 0.03727939 3.362305
    11470 2.201258e-05 0.09075900 2.586146
    4910  2.495165e-05 0.09075900 2.501700
    1755  3.124884e-05 0.09075900 2.347645
    7832  3.131780e-05 0.09075900 2.346120
    11987 5.008082e-05 0.12764486 2.014907
    21859 5.505731e-05 0.12764486 1.946501
    #  两者有些不同,原因可能是第一种方法是专门为“对称”的cDNA芯片设计,而第二种是为非对称设计。
    ===========================================================================
    4.2 类型2的芯片
    参考文献8.4和8.5节。其中用到了model.matrix()这个函数,什么意思?
    4.3 正交实验
    参考文献8.6和8.7节。
    4.4 将两个通道分开分析
    参考文献8.8节

    5. 做图
    5.1 背景和前景
    imageplot(z, layout, low = NULL, high = NULL, ncolors = 123, zerocenter = NULL, zlim = NULL, mar=c(2,1,1,1), legend=TRUE, ...)
    z:可以定义R/Rb/G/Gb,也可以是导数;
    low/high:规定颜色,比如low = "white" , high = "red"
    下图是一张芯片的绿色前景图、红色背景图和log-ratio图(M值图):
    =======================================================
    imageplot(log2(beta7$G[, 1]), beta7$printer, low = "white", high = "green")
    imageplot(log2(beta7$Rb[, 1]), beta7$printer, low = "white", high = "red")
    imageplot(beta7q$M[, 1], beta7$printer)
    =======================================================
    limma包的使用技巧limma包的使用技巧limma包的使用技巧
    也可以使用imageplot3by2(RG, z="Gb", prefix=paste("image",z,sep="-"), path=NULL, zlim=NULL, common.lim=TRUE, ...),将图以3*2的格式六个一组保存成图片。
    5.2 M-A图
    使用plotMA(MA, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)画单个MA图。但这个函数有些问题,有时画不出。所以,完全可以自己来画,比如:
    ==========================================================
    # plot the MAplot befor normalization
    M <- log2((beta7$R[, 1]-beta7$Rb[, 1])/(beta7$G[, 1]-beta7$Gb[, 1]))
    A <- (log2((beta7$R[, 1]-beta7$Rb[, 1]))+log2((beta7$G[, 1]-beta7$Gb[, 1])))/2
    plot(A, M, cex = 0.5)
    # we can also plot the MAplot above using marray package
    beta7ma <- as(beta7, "marrayRaw")
    plot(maA(beta7ma)[,1], maM(beta7ma)[, 1])
    # plot the MAplot after normalization
    plot(beta7q$A[, 1], beta7q$M[, 1], cex = 0.5)
    ===========================================================
    limma包的使用技巧limma包的使用技巧

    使用plotPrintTipLoess(object,layout,array=1,span=0.4,...)画print-tip的MA图,下图分别是归一化前和后的print-tip的MA图。
    limma包的使用技巧

    limma包的使用技巧
    使用boxplot()画盒箱图,使用plotDensity()画密度曲线图。下图为原始数据、组内归一化(method = "normexp", normexp.method="mle")和组建归一化(method = "scale")的盒箱图。
    limma包的使用技巧

    limma包的使用技巧

    limma包的使用技巧
    最后,还可以使用 volcanoplot(fit, coef=1, highlight=0, names=fit$genes$ID, xlab="Log Fold Change", ylab="Log Odds", pch=16, cex=0.35, ...)画“火山图”(highlight标记前N个差异基因);
    使用vennDiagram(object, include="both", names, mar=rep(1,4), cex=1.5, lwd=1, circle.col, counts.col, show.include, ...)画“韦恩图”;
    使用plotMDS()绘制多纬标度图(multidimensional scaling plot, MDS)。




    参考文献:limma: Linear Models for Microarray Data User’s Guide

    转载于:https://www.cnblogs.com/huzs/p/3741979.html

    展开全文
  • 使用limma、Glimma和edgeR,RNA-seq数据分析 文章目录使用limma、Glimma和edgeR,RNA-seq数据分析六、差异表达分析1、创建设计矩阵和对比2、从表达计数数据中删除异方差3、拟合线性模型以进行比较4、检查DE基因数量5、...

    使用limma、Glimma和edgeR,RNA-seq数据分析


    本篇继续上一篇数据预处理之后做差异表达分析。

    六、差异表达分析

    1、创建设计矩阵和对比

    在此研究中,我们想知道哪些基因在我们研究的三组细胞之间以不同水平表达。我们的分析中所用到的线性模型假设数据是正态分布的。首先,我们要创建一个包含细胞类型以及测序泳道(批次)信息的设计矩阵。

    #创建一个包含细胞类型和测序泳道信息的设计矩阵
    design<-model.matrix(~0+group+lane)
    colnames(design)<-gsub("group","",colnames(design))
    

    design:
    在这里插入图片描述
    对于一个给定的实验,通常有多种等价的方法都能用来创建合适的设计矩阵。 比如说,~0+group+lane去除了第一个因子group的截距,但第二个因子lane的截距被保留。 此外也可以使用~group+lane,来自group和lane的截距均被保留。 理解如何解释模型中估计的系数是创建合适的设计矩阵的关键。 我们在此分析中选取第一种模型,因为在没有group的截距的情况下能更直截了当地设定模型中的对比。用于细胞群之间成对比较的对比可以在limma中用makeContrasts函数设定。

    contr.matrix<-makeContrasts(
      BasalvsLP=Basal-LP,
      BasalvcML=Basal-ML,
      LPvsML=LP-ML,
      levels = colnames(design))
    

    在这里插入图片描述
    limma线性模型方法的核心优势之一便是其适应任意实验复杂程度的能力。简单的实验设计,比如此流程中关于细胞类型和批次的实验设计,直到更复杂的析因设计和含有交互作用项的模型,都能够被较简单地处理。当实验或技术效应可被随机效应模型(random effect model)模拟时,可以使用limma中的duplicateCorrelation函数来估计交互作用,这需要在此函数以及lmFit的线性建模步骤均指定一个block参数。

    2、从表达计数数据中删除异方差

    对于RNA-seq计数数据而言,原始计数或log-CPM值的方差并不独立于均值(Law et al. 2014)。有些差异表达分析方法使用负二项分布模型,假设均值与方差间具有二次的关系。而在limma中,假设log-CPM值符合正态分布,因此我们在对RNA-seq的log-CPM值进行线性建模时,需要使用voom函数计算每个基因的权重从而调整均值与方差的关系,否则分析得到的结果可能是错误的。

    当操作DGEList对象时,voom从x中自动提取文库大小和归一化因子,以此将原始计数转换为log-CPM值。在voom中,对于log-CPM值额外的归一化可以通过设定normalize.method参数来进行。

    下图左侧展示了这个数据集log-CPM值的均值-方差关系。通常而言,方差是测序实验操作中的技术差异和来自不同细胞类群的重复样本之间的生物学差异的结合,而voom图会显示出一个在均值与方差之间递减的趋势。 生物学差异高的实验通常会有更平坦的趋势,其方差值在高表达处稳定。生物学差异低的实验更倾向于急剧下降的趋势。

    不仅如此,voom图也提供了对于上游所进行的过滤水平的可视化检测。如果对于低表达基因的过滤不够充分,在图上表达低的一端,受到非常低的表达计数的影响,会出现方差的下降。如果观察到了这种情况,应当回到最初的过滤步骤并提高用于该数据集的表达阈值。

    当先前绘制的MDS图中发现组内重复样本的聚集与分离程度出现明显的差异时,可以用voomWithQualityWeights函数(Liu et al. 2015)来代替voom,在计算基因权重值以外还能计算每个样本的权重值。

    par(mfrow=c(1,2))
    v<-voom(x,design,plot=TRUE)
    vfit<-lmFit(v,design)
    vfit<-contrasts.fit(vfit,contrasts = contr.matrix)
    efit<-eBayes(vfit)
    plotSA(efit,main="Final model:Mean-variance trend")
    

    在这里插入图片描述

    • 图中绘制了每个基因的均值(x轴)和方差(y轴),显示了在该数据上使用voom前它们之间的相关性(左),以及当运用voom的权重后这种趋势是如何消除的(右)。左侧的图是使用voom函数绘制的,它为log-CPM转换后的数据拟合线性模型并提取残差方差。然后,对方差取四次方根(或对标准差取平方根),并相对每个基因的平均表达作图。均值通过平均计数加上2再进行log2转换计算得到。右侧的图使用plotSA绘制了log2残差标准差与log-CPM均值的关系。在这两幅图中,每个黑点表示一个基因。左侧图中,红色曲线展示了用于计算voom权重的估算所得的均值-方差趋势。右侧图中,由经验贝叶斯算法得到的平均log2残差标准差由水平蓝线标出。

    3、拟合线性模型以进行比较

    limma的线性建模使用lmFit和contrasts.fit函数进行,它们原先是为微阵列而设计的。这些函数不仅可以用于微阵列数据,也可以用于使用voom计算了基因权重后的RNA-seq数据。每个基因的表达值都会单独拟合一个模型。然后通过借用全体基因的信息来进行经验贝叶斯调整(empirical Bayes moderation),这样可以更精确地估算各基因的差异性(Smyth 2004)。图4为此模型的残差关于平均表达值的图。从图中可以看出,方差不再与表达水平均值相关。

    4、检查DE基因数量

    为快速查看表达差异水平,显著上调或下调的基因可以汇总到一个表格中。显著性的判断使用默认的校正p值阈值,即5%。在basal与LP的表达水平之间的比较中,发现了4648个在basal中相较于LP下调的基因,和4863个在basal中相较于LP上调的基因,即共9511个差异表达基因。在basal和ML之间发现了一共9598个差异表达基因(4927个下调基因和4671个上调基因),而在LP和ML中发现了一共5652个差异表达基因(3135个下调基因和2517个上调基因)。在basal类群所参与的比较中皆找到了大量的差异表达基因,这与我们在MDS图中观察到的结果相吻合。

    summary(decideTests(efit))
    

    在这里插入图片描述
    在某些时候,不仅仅需要用到校正p值阈值,还需要差异倍数的对数(log-FCs)也高于某个最小值来更为严格地定义显著性。treat方法(McCarthy and Smyth 2009)可以按照对最小log-FC值的要求,使用经过经验贝叶斯调整的t统计值计算p值。当我们要求差异表达基因的log-FC显著大于1(等同于不同细胞类群之间表达量差两倍)并对此进行检验时,差异表达基因的数量得到了下降,basal与LP相比只有3648个差异表达基因,basal与ML相比只有3834个差异表达基因,LP与ML相比只有414个差异表达基因。

    tfit<-treat(vfit,lfc = 1)
    dt<-decideTests(tfit)
    summary(dt)
    

    在这里插入图片描述
    在多个对比中皆差异表达的基因可以从decideTests的结果中提取,其中的0代表不差异表达的基因,1代表上调的基因,-1代表下调的基因。共有2784个基因在basal和LP以及basal和ML的比较中都差异表达,其中的20个于下方列出。write.fit函数可用于将三个比较的结果一起提取并写入一个输出文件。

    de.common<-which(dt[,1]!=0 & dt[,2]!=0)
    length(de.common)
    head(tfit$genes$SYMBOL[de.common],n=20)
    vennDiagram(dt[,1:2],circle.col =c("turquoise","salmon"))
    

    在这里插入图片描述

    • 韦恩图展示了仅basal和LP(左)、仅basal和ML(右)的对比的DE基因数量,还有两种对比中共同的DE基因数量(中)。在任何对比中均不差异表达的基因数量标于右下。
    write.fit(tfit,dt,file = "result.txt")#保存处理结果
    

    5、检查单个DE基因

    使用topTreat函数可以列举出使用treat得到的结果中靠前的DE基因(对于eBayes的结果则应使用topTable函数)。默认情况下,topTreat将基因按照经过多重假设检验校正的p值从小到大排列,并为每个基因给出相关的基因信息、log-FC、平均log-CPM、校正t值、原始及校正p值。列出前多少个基因的数量可由用户设定,如果设为n=Inf则会包括所有的基因。基因Cldn7和Rasef在basal与LP和basal于ML的比较中都位于DE基因的前几名。

    basal.vs.lp<-topTreat(tfit,coef = 1,n=Inf)
    basal.vs.ml<-topTreat(tfit,coef = 2,n=Inf)
    head(basal.vs.lp)
    head(basal.vs.ml)
    

    在这里插入图片描述

    6、差异表达结果可视化

    为可视化地总结所有基因的结果,可使用plotMD函数绘制均值-差异关系(MD)图,其中展示了线性模型拟合所得到的每个基因log-FC与log-CPM平均值间的关系,而差异表达的基因会被重点标出。
    limma:

    #差异表达结果可视化
    plotMD(tfit,column = 1,status = dt[,1],main=colnames(tfit)[1],xlim=c(-8,13))
    

    在这里插入图片描述
    Glimma的glMDPlot函数提供了交互式的均值-差异图,拓展了这种图表的功能。此函数的输出为一个html页面,左侧面板为结果的总结性图表(与plotMD的输出类似),右侧面板展示了各个样本的log-CPM值,而下方面板为结果的汇总表。这使得用户可以使用提供的注释中的信息(比如基因符号)搜索特定基因,而这在R统计图中是做不到的。
    Glimma:

    glMDPlot(tfit,coef = 1,status = dt,main = colnames(tfit)[1],side.main = "ENTREZID",counts = lcpm,groups = group,launch = TRUE)
    

    在这里插入图片描述
    使用Glimma生成的均值-差异图。左侧面板显示了总结性数据(log-FC与log-CPM值的关系),其中选中的基因在每个样本中的数值显示于右侧面板。下方为结果的表格,其搜索框使用户得以使用可行的注释信息查找某个特定基因,如基因符号Clu。

    • Glimma提供的交互性使得单个图形窗口内能够呈现出额外的信息。 Glimma是以R和Javascript实现的,使用R代码生成数据,并在之后使用Javascript库D3(https://d3js.org)转换为图形,使用Bootstrap库处理界面并生成互动性可搜索的表格的数据表。这使得图表可以在任意当代浏览器中直接查看而无需后段服务器时刻运行,对于将其作为关联文件附加在Rmarkdown分析报告中而言非常方便。

    前文所展示的图表中,一些仅展示了在任意一个条件下表达的所有基因而缺少单个基因的具体信息(比如不同对比中共同差异表达基因的韦恩图或均值-差异图),而另一些仅展示单个基因(交互式均值-差异图右边面板中所展示的log-CPM值)。而热图使用户得以查看某一小组基因的表达情况,既便于查看单个组或样本的表达,又不至于在关注于单个基因时失去对于整体的注意,也不会造成由于对上千个基因取平均值而导致的分辨率丢失。

    使用gplots包的heatmap.2函数,我们为basal与LP的对照中前100个差异表达基因(按校正p值排序)绘制了一幅热图。此热图中正确地将样本按照细胞类型聚类,并重新排序了基因,表达相似的基因形成了块状。从热图中,我们发现basal与LP之间的前100个差异表达基因在ML和LP样本中的表达非常相似。

    #gplots绘制热图
    library(gplots)
    basal.vs.lp.topgenes<-basal.vs.lp$ENTREZID[1:100]#取100个即可,没有对某一特征进行排序
    i<-which(v$genes$ENTREZID %in% basal.vs.lp.topgenes)
    mycol<-colorpanel(1000,"blue","white","red")
    heatmap.2(lcpm[i,],scale = "row",labRow = v$genes$SYMBOL[i],labCol = group,
              col = mycol,trace = "none",density.info = "none",margins = c(8,6),
              lhei = c(2,10),dendrogram = "column")
    

    在这里插入图片描述

    • 在basal和LP的对比中前100个DE基因log-CPM值的热图。经过缩放调整后,每个基因(每行)的表达均值为0,并且标准差为1。给定基因相对高表达的样本被标记为红色,相对低表达的样本被标记为蓝色。浅色和白色代表中等表达水平的基因。样本和基因已通过分层聚类的方法重新排序。图中显示有样本聚类的树状图。

    7、使用camera的基因检验

    在此次分析的最后,我们要进行一些基因集检验。为此,我们将camera方法(Wu and Smyth 2012)应用于Broad Institute的MSigDB c2中的(Subramanian et al. 2005)中适应小鼠的c2基因表达特征,这可从http://bioinf.wehi.edu.au/software/MSigDB/以RData对象格式获取。 此外,对于人类和小鼠,来自MSigDB的其他有用的基因集也可从此网站获取,比如标志(hallmark)基因集。C2基因集的内容收集自在线数据库、出版物以及该领域专家,而标志基因集的内容来自MSigDB,从而获得具有明确定义的生物状态或过程。

    load(system.file("extdata", "mouse_c2_v5p1.rda", package = "RNAseq123"))
    idx <- ids2indices(Mm.c2,id=rownames(v))
    #Basal 与LP做对比
    cam.BasalvsLP <- camera(v,idx,design,contrast=contr.matrix[,1])
    head(cam.BasalvsLP,5)
    #Basal 与ML做对比
    cam.BasalvsML <- camera(v,idx,design,contrast=contr.matrix[,2])
    head(cam.BasalvsML,5)
    #LP与ML做对比
    cam.LPvsML <- camera(v,idx,design,contrast=contr.matrix[,3])
    head(cam.LPvsML,5)
    

    在这里插入图片描述
    camera函数通过比较假设检验来评估一个给定基因集中的基因是否相对于不在集内的基因而言在差异表达基因的排序中更靠前。 它使用limma的线性模型框架,并同时采用设计矩阵和对比矩阵(如果有的话),且在检验的过程中会运用到来自voom的权重值。 在通过基因间相关性(默认设定为0.01,但也可通过数据估计)和基因集的规模得到方差膨胀因子(variance inflation factor),并使用它调整基因集检验统计值的方差后,将会返回根据多重假设检验进行了校正的p值。

    绘制LP与ML的条码图

    #绘制条码图
    barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP, 
                index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")
    

    process=image/watermark,type_ZHJvaWRzYW5zZmFsbGJhY2s,shadow_50,text_Q1NETiBA5p-a5a2Q5ZGz55qE576K,size_20,color_FFFFFF,t_70,g_se,x_16)
    今日告一段落啦,明天见,加油呀~

    展开全文
  • limma

    千次阅读 2019-07-10 14:34:00
    limma:Linear Models for Microarray and RNA-Seq Data http://www.bioconductor.org/packages/release/bioc/html/limma.html 安装 if (!requireNamespace("BiocManager", quietly = TRUE)) install.package...

    limma:Linear Models for Microarray and RNA-Seq Data

    http://www.bioconductor.org/packages/release/bioc/html/limma.html

    安装

    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    
    BiocManager::install("limma")

    使用

    library("limma")

    usersguide:http://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf

    limma中的函数查询:http://www.bioconductor.org/packages/release/bioc/manuals/limma/man/limma.pdf

    
    

    转载于:https://www.cnblogs.com/0820LL/p/11163776.html

    展开全文
  • 生信入门(二)——使用limma、Glimma和edgeR,RNA-seq数据分析 文章目录生信入门(二)——使用limma、Glimma和edgeR,RNA-seq数据分析一、简介 一、简介 简单且高效地分析RNA 测序数据的能力是生信的核心优势之一。...
  • 在基因差异分析的过程中常见的几大差异化常用的R包 主要是edgeR Deseq2 和limma(其中limma 包主要内置于edgeR) 那么在分析的过程中呢,每个包有他们各自的一些数据处理方式。 在接受的数据格式上和样本数目上呢也...
  • 最近工作关系,需要重现一个文章的基因芯片数据分析,查找差异基因,花了一天时间跑了limma流程,供大家参考。 论文名字为 Identification of inflammatory mediators in patients with Crohn’s disease ...
  • TCGA_RNA-seq_limma分析

    2021-05-15 11:01:44
    limma分析的大致流程 导入read count, 保存为专门的对象用于后续分析 原始数据过滤,根据标准化read count 或者 raw count 作为筛选标准 raw read count 标准化 通过各种算法(如经验贝叶斯,EM)预测dispersion...
  • limma包和WGCNA包进行RNA-seq数据分析 #数据提取# GE('TCGA-COAD.htseq_counts.tsv',header=T,sep='\t',stringsAsFactors = F) #60488*513 512个样本,其中对照组41个 # group_data(colnames(GE)[-1]) group('tumor...
  • 好吧具体它怎么工作的咱也看不懂(voom参考文献来源)……不过搞懂它的分析流程,以及结果怎么解读,还是可以的 limma_voom (dgelist_norm, design, plot = TRUE) limma_voom fit (limma_voom, design) #拟合 fit ...
  • 1.两组间阐述较清楚 2.大致相似,多组间差异表达看这里 3.生信技能树,含topTreat函数+韦恩图 4.生信星球,差异表达分子分母弄错了但解释较清楚
  • limma包进行多组差异表达分析

    千次阅读 多人点赞 2020-03-15 00:20:16
    写在前面:最近在使用limma包进行差异表达分析,参考了网上许多教程都觉得说的云里雾里,很不清楚。经过我自己一段时间非常痛苦的钻研,弄明白了,解决了我的实际需求。于是决定将我的分析经验写下来,分享给需要的...
  • 差异分析流程(二)limma

    千次阅读 2019-12-14 11:16:34
    差异分析流程(二)limma1. 制作三个矩阵2. 建模前归一化3. 建模及结果4.检查 参考链接: http://www.freesion.com/article/752576024/ ......
  • # 清空环境 rm(list = ls()) options(stringsAsFactors = F) # 读取counts矩阵 counts ('./UUO-Sham_count.csv', ...ctrl', n=Inf) DEG_limma_voom = na.omit(tempOutput) write.csv(tempOutput, './DEGs_voom.csv')
  • 简单使用limma做差异分析

    万次阅读 2019-04-10 21:42:31
    首先需要说明的是,limma是一个非常全面的用于分析芯片以及RNA-Seq的差异分析,按照其文章所说: limma is an R/Bioconductor software package that provides an integrated solution for...
  • write.csv(DE, "limma_genes_test.csv") selected = rownames(DE[DE$P.Value ,]) # 默认是holm方法控制FWER esetSel = y[selected, ] heatmap(esetSel@.Data[[1]]) #p[DE$P.Value ,]$P.Value #esetSela [] #p....
  • TCGA的差异分析(limma和edgeR)

    千次阅读 2020-04-19 16:41:37
    这篇文章简要介绍一下limma 包和edgeR包找TCGA差异基因的代码,主要是备忘
  • limma包进行差异分析

    千次阅读 2019-09-22 17:23:44
    limma包进行差异分析 载入limma包,利用edgeR的标准化方法 library(edgeR) library(limma) 读取两组对比的数据 a<-read.table("E3和E4.txt") E3有81个样本,E4有190个样本,所以创建一个dataframe来统计样本的...
  • 典型医学设计实验GEO数据分析 (step-by-step) - 数据获取到标准化介绍了实验的设计、数据获取、数据标准化和注释,下面是如何利用Limma和线性模型鉴定差异基因,并进行GO富集分析。 线性模型 为了分析发炎和未发炎...
  • GSVA基因集变异分析结合limma包分析差异基因集 基因集变异分析即GSVA(Gene set variation analysis),是一种非参数的无监督分析方法,主要用来评估芯片和转录组的基因集富集结果,适合单细胞转录组测序基因集的差异...
  • 1. 加载R包和输入数据 ...library("limma") library("edgeR") expr = read.csv("mRNA_exprSet.csv",sep = ',',header=T) head(expr) TCGA-mRNA数据链接 链接:https://pan.baidu.com/s/1b6l4qg4...
  • 使用limma包进行差异基因分析时,做最多的是两分类的,例如control组和disease组,但也会碰到按照序列进行的分组。这时,如果逐一使用两两比较求差异基因则略显复杂。其实开发limma包的大神们已经替我们考虑到。我...
  • limma:RNA-Seq Data

    2020-01-06 12:41:51
    生成计数矩阵 标准化和过滤 第一步先创建DGEList对象。 dge = DGEList(counts = counts) 移除0和低计数值得行。 keep = filterByExpr(dge,design) dge = dge[keep,,keep.lib....差异表达:limma-trend ...
  • 我用beadarry和limma包对GSE49454数据集进行差异基因表达分析时,最后卡在这个group name上面了,尝试了好多方法都不行 代码:library(GEOquery) library(beadarray) library(illuminaHumanv4.db) #下载数据 ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 707
精华内容 282
关键字:

limma