对样本的无监督聚类
在我们看来，用于检查基因表达分析的最重要的探索性图表之一便是MDS图或其余类似的图。这种图表使用无监督聚类方法展示出了样品间的相似性和不相似性，能让我们在进行正式的检验之前对于能检测到多少差异表达基因有个大致概念。理想情况下，样本会在不同的实验组内很好的聚类，且可以鉴别出远离所属组的样本，并追踪误差或额外方差的来源。如果存在技术重复，它们应当互相非常接近。
这样的图可以用limma中的plotMDS函数绘制。第一个维度表示能够最好地分离样品且解释最大比例的方差的引导性的倍数变化（leading-fold-change），而后续的维度的影响更小，并与之前的维度正交。当实验设计涉及到多个因子时，建议在多个维度上检查每个因子。如果在其中一些维度上样本可按照某因子聚类，这说明该因子对于表达差异有影响，在线性模型中应当将其包括进去。反之，没有或者仅有微小影响的因子在下游分析时应当被剔除。
suppressMessages(library(RColorBrewer))
col.group <- pheno
levels(col.group) <- brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
plotMDS(lcpm, labels=pheno,col = col.group)
title(main="Sample groups")
差异表达分析
主要说明一下edgeR中的glmQLFTest,exactTest以及limma中的voom这几种获取差异基因的不同方式
1基于limma的差异分析
Limma包基于线性模型建模。 它最初设计用于分析微阵列数据，但最近已扩展到RNA-seq数据。 根据limma用户指南的当前建议是使用edgeR包的TMM标准化和“voom”转换，其本质上将标准化数据取对数并估计它们的均值 - 方差关系以确定在线性建模之前每次观察的权重。 默认情况下，Benjamini-Hochberg程序用于估计FDR 。
首先先建立design矩阵,在此研究中，我们想知道哪些基因在我们研究的两组之间以不同水平表达。在我们的分析中，假设基础数据是正态分布的，为其拟合一个线性模型。在此之前，需要创建一个包含细胞类型design矩阵。
design <- model.matrix(~0+pheno)
colnames(design) <- gsub("pheno", "", colnames(design))
rownames(design) <- colnames(y)
#注意使用colnames()，而不能使用names(),后者只会得到：[1] "counts" "samples"
design
据显示对于RNA-seq计数数据而言，当使用原始计数或当其被转换为log-CPM值时，方差并不独立于均值。使用负二项分布来模拟计数的方法假设均值与方差间具有二次的关系。在limma中，假设log-CPM值符合正态分布，并使用由voom函数计算得到的精确权重来调整均值与方差的关系，从而对log-CPM值进行线性建模。
当操作DGEList对象时，voom从x中自动提取文库大小和归一化因子，以此将原始计数转换为log-CPM值。在voom中，对于log-CPM值额外的归一化可以通过设定normalize.method参数来进行。
下图左侧展示了这个数据集log-CPM值的均值-方差关系。通常而言，方差是测序实验中的技术差异和不同细胞类型的重复样本之间的生物学差异的结合，而voom图会显示出一个在均值与方差之间递减的趋势。 生物学差异高的实验通常会有更平坦的趋势，其方差值在高表达处稳定。 生物学差异低的实验更倾向于急剧下降的趋势。
不仅如此，voom图也提供了对于上游所进行的过滤水平的可视化检测。如果对于低表达基因的过滤不够充分，在图上表达低的一端，受到非常低的表达计数的影响，可以观察到方差水平的下降。如果观察到了这种情况，应当回到最初的过滤步骤并提高用于该数据集的表达阈值。
par(mfrow=c(1,2))
v <- voom(y, design, plot=TRUE)
fit <- lmFit(v, design)
head(coef(fit),3)
contr <- makeContrasts(trt - untrt, levels = colnames(coef(fit)))
contr
# Contrasts
#Levels trt - untrt
# trt 1
# untrt -1
diff <- contrasts.fit(fit, contr)
diff <- eBayes(diff)
plotSA(diff, main="Final model: Mean-variance trend")
top.table <- topTable(diff, sort.by = "P", n = Inf)
DEG <- na.omit(top.table)
summary(decideTests(diff))
2基于edgeR的差异基因分析
edgeR使用经验贝叶斯估计和基于负二项模型的精确检验来确定差异基因。 特别地，经验贝叶斯用于通过在基因之间来调节跨基因的过度离散程度。 使用类似于Fisher精确检验但适应过度分散数据的精确检验用于评估每个基因的差异表达。edgeR 在默认情况下，执行TMM标准化程序以考虑样本之间的不同测序深度，通过Benjamini-Hochberg用于控制FDR 。
利用exactTest估计差异基因
精确检验适用于多组实验的精确统计法，使用函数exactTest估计差异基因，人们将其称为classic edgeR。
estimateDisp函数在一组离散网格点上为每个标签(基因)计算一个似然矩阵，然后应用加权似然经验贝叶斯方法获得后验离散度估计。如果没有设计矩阵，它计算每个标签的分位数的条件似然，然后将其最大化。在这种情况下，它类似于函数estimateCommonDisp和estimateTagwiseDisp。
同样首先利用calcNormFactors对因子进行矫正
#对因子矫正
y <- calcNormFactors(y)
#对每个基因估测了一个经验贝叶斯文件离散值，一个公共离散值，一个趋势离散值
y <- estimateDisp(y,design = design)
### 等同于下面两步
# 估计变异系数，即估计方差；估计内部差异程度，看组间差异是否比内部差异大，如果大，可选为差异基因
y <- estimateCommonDisp(y)
# 估计逐段分布
y <- estimateTagwiseDisp(y)
# 绘制基因的趋势图
plotBCV(y)
# 比较获得差异基因
et <- exactTest(y,pair = c("normal","tumor"))
利用glmQLFTest估计差异基因
似然比检验是用广义线性模型（glms）的统计学方法，适用于不同复杂程度的多因素实验，而QLFTest则是首选，因为它反映了估计每个基因离散度的不确定性。当重复次数较少时，它可以提供更可靠的错误率控制。
y <- calcNormFactors(y)
y <- estimateDisp(y, design)
plotBCV(y)
fit <- glmQLFit(y, design, robust=TRUE)
plotQLDisp(fit)
con <- makeContrasts(trt- untrt, levels = design)
res <- glmQLFTest(fit, contrast=con)
summary(decideTestsDGE(res))
plotMD(res)
DEmarkers <- topTags(res, n=Inf)$table
head(DEmarkers)
is.sig <- DEmarkers$FDR <= 0.05
plotSmear(res, de.tags=rownames(DEmarkers)[is.sig], cex=0.1)
