使用Excel对read counts表达矩阵进行简单的差异分析

2022-10-01 本文已影响0人就是大灯不用装

获得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进行快速填充。

批量填充的公式为“=LN(X)”，这里counts为0的数据转换后结果为-Inf。

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值

使用exp()转换回e次幂值

计算出每个样本的sizefactor之后，将该样本原始的表达量除以该样本的sizefactor, 就得到了归一化之后的表达量。

使用=E2/$AA$55338，让每列counts除以固定的单元格——即之前算得的sizefactor。

得到的归一化的counts值

4. 计算basemean和FoldChange

随后，依照对照组和处理组的分组，计算归一化后的counts的均值basemean，并进行组间的均值比较。

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,为显著性升高。

得到的FoldChange比值

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”的标准来筛选显著性差异的基因。

得到的火山图