edgeR基因表达差异分析

edgeR基因表达差异分析

文章目录

• edgeR进行基因表达差异分析

  • 官方文档概述
    • 读取read值并初始化数据结构
      • 创建DGEList对象并设置相应的分组信息
      • 施行数据清洗步骤以去除低表达基因样本
      • 进行CPM尺度转换处理以标准化数据分布特征
    • 设计手动筛选流程以辅助数据预处理工作流程
    • 配置自动筛选机制以实现高效的数据质量监控功能

• 标准化

  • • 测序深度

• 有效数据集规模

• GC比例

• 基因长度

• MDS图形展示用于样本的无监督聚类分析。

  复制代码

  * 负二项式模型
  * * 计算生物变异系数

  • 计算差异基因

    • 广义线性模型（Glm）
    • • 计算离散度

  • 计算DE基因

    • 如果没有重复样本？
    • 输出结果
    • 查看统计

参考文献：提供详细说明的示例。
官方主页：访问边距分析工具包的页面。
RNAseq123资源：查看关于使用Limma工具包分析RNASeq数据的详细指南。

官方文档总结

注意⚠️：

edgeR被设计用于处理实际读取到的计数数据。我们强烈推荐避免直接导入预期的转录本数量到边距R中

读取read数

其中一个方法是使用read.table()或read.delim()；当文件数量较多时，则采用readDGE(files, columns=)来一次性读取多个文件。注意，请确保readDGE()函数被配置为指定两列：其中一列为计数数据，另一列为基因标识符。

    x <- readDGE(files, columns=c(1,3))

DGEList对象、构建分组

DGEList是一种能够整合多种数据与统计分析功能的数据容器。该列表的基本构成要素包括两个核心组件：计数矩阵（\texttt{counts}）与样本信息表（\texttt{samples}），其中样本信息表进一步细分为分组标识符（group）与各实验样本的大小参数（lib.size）。计数矩阵（\texttt{counts}）用于存储基因表达水平的数据矩阵，在此框架下可进行多样化的统计分析操作。样本信息表（\texttt{samples}）则用于记录每个样本的身份标记以及对应的测序库大小，并在此基础上完成标准化处理及差异分析计算过程。

    # 构建一个含有组别标记的DGEList
    y<- DGEList(counts=x)
    group <-  c(1,1,2,2)
    y<- DGEList(counts=x,group=group)

在DGEList中进行分类分组是必要步骤。此外还可以添加额外的信息进去。例如lane道、基因长度以及基因注释等信息。

过滤，删除低表达基因

然而由于这些低表达基因既没有实际应用价值又会对实验结果造成干扰，在进行 downstream分析时应当筛选出所有样本中序列片段数量不足的基因并予以剔除

1. 在生物学术语中, 低表达指标在生物学术语中没有实际意义。
2. 通过去除那些在生物学术语中无显著意义的数据点, 在统计学分析中可以更准确地评估均值与方差之间的关系。
3. 通过筛选掉那些在生物学术语中无显著意义的数据点, 在后续差异基因分析中将计算开销降低。

该R包中的filterByExpr函数提供了一种自动筛选基因的方法，在此过程中能够尽可能保留那些具有足够高表达计数的基因。该函数基于最小组内样本数量作为基准，并筛选出在至少达到这个基准数量的样本中具有10个或更多序列片段计数的基因。具体而言，在所有可能的情况下，默认要求每个被保留下来的基因必须在其所属组内的所有样本中都满足最低计数值的要求（此处最低计数值等于该组内的最小样本数量），换算为cpm即cut.off.cpm=10/<sum_col_counts>。

CPM标度转换

常用的标度转换有CPM（counts per million）、log-CPM、FPKM、RPKM.
CPM 是将counts转变为CPM指数counts per million，消除测序深度影响。
log-CPM 是将CPMlog2化。cpm函数会在进行log2转换前给CPM值加上一个弥补值。默认的弥补值是2/L其中2是“预先计数”，而L是样本总序列数（以百万计）的平均值，所以log-CPM值是根据CPM值通过log2(CPM + 2/L)计算得到的。
对于一个基因，CPM值为1相当于在序列数量约2千万的样品中，有20个计数，或者在序列数量约7.6千万有76个计数。

    cpm<- edgeR::cpm(x)
    lcpm <- edgeR::cpm(x, log=TRUE)

手动过滤

经验设定为cpm=1作为分界线。然而这并非最优选择 因为当测序深度提升时 例如达到2千万时 cpm=1意味着counts=20 此时阈值可能偏高。而当测序深度较低时 例如仅达到两百万时 cpm=1则对应counts=1 此时阈值可能偏低。因此建议采用自动筛选功能 并根据cut.off.cpm设定为10与<sum_col_counts>的比例来进行筛选

    cpm(y)
    keep_cpm <- rowSums(cpm(y)>1) >= 3      #此例设置最小组内样本数为3
    x<- x[keep_cpm,,keep.lib.sizes=FALSE]

自动过滤

使用edgeR::filterByExpr()自动过滤

    keep.exprs<- edgeR::filterByExpr(x,group=group)
    x<- x[keep.exprs,,keep.lib.sizes=FALSE]

归一化

