edgeR基因表达差异分析
文章目录 edgeR基因表达差异分析官方文档总结读取read数DGEList对象、构建分组过滤，删除低表达基因CPM标度转换手动过滤自动过滤 归一化测序深度有效库大小GC含量基因长度MDS图形展示 样本无监督聚类 负二项式模型计算生物变异系数计算差异基因 广义线性模型（Glm）计算离散度计算DE基因 如果没有重复样本？输出结果查看统计
参考：
一个比较详细的例子：http://www.iwhgao.com/edger%E7%AE%80%E5%8D%95%E4%BD%BF%E7%94%A8%E6%89%8B%E5%86%8C/
bioconductor主页：http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
RNAseq123:http://www.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html
官方文档总结
注意⚠️：
edgeR设计用于实际读取计数。我们不建议将预测的转录本丰度输入到edgeR而不是实际计数中。 读取read数
可以使用read.table)，read.delim),如果文件比较多，readDGEfiles, columns=)来一次性读取多个文件。注意，readDGE)需要指定两列，一列用于计数，一列用于基因标识符
x <- readDGEfiles, columns=c1,3)) DGEList对象、构建分组
DGEList是一个可以包含多种内容和统计的列表。DGEList至少需要的元素：counts、samples（包含group分组信息和lib.size文库大小），counts用来存放表达矩阵，samples用来标记样本信息和库的大小，group声明组别 。
# 构建一个含有组别标记的DGEListy<- DGEListcounts=x)group <- c1,1,2,2)y<- DGEListcounts=x,group=group)
DGEList中分组是必要的，也可以添加其他信息进去，例如lane道，基因长度，基因注释等信息
一个注释过的DGEList
过滤，删除低表达基因
低表达基因不仅用不到，而且会干扰结果，所以要去除在任何样本中都没有足够多的序列片段的基因应该从下游分析中过滤掉，因为：
低表达没有生物学意义去除低表达数据可以对数据中均值-方差关系有更精确的估计减少了观察差异表达下游分析中的运算量
edgeR包中的filterByExpr函数提供了自动过滤基因的方法，可保留尽可能多的有足够表达计数的基因。此函数默认选取最小的组内的样本数量为最小样本数，保留至少在这个数量的样本中有10个或更多序列片段计数的基因。
过滤标准是，以最小对组内样本数为标准，（此例最小组内样本为3），如果有基因在所有样本中表达数（count）小于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值通过log2CPM + 2/L)计算得到的。
对于一个基因，CPM值为1相当于在序列数量约2千万的样品中，有20个计数，或者在序列数量约7.6千万有76个计数。
cpm<- edgeR::cpmx)lcpm <- edgeR::cpmx, log=TRUE) 手动过滤
经验设置为cpm=1位为cutoff点。但是，这并不是最精准。因为随着测序深度增加，例如20million（2 千万），cpm=1 意味着 counts=20。阈值可能会有点高。测序深度低的话，例如2million（2百万），cpm=1 意味着counts=1。阈值可能会太低。此时可以使用自动过滤，或者根据cut.off.cpm=10/<sum_col_counts>来计算。
cpmy)keep_cpm <- rowSumscpmy)>1) >= 3 #此例设置最小组内样本数为3x<- x[keep_cpm,,keep.lib.sizes=FALSE] 自动过滤
使用edgeR::filterByExpr)自动过滤
keep.exprs<- edgeR::filterByExprx,group=group)x<- x[keep.exprs,,keep.lib.sizes=FALSE] 归一化
归一化不是绝对必要的，但是推荐进行归一化。
有重复的样本中，应该不具备生物学意义的外部因素会影响单个样品的表达，例如中第一批制备的样品会总体上表达高于第二批制备的样品，假设所有样品表达值的范围和分布都应当相似，需要进行归一化来确保整个实验中每个样本的表达分布都相似。可以使用edgeR::calcNormFctors)函数来用TMM方法进行归一化,通过计算归一化系数来决定文库大小的缩放系数当进行归一化时，如果对象是DGElist，那么这些归一化系数会被自动存在x$samples$norm.factors
测序深度
不同文库大小代表不同测序深度。这是基本建模过程的一部分，并自动进入倍数变化或p值计算。它始终存在，不需要任何用户干预。
有效库大小
在某些情况下，如果一小部分高表达基因消耗了一个特定样本文库大小的很大一部分，这将导致该样本剩余基因的采样不足。除非对这种效应进行调整，否则在该样本中其余的基因可能会错误地表现为下调。（因为高表达基因比较多，消耗了大量的测序资源，导致非高表达基因会相对偏低。）
使用函数edgeR::calcNormFactors)，默认使用TMM方法进行归一化，归一化后，会给样品分配缩放系数。
将原始库大小与缩放因子的乘积称为有效库大小。有效的库大小会在所有下游分析中替换原始库大小。
y <- calcNormFactorsy)
注意⚠️：
归一化并不会直接在counts数值上修改，而是归一化系数会被自动存在x$samples$norm.factors。
GC含量 基因长度
因为GC含量和基因长度在每个基因间不会改变，因此这些都是相对的，因此对差异分析影响很小。
MDS图形展示 样本无监督聚类
这种图表使用无监督聚类方法展示出了样品间的相似性和不相似性，能让我们在进行正式的检验之前对于能检测到多少差异表达基因有个大致概念。
# 图形展示libraryRColorBrewer)lcpm <- cpmy, log=TRUE)col.group <- grouplevelscol.group) <- brewer.palnlevelscol.group), “Set1”)col.group <- as.charactercol.group)plotMDSlcpm,labels=group,col=col.group)#或者直接使用 plotMDSy)
负二项式模型 计算生物变异系数
edgeR中对于离散度的估计主要有两种模型，其中比较经典的是qCML。但qCML方法仅适用于单因素的数据集。例如，比较癌和正常样本间的基因表达差异，而不考虑其他因素（如年龄、性别等）的影响。qCML后续会使用与fisher精确检验比较类似的 exact test 进行差异表达分析。
如果使用qCML估计离散度，则不需要预先设定实验矩阵！因为就只看一个因素，不必大费周折)
而对于更复杂的实验设计，则推荐使用GLMs模型，此时就需要提供实验矩阵design。
qCML方法：仅适用于具有单因素设计的数据集，但是有着更可靠的性能，尤其在小样本下有很好的表现，比较切合NGS测序。
通过estimateDisp)来估计公共离散度和 tagwise离散度（一个命令来运行）
y < -estimateDisp（y）
或者先计算common离散度，再计算tagwise离散度
y <- estimateCommonDispy)y <- estimateTagwiseDispy) 计算差异基因 et<- exactTesty)topTagset) 广义线性模型（Glm）
对于更复杂的实验设计（有多个因素），可以通过广义线性模型来及性能你和。
计算离散度
通过下面来估计common离散度、trended离散度、tagwise离散度。
y<- estimateDispy,design)
或者分开依次进行
y <- estimateGLMCommonDispy, design)y <- estimateGLMTrendedDispy, design) #估计trended离散度y <- estimateGLMTagwiseDispy, design)
关于design的构建?)：
#group是区别分组的factor对应不同样本，lane是不同样本的的lane道design <- model.matrix~0+group+lane) colnamesdesign) <- gsub”group”, “”, colnamesdesign)) 计算DE基因
edgrR涉及到差异表达分析的函数有很多： exactTest、glmFit、glmLRT、glmQLFit、glmQLFTest。
qCML估计离散度需要搭配 exact test 进行差异表达分析，对应 exactTest 函数。
而其他四个glm*都是与GLM模型搭配使用的函数。其中，glmFit 和 glmLRT 函数是配对使用的，用于 likelihood ratio test 似然比检验)，而 glmQLFit和 glmQLFTest则配对用于 quasi-likelihood F test 拟极大似然F检验)。
有两个方法glmQLFit)和glmQLFTest),在两个检验方法中，首选QLFit，因为它反映了估计每个基因的离散度时的不确定性。当重复次数较少时，它提供了更强大和可靠的错误率控制
# 例子group <- factorc1,1,2,2,3,3))design <- model.matrix~group)fit <- glmQLFity, design) 如果没有重复样本？
首先确定一个良好的BCV，人类数据一般设定为0.4
一个例子：
bcv <- 0.2counts <- matrix rnbinom40,size=1/bcv^2,mu=10), 20,2)y <- DGEListcounts=counts, group=1:2)et <- exactTesty, dispersion=bcv^2)
如果是人类数据：
y_bcv <- ybcv <- 0.4et <- exactTesty_bcv,dispersion = bcv ^2)
注：
如果差异基因过少，可以尝试下调BCV
输出结果 result = topTagset, n = nrowet$table))$table 查看统计 summaryde <- decideTestsDGEet))
我想建立并管理一个高质量的生信&统计相关的微信讨论群，如果你想参与讨论，可以添加我微信： veryqun 。我会拉你进群，当然有问题也可以微信咨询我。
因为我的笔记分布CSDN、简书、知乎专栏等比较零散，管理起来比较麻烦，思考再三申请了一个 微信公众号，会更加方便地发布更多有关生信息、统计方面内容，如果你觉得有需要可以尝试性和我缔结契约 ^ . ^ 。公众号如下：
对我唯一的精神鼓励可能就是下方的点赞了吧 ^ ^