获得read counts的表达矩阵后，一般来说可以用DESeq2、limma和edgeR等R包进行差异比较分析。
本文将展示利用Excel进行简单的RNA-seq的差异表达分析。
本文使用的差异方法主要基于DESeq2的标准化和差异比较方法，辅以简单的ttest检验，来实现简单的Excel差异表达分析。
以GSE158197为例。
1. Log转换
首先对counts数据进行log转换。DESeq2默认是使用自然对数，对应到Excel即为LN()，也可以使用log2或log10。
单击选中要填充列的第一行单元格后，移动表格到最底部，SHIFT+左键单击，即可选中要填充的整个区域，然后CTRL+D进行快速填充。
2. DESeq2的原理文档说明（可忽略）
需要对数据进行均一化，方法参照DESeq2的estimateSizeFactors
模块介绍：https://bioconductor.org/packages/release/bioc/manuals/DESeq2/man/DESeq2.pdf
########normalizationFactors##########
dds <- makeExampleDESeqDataSet(n=100, m=4)
normFactors <- matrix(runif(nrow(dds)*ncol(dds),0.5,1.5),
ncol=ncol(dds),nrow=nrow(dds),
dimnames=list(1:nrow(dds),1:ncol(dds)))
# the normalization factors matrix should not have 0's in it
# it should have geometric mean near 1 for each row
normFactors <- normFactors / exp(rowMeans(log(normFactors)))
normalizationFactors(dds) <- normFactors
dds <- DESeq(dds)
########estimateSizeFactors##########
dds <- makeExampleDESeqDataSet()
estimateSizeFactorsForMatrix(counts(dds))
geoMeans <- exp(rowMeans(log(counts(dds))))
estimateSizeFactorsForMatrix(counts(dds),geoMeans=geoMeans)
##########
S4 method for signature 'DESeqDataSet'
estimateSizeFactorsForMatrix(
counts,
locfunc = stats::median,
geoMeans,
controlGenes,
type = c("ratio", "poscounts")
)
可以看到，在DESeq2的文档中，标准化的处理要求数据矩阵中不能有0，处理数据是先对数据进行log()，然后获得该行的Mean，随后用exp转换回counts，最后使用原数据normFactors比上exp(rowMeans(log(normFactors)，即得到了标准化的数据。
3. 计算sizefactor
计算单个样本的sizafactor时，将该样本中每个基因的表达量减去对应的所有样本中的均值，然后取中位数。由于开始进行了log转换，最后在转换回来。
回到Excel的操作，则是我们利用AVERAGE函数对所有样本的基因counts均值进行计算，操作同上，键入公式后进行快速填充：
average函数计算10例样本的基因表达均值
随后计算counts数与均值的差值：
利用公式计算Δcounts
随后，根据每列的Δcounts数值，使用MEDIAN计算得到中位数，同时因为存在0值log后产生的“#NUM!”，故使用IF(ISNUMBER())只获取列中的数值。
以本数据集为例，即为 "=MEDIAN(IF(ISNUMBER(AA2:AA55336),AA2:AA55336))"，得到中位数后同样进行快速填充。
提取每个样本的Δcounts的中位数
随后用exp()函数转换回counts值
计算出每个样本的sizefactor之后，将该样本原始的表达量除以该样本的sizefactor, 就得到了归一化之后的表达量。
使用=E2/$AA$55338，让每列counts除以固定的单元格——即之前算得的sizefactor。
4. 计算basemean和FoldChange
随后，依照对照组和处理组的分组，计算归一化后的counts的均值basemean，并进行组间的均值比较。
中位数算法可以进一步降低高差异表达基因的影响，而提高中等表达的基因的地位。
标准化因子的生物学意义
其实这个标准化因子算法就是选出一个有代表性的gene X（其实是每个样本一个代表性gene X），而这个gene X的reads for gene X/average for gene X比值就是标准化因子。只不过选取gene X的时候，通过对数变换和中位数的方法，更多的参考了中等表达基因和管家基因的数据趋势，而剔除了特异性表达基因和高差异表达基因的影响。相比较RPKM,FPKM,TPM标准化方法是除以总Read数，DESeq2标准化方法是除以一个有代表性基因的Read数，只不过这个Read数进行了变换（它除以了几何平均Read数， reads for gene X/average for gene X）。因为更能处理存在特异性表达基因和高差异表达基因的数据。
两组的basemean相除，即得到FoldChange。注：因为我们直接使用的是Counts数据，因此比较得到的是FoldChange而非logFoldChange。
故FoldChange>1即为有表达升高，FoldChange>2 才能被等同视为log2FoldChange>1,为显著性升高。FoldChange<1即为有表达降低，FoldChange<0.5 才能被等同视为log2FoldChange<-1,为显著性升高。
5. 计算Pvalue
DESeq2实际上使用的是负二项分布检验（LRT或Wald），更好的符合RNA-seq样本的实际分布情况。但其过程中有较多中间参数需要计算，故仅使用Excel进行简单的T检验进行验证。
ttest函数对两组数据进行假设检验
TTEST 函数语法具有下列参数：
Array1必需。 第一个数据集。
Array2必需。 第二个数据集。
tails必需。 指定分布尾数。 如果 tails = 1，则 TTEST 使用单尾分布。 如果 tails = 2，则 TTEST 使用双尾分布。
Type必需。 要执行的 t 检验的类型。
因为我们不知道比较的方向性，只检测组间比较间是否存在差异，故选择双尾分布。因为Excel分析简单的测序样本量较少，默认方差不齐，选用双样本异方差假设，即Type=3。
若希望进行方差齐性检验后再确认t检验类型，则可以用以下参数
F=F.TEST(AK2:AO2,AP2:AT2) ######新建一列F检验的结果数值
pvalue=if(e>0.05,T.TEST(AK2:AO2,AP2:AT2,2,2),T.TEST(AK2:AO2,AP2:AT2,2,3))) ####依照新建列的F检验结果，确定选用T检验的类型
也可使用公式“=LOG(FoldChange, 2)”，将其转换为更常见的数据格式。
同时可以使用“=-LOG10(pvalue)”来生成值，以此为Y值，LOG2FoldChange为X值，进行散点图的绘制，即为火山图。
并且可依照“pvalue<0.05”和“log2FoldChange >1 or log2FoldChange <-1”的 标准来筛选显著性差异的基因。
