2025年转录组入门（7）：差异表达分析

大家好，我是讯享网，很高兴认识大家。

数据填坑
理论基础：线性模型，设计矩阵和比较矩阵
标准化一二事
探索性分析一二事
使用DESeq2进行差异基因分析
使用edgeR进行差异基因分析
使用limma进行差异基因分析
- 不同软件包分析结果比较
使用GFOLD进行无重复样本的差异基因分析
不同差异表达分析的比较

数据填坑

原先三个样本的HTSeq-count计数的数据可以在我的GitHub中找到，但是前面已经说过Jimmy失误让我们分析的人类就只有3个样本，另外一个样本需要从另一批数据获取（请注意batch effect)，所以不能保证每一组都有两个重复。

我一直坚信”你并不孤独“这几个字，遇到这种情况的人肯定不止我一个，于是我找到了几种解决方法

使用edgeR，指定dispersion值
无重复转录组数据推荐用同济大学的GFOLD

以上方法都会在后续进行介绍，但是我们DESeq2必须得要有重复的问题亟待解决，没办法我只能自己瞎编了。虽然是编，我们也要有模有样，不能直接复制一份，要考虑到高通量测序的read是默认符合泊松分布的。我是这样编的。

计算KD重复组的均值差，作为泊松分布的均值
使用概率函数rpois()随机产生一个数值，前一步的均值作为lambda，
对一些read count 低于均值的直接加上对应KD重复组之间的差值

# import data if sample are small options(stringsAsFactors = FALSE) control <- read.table("F:/Data/RNA-Seq/matrix/SRR.count", sep="\t", col.names = c("gene_id","control")) rep1 <- read.table("F:/Data/RNA-Seq/matrix/SRR.count", sep="\t", col.names = c("gene_id","rep1")) rep2 <- read.table("F:/Data/RNA-Seq/matrix/SRR.count", sep="\t",col.names = c("gene_id","rep2")) # merge data and delete the unuseful row raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id") raw_count_filt <- raw_count[-1:-5,] ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id) row.names(raw_count_filt) <- ENSEMBL the sample problem delta_mean <- abs(mean(raw_count_filt$rep1) - mean(raw_count_filt$rep2)) sampleNum <- length(raw_count_filt$control) sampleMean <- mean(raw_count_filt$control) control2 <- integer(sampleNum) for (i in 1:sampleNum){ if(raw_count_filt$control[i] < sampleMean){ control2[i] <- raw_count_filt$control[i] + abs(raw_count_filt$rep1[i] - raw_count_filt$rep2[i]) } else{ control2[i] <- raw_count_filt$control[i] + rpois(1,delta_mean) } } # add data to raw_count raw_count_filt$control2 <- control2

讯享网

这仅仅是一种填坑的方法而已，更好模拟数据的方法需要参阅更加专业的文献， 有生之年 我希望能补上这一个部分。

理论基础：线性模型，设计矩阵和比较矩阵

这部分内容最先在 RNA-Seq Data Analysis 的8.5.3节看到，刚开始一点都不理解，但是学完生物统计之后，我认为这是理解所有差异基因表达分析R包的关键。

讯享网

基本上，统计课都会介绍如何使用t检验用来比较两个样本之间的差异，然后在样本比较多的时候使用方差分析确定样本间是否有差异。当然前是样本来自于正态分布的群体，或者随机独立大量抽样。

对于基因芯片的差异表达分析而言，由于普遍认为其数据是服从正态分布，因此差异表达分析无非就是用t检验和或者方差分析应用到每一个基因上。高通量一次性找的基因多，于是就需要对多重试验进行矫正，控制假阳性。目前在基因芯片的分析用的最多的就是limma。

但是，高通量测序(HTS)的read count普遍认为是服从泊松分布（当然有其他不同意见），不可能直接用正态分布的t检验和方差分析。当然我们可以简单粗暴的使用对于的非参数检验的方法，但是统计力不够，结果的p值矫正之估计一个差异基因都找不到。老板花了一大笔钱，结果却说没有差异基因，是个负结果，于是好几千经费打了水漂，他肯定是不乐意的。因此，还是得要用参数检验的方法，于是就要说到方差分析和线性模型之间的关系了。

线性回归和方差分析是同一时期发展出的两套方法。在我本科阶段的田间统计学课程中就介绍用方差分析（ANOVA）分析不同肥料处理后的产量差异，实验设计如下

肥料	重复1	重复2	重复3	重复4
A1	...	...	...	...
A2	...	...	...	...
A3	...	... ...	...	...

这是最简单的单因素方差分析，每一个结果都可以看成 yij = ai + u + eij，其中u是总体均值，ai是每一个处理的差异，eij是随机误差。

2025年转录组入门（7）：差异表达分析

目录

数据填坑

理论基础：线性模型， 设计矩阵和比较矩阵

相关推荐

理论基础：线性模型，设计矩阵和比较矩阵