归一化并非绝对必要但建议仍需执行归一化处理。
在存在重复样本的情况下 由于生物学背景合理的外部因素可能会影响单个样品的表现 因此我们通常会观察到第一批制备得到的产品在全局上表现出更高的表达水平与第二批结果相比 假设所有样品表达值范围与分布应保持一致 必须对实验数据进行归一化处理以确保每个样本的表现一致性 可以使用edgeR::calcNormFactors()函数来采用TMM方法执行归一化处理；归一化过程中所计算出的归一化系数将被自动存储在DGElist对象对应的samples属性下的norm.factors字段中。

测序深度

差异性文库大小反映了不同的测序深度。
这属于基本建模流程的一部分，并自然地纳入到倍数变化或p值计算的过程中。
这一现象持续存在，并无需任何用户主动干预。

有效库大小

在特定情境下

    y <- calcNormFactors(y)

请注意：归一化不会直接修改counts数值，并非会直接进行数值调整；相反地，在系统内部会预先计算并存储归一化系数于变量x$samples$norm.factors中以供后续使用。

GC含量

基因长度

由于每个基因间的GC含量及基因长度保持恒定，在这种情况下这些数值被视为相对指标；从而对差异分析的影响较小。

MDS图形展示 样本无监督聚类

这种图表采用了无监督聚类方法来直观地呈现样本之间的相似性和差异性，并有助于我们在开展正式分析前对可识别差异表达基因的数量形成大致了解。

    # 图形展示
    library(RColorBrewer)
    lcpm <- cpm(y, log=TRUE)
    col.group <- group
    levels(col.group) <-  brewer.pal(nlevels(col.group), "Set1")
    col.group <- as.character(col.group)
    plotMDS(lcpm,labels=group,col=col.group)
    
    #或者直接使用 
    plotMDS(y)

负二项式模型

计算生物变异系数

在edgeR软件包中，默认情况下用于估计离散度的主要有两种模型，在生物信息学领域内被广泛采用的就是qCML方法；然而该方法仅适用于涉及单一因素的数据集分析；例如当研究癌细胞与健康细胞之间的基因表达差异时无需考虑其他潜在影响因素如年龄或性别等；随后该方法将采用类似于Fisher精确检验的精确检验程序来进行差异表达分析；如果采用qCML方法进行离散度估计，则无需预先定义实验矩阵的设计矩阵！这是因为该方法仅考虑单一因素的影响情况；对于更为复杂的实验设计方案，则建议采用广义线性模型（GLMs）框架进行建模与分析；QCMLEstimation Methodology：主要应用于具有单一研究因素的数据集结构；通过estimateDisp()函数可一次性完成公共离散度与tagwise离散度的估算

    y < -estimateDisp（y）

或者先计算common离散度，再计算tagwise离散度

    y <- estimateCommonDisp(y)
    y <- estimateTagwiseDisp(y)

计算差异基因

    et<- exactTest(y)
    topTags(et)

广义线性模型（Glm）

对于更为复杂的实验设计（包含多个变量），广义线性模型可用于评估其性能水平

计算离散度

通过下面来估计common离散度、trended离散度、tagwise离散度。

    y<- estimateDisp(y,design)

或者分开依次进行

    y <- estimateGLMCommonDisp(y, design)
    y <- estimateGLMTrendedDisp(y, design)  #估计trended离散度
    y <- estimateGLMTagwiseDisp(y, design)

关于design的构建(?)：

    #group是区别分组的factor对应不同样本，lane是不同样本的的lane道
    design <- model.matrix(~0+group+lane) 
    colnames(design) <- gsub("group", "", colnames(design))

计算DE基因

该包包含了多种用于差异表达分析的函数。
为了准确估计离散度并进行差异性分析,qCML方法需配合exact test使用,其对应的函数是exactTest。
其他四个与GLM模型相关的函数包括(glmFit与glmlrt)以及它们各自的作用:glmmtest用于似然比检验/glmmftest则应用于拟极大似然F检验。
该包提供了两个核心测试方法: glmQLFit()和 glmQLFTest().推荐优先选择glq_fit方法,因为它能够有效评估每个基因离散度下的不确定性.在样本量较少的情况下,该方法能提供更为稳健和精确的错误率控制。

    # 例子
    group <- factor(c(1,1,2,2,3,3))
    design <- model.matrix(~group)
    fit <- glmQLFit(y, design)

如果没有重复样本？

首先确定一个良好的BCV，人类数据一般设定为0.4
一个例子：

    bcv <- 0.2
    counts <- matrix( rnbinom(40,size=1/bcv^2,mu=10), 20,2)
    y <- DGEList(counts=counts, group=1:2)
    et <- exactTest(y, dispersion=bcv^2)

如果是人类数据：

    y_bcv <- y
    bcv <- 0.4
    et <- exactTest(y_bcv,dispersion = bcv ^2)

注：
如果差异基因过少，可以尝试下调BCV

输出结果

    result = topTags(et, n = nrow(et$table))$table

查看统计

    summary(de <- decideTestsDGE(et))

我希望创建并维护一个高质量的生信与统计相关的微信交流群。如需参与讨论的朋友，请添加我的微信： veryqun 。我会将您邀请进群，并且如果有任何问题或疑问也可以通过微信与我联系。

因为我的笔记主要分布在博客、简书平台以及知乎专栏等分散的位置,日常维护较为繁琐.经过深思熟虑后开通了一个 微信公众号账号,这将显著提升内容的发布效率,更方便地分享更多关于生物信息及统计分析的内容.如需进一步了解或合作意向,请随时与我联系 ^ . ^ 。

对我唯一的精神鼓励可能就是下方的点赞了吧 ^ ^