转录组入门（7）：差异表达分析-CSDN博客

这个步骤推荐在R里面做，载入表达矩阵，然后设置好分组信息，统一用DEseq2进行差异分析，当然也可以走走edgeR或者limma的voom流程。
基本任务是得到差异分析结果，进阶任务是比较多个差异分析结果的异同点。

数据填坑
理论基础：线性模型，设计矩阵和比较矩阵
标准化一二事
探索性分析一二事
使用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/SRR3589956.count",
                       sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589957.count",
                    sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589958.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）分析不同肥料处理后的产量差异，实验设计如下