DEseq的问题

单细胞转录组数据校正批次效应实战（上）
简而言之，不同时间、不同操作者、不同试剂、不同仪器导致的实验误差，反映到细胞的表达量 上就是批次效应，这个很难去除但可以缩小。如果效应比较还可以接受，但是批次效应很严重， 就可能会和真实的生物学差异相混淆，让结果难以捉摸。我们需要辨别到底存在多大程度的批次 效应，对我们真实的生物学样本会不会产生影响。
校正批次效应的目的就是：减少batch之间的差异，尽量让多个batch的数据重新组合在一起，这 样下游分析就可以只考虑生物学差异因素

专用的标准化方法：之前听过的CPM/FPKM/TPM/TMM都是适用于bulk RNA-seq，分析时也经 常移植到单细胞数据。但是单细胞数据有自己的特点，例如存在细胞差异和基因差异两类系统 偏差，而上述方法主要考虑了基因差异。为了校正细胞间的差异， scran 包特意利用去卷积法 (deconvolution) 开发了 computeSumFactors 函数(Lun, Bach, and Marioni 2016)。它将聚类后的多 组细胞合并在一起屏蔽0值分散的问题，并采用类似CPM的方法计算标准化因子(size factor)
补充：CPM为原始reads除以一个样品总的可用reads数乘以 1,000,000 ，但这种方法容易 受到极高表达且在不同样品中存在差异表达的基因的影响，有点"牵一发动全身"的感觉

scran package implements a variant on CPM specialized for single-cell data (L. Lun, Bach, and Marioni 2016). Briefly this method deals with the problem of vary large numbers of zero values per cell by pooling cells together calculating a normalization factor (similar to CPM) for the sum of each pool. Since each cell is found in many different pools, cell-specific factors can be deconvoluted from the collection of poolspecific factors using linear algebra.

DESeq2分析转录组之预处理
前面导入了DESeq2需要的数据 （未归一化）

导入数据后首先初步过滤(Pre-filtering)
这不是必须，只是为了减少后面计算量，主要就是去除表达量非常少的行，比如设置阈值为每行表达量为10

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
告诉DESeq2你想怎么比较
默认情况下，R会根据字母表顺序排列因子型变量。如果不告诉DESeq2哪组作为对照组的话，它会自动根据字母表顺序设置，这里有两种解决方案：

明确指出contrast 是怎么设置(下面进行详细探讨)

在design中设置factor level，必须在分析之前就设置好，不能分析进行中或完成后再设置，它的设置方法也有两种

利用factor：

dds$condition<-factor(dds$condition, levels = c(“untreated”,“treated”))

看下前后对比

ddsKaTeX parse error: Expected 'EOF', got '#' at position 11: condition #̲ 之前不设置时，默认contr…condition, levels = c(“untreated”,“treated”)) # 设置后，control组(untreated)排在了前面
[1] treated treated treated untreated untreated untreated untreated
Levels: untreated treated
利用relevel

dds$condition<-relevel(dds$condition, ref = “untreated”)
另外，这个因子型的比较公式不是一成不变的，如果我们从原来的表达矩阵中去除了某些样本，那么比较公式中也应该去除，利用droplevels()可以很方便去掉没有样本的factor level

比如这里去掉了第一个样本，那么最后的比较公式就是

droplevels(dds[,-1]$condition)
[1] treated treated untreated untreated untreated untreated
Levels: untreated treated
合并技术重复到同一个表达矩阵中[可选]
技术重复technical replicate就是指一个文库的多个测序的run

提供了collapseReplicates函数可以做，但是不要用这个去合并生物学重复

差异表达分析
标准的差异分析流程被写进了单独的函数DESeq ，然后结果利用results函数生成，提取出log2FC、P值、校正后的p值等6列。

如果没有额外的参数，log2FC和P值是默认对design公式中的最后一个变量或者最后一个因子与参考因子进行比较，举个例子

之前构建的时候，选择的design是condition这个因子型变量

dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = ~ condition)

condition又被手动设置成了参考在前，比较在后的样子

dds$condition
[1] treated treated treated untreated untreated untreated untreated
Levels: untreated treated

这里results利用的就是最后一个因子(treated)与参考因子(untreated)进行比较，算出的log2FC也就是log2(treated/untreated)

!!!这里不能反，因为即使反过来也能得到结果，但与事实就是完全相反的

dds <- DESeq(dds)
res <- results(dds)

截取头信息看一看，的确是treated vs untreated

res
log2 fold change (MLE): condition treated vs untreated
Wald test p-value: condition treated vs untreated
DataFrame with 9921 rows and 6 columns

如果还不放心或者进行更正，可以手动设置(方法二选一)

res <- results(dds, name=“condition_treated_vs_untreated”)
res <- results(dds, contrast=c(“condition”,“treated”,“untreated”))
log2FC(LFC) vs lfcShrink
The shrunken fold changes are useful for ranking genes by effect size and for visualization. --Michael Love (the developer of deseq2)

关于lfcShrink：New function lfcShrink() in DESeq2 这个函数2年前就已经开发出来了

为何采用lfcShrink？log2FC estimates do not account for the large dispersion we observe with low read counts.因此，两种数据特别需要：低表达量占比高的；数据特别分散的

As with the shrinkage of dispersion estimates, LFC shrinkage uses information from all genes to generate more accurate estimates.

举一个来自 DESeq2 paper的例子

左图中相同两组（6J和2J）中绿色和紫色基因的平均值是一样的，但是绿色基因的方差更小，离散程度更小，因此右图中unshrunken LFC estimate (绿色实线) 与 shrunken LFC estimate (绿色虚线)是基本一致的；但是紫色基因由于高离散度导致这两者差别较大。另外可以看出，即使两个基因的标准化count值相似，但也可能存在不同的LFC shrinkage

注意：Shrinking the log2 fold changes不会改变显著差异的基因总数

例如，如果要根据LFC值提取差异基因，需要shrunken values

另外，进行功能分析例如GSEA时，需要提供shrunken values

作者还是非常推荐使用lfcShrink的https://support.bioconductor.org/p/95695/

DESeq2的老版本1.16.0之前是默认进行shrink的，但是并不是所有数据都是需要shrink的，因此作者把这个功能设成了默认关闭 。把原来DESeq中自带的 log2 fold change shrinkage移到了lfcShrink中，可以手动添加[最新版的DESeq2是1.22]

如果使用的是旧版本的，并且不需要shrink，可以这样关闭：

DESeq(dds, betaPrior=FALSE)

这个结果和下面直接计算的结果是一样的

log2 (normalized_counts_group1 / normalized_counts_group2)

对于新版本而言，可以这样手动操作

resultsNames(dds)

[1] “Intercept” “condition_treated_vs_untreated”

resLFC <- lfcShrink(dds, coef=“condition_treated_vs_untreated”, type=“apeglm”)

选择的apeglm参数(Zhu, Ibrahim, and Love 2018)进行effect size shrinkage,which improves on the previous estimator

resLFC相比res数据更加紧凑

发现shrink前后少了一列stat

names(resLFC)
[1] “baseMean” “log2FoldChange” “lfcSE” “pvalue”
[5] “padj”
names(res)
[1] “baseMean” “log2FoldChange” “lfcSE” “stat”
[5] “pvalue” “padj”
开启多线程
上面的分析步骤大约只需要30s，但是对于复杂的design公式、大样本数据，就可以采用多线程加速。DESeq、results、lfcShrink 都可以通过加载BiocParallel实现

library(“BiocParallel”)

先预定4个核，等需要的时候直接使用parallel=TRUE来调用

register(MulticoreParam(4))
探索p值与adjusted p值
作者给出的建议是：Need to filter on adjusted p-values, not p-values, to obtain FDR control. 10% FDR is common because RNA-seq experiments are often exploratory and having 90% true positives in the gene set is ok

先从小到大排个序：

resOrdered <- res[order(res$pvalue),]

order函数是给出从小到大排序后的位置，默认decreasing = FALSE

利用summary可以探索一些基本的统计量

summary(res)

out of 9921 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 518, 5.2%
LFC < 0 (down) : 536, 5.4%
outliers [1] : 1, 0.01%
low counts [2] : 1539, 16%
(mean count < 6)
[1] see ‘cooksCutoff’ argument of ?results
[2] see ‘independentFiltering’ argument of ?results
看看有多少p值小于0.1的？

sum(res$padj < 0.1, na.rm=TRUE) # 有多少p值小于0.1的？
[1] 1054

RNA-seq中的基因表达量计算和表达差异分析
差异分析的步骤：
1）比对；
2） read count计算；
3） read count的归一化；
4）差异表达分析；
表达量计算的本质
目标基因表达量相对参照系表达量的数值。
参照的本质：
（ 1）假设样本间参照的信号值应该是相同的；
（ 2）将样本间参照的观测值校正到同一水平；
（ 3）从参照的数值，校正并推算出其他观测量的值。
例如：Qpcr:目标基因表达量（循环数）相对看家基因表达量（循环数）；RNA-seq:目标基因的表达量（测序reads数），相对样本RNA总表达量（总测序量的reads数），这是最常用的标准。
归一化的原因及处理原则：
1）基因长度
2）测序量
3）样本特异性（例如，细胞mRNA总量，污染等）前两者使用普通的RPKM算法就可以良好解决，关键是第三个问题，涉及到不同的算法处理。
RNA-Seq归一化算法的意义：
基因表达量归一化：在高通量测序过程中，样品间在数据总量、基因长度、基因数目、高表达基因分布甚至同一个基因的不同转录本分布上存在差别。因此不能直接比较表达量，必须将数据进行归一化处理。
RNA-seq差异表达分析的一般原则
1）不同样品的基因总表达量相似
2）上调差异表达与下调差异表达整体数量相似（上下调差异平衡）
3）在两组样品中不受处理效应影响的基因， 表达量应该是相近的（差异不显著）。
4）看家基因可作为表达量评价依据（ 待定）

不同的算法比较：
以什么数值来衡量表达量：RPKM、FPKM、TPM
以什么作为参照标准：TMM（edgeR软件）、De seq矫正
RPKM：是Reads Per Kilobase per Million mapped reads的缩写，代表每百万reads中来自于某基因每千碱基长度的reads数。

本质：
1）以reads数为计算单位；
2）对基因长度（基因间的比较）和总数据量（样本间的比较）做矫正；
RPKM的弊端
1）由于可变剪切，同一基因有效转录区域长度未必相同（这个一般情况下可以不考虑，了解一下：Cufflinks软件考虑了这个问题）优化策略：外显子或转录本水平的表达量分析。
2） 使用reads数计算基因表达量有轻微误差（这里暂不展开，主要了解一下定义）优化策略：FPKM或 TPM
3） mRNA的总量未必相等。
RPKM的优化：FPKm
F = Fragment，即测序片段数量。这些片段都是从完整的cDNA打碎而来的；
本质：以文库中的片段数量为计算单位在Paired-end测序中，一个fragment就是两条PE reads构成的片段。由于是PE比对，理论上比SE比对更可靠。
RPKM的优化：TPM
T = Transcripts
本质：以转录本的条数为计算单位。使用转录本的条数（或者说：转录本的测序深度），代替reads数，在一定条件下定量更准，尤其样本间表达基因总数差异很大的时候（例如，对照样本有1万个基因表达，另外处理组仅有4000个基因表达）。
mRNA总量未必相等
mRNA总量不等——细胞本身不同
例如：活跃组织vs休眠的组织；癌细胞vs正常细胞
mRNA总量不等——污染
例如：核糖体污染外源RNA污染
解决方法——不同算法比较
其中归一化算法介绍：
1）Total Count（TC）：总reads数矫正
2）Upper Quartile（UQ）：上四分之一分位数（总reads）
矫正
3）Median（Med）；中位数（总reads数）矫正
4）Quantile (Q)：基因芯片软件limma中的校正算法；
5）RPKM：总reads数，但引入了基因长度
6）几何平均数：Deseq软件中的算法；
7）TMM：edgeR软件中的算法；
8）RPKM
逻辑1：不同位置数值的稳定性不同

四分位数quartile:将数据按从小到大排列，并分成四等分，这样得到3个分割点，第一个分割点叫做lowerquartile，第二个叫Media，第三个叫Upper quartile
很显然，极大值具有极大不稳定性，而且可能会显著影
响总体之和（假设，我们之中有个马云，我们的总收入
有什么变化？）
所以，Upper quartile和Median的数值，比总表达量之
和更加稳定，更适合作为参照。
逻辑2：表达量居中的基因的表达量值，其数值应该是相似的。
DESeq与edgeR，默认情况下都使用这一的逻辑校正。（DESeq and edgeR Bioconductor packages）
Deseq：异常高表达的基因，会显著影响细胞中的总mRNA的数量。类似的，如果样本中受到不同程度的外源RNA，如病毒、真菌等的污染，也会显著影响样本总mRNA数，导致RPMK值的误差。对于这样的问题，Deseq尝试对数据进行矫正（矫正因子），使表达量处于中间位置的基因表达量应该是基本相同的（即使用表达量处于中间的基因表达量值作为参照，而减少高表达基因的作用）。
Deseq： 校正因子=样本表达中位数/所有样本表达量中位数：回答了一个关键的问题：Deseq不同差异比较组间，计算得到的表达量值不同。因
为样本在变化，“所有样本表达量的中位数”也在变动。RPKM：总表达量为参照
Deseq：中位数为参照

TMM（edgeR）：与Deseq类似，在去除高表达基因和差异最大的基因后，TMM也是要找到一个加权系数，使剩余的基因在被矫正后差异倍数可能小。TMM的加权系数是基于两两样本比较后推算获得的（也就是两组样本的比较，将产生与这次比较相关的加权系数）。然后将所有基因除以这个加权系数，从而保证大部分表达量居中的基因表达量最相似。
不同RNA-seq表达量归一化算法的区别
Deseq类的校正算法：理论上更加稳定；但不同批次的比较会得到不同的表达量值，不利于进行多处理组/批次数据的统一分析（例如，趋势分析、共表达分析）校正会掩盖一些问题（例如：样本污染）
RPKM类的算法： 容易受异常高表达基因、外源污染等的干扰；但也更容易从结果的异常中，发现潜在问题；得到的表达量值是恒定的，多处理组/批次的数据可以合并分析。折中的方法：使用RPKM类的算法，但需要人工检查数据是否
异常。备注： Deseq软件也可以关闭校正的功能。

实际经验总结
总之：从多方面考虑，RPKM类算法，如果合理使用，依然是最优的。具体问题具体分析：在遇到问题的时候，找到问题的来源，从而给出解决方案（没有完美的流程，只有最佳解决方案）

DESeq2和EdgeR都可用于做基因差异表达分析，主要也是用于RNA-Seq数据，同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据。

这两个都属于R包，其相同点在于都是对count data数据进行处理，都是基于负二项分布模型。因此会发现，用两者处理同一组数据，最后在相同阈值下筛选出的大部分基因都是一样的，但是有一部分不同应该是由于其估计离散度的不同方法所导致的。

5数据预处理
5.1原始数据尺度转换
由于更深的测序总会产生更多的序列，在差异表达相关的分析中，我们很少使用原始的序列数。在实践中，我们通常将原始的序列数进行归一化，来消除测序深度所导致的差异。通常被使用的方法有基于序列的CPM（counts per million）、log-CPM、FPKM（fragments per kilobase of transcript per million），和基于转录本数目的RPKM（reads per kilobase of transcript per million）。

尽管CPM和log-CPM转换并不像RPKM和FPKM那样考虑到基因长度区别的影响，但在我们的分析中经常会用到它们。虽然也可以使用RPKM和FPKM，但CPM和log-CPM只使用计数矩阵即可计算，且已足以满足我们所感兴趣的比较的需要。假设不同条件之间剪接异构体（isoform）的使用没有差别，差异表达分析研究同一基因在不同条件下的表达差异，而不是比较多个基因之间的表达或测定绝对表达量。换而言之，基因长度在我们感兴趣的比较中始终不变，且任何观测到的差异是来自于条件的变化而不是基因长度的变化。

在此处，使用edgeR中的cpm函数将原始计数转换为CPM和log-CPM值。其中，在进行对数转换前会给所有基因的计数加上2，以避免对零取对数，且可减小低表达基因之间的差异。 如果可以提供基因长度信息，可以使用edgeR中的rpkm函数计算RPKM值，就像计算CPM值那样简单。
5.2删除低表达基因
所有数据集中都混有表达的基因与不表达的基因。尽管我们想要检测在一种条件中表达但再另一种条件中不表达的基因，也有一些基因在所有样品中都不表达。实际上，这个数据集中19%的基因在所有九个样品中的计数都是零。

我们应当忽略在任何条件都没有表达到具有生物学意义的程度的基因，这样可以减小基因的子集而保留感兴趣的基因，且可以减小下游在进行差异表达测试时的计算量。 通过检查log-CPM值，可以看出每个样本中很大一部分基因都是不表达或低表达的（如下图中A部分所示）。在此使用CPM值为1作为标准（相当于log-CPM值为0），将表达量高于此阈值的基因看作表达，反之看作不表达。我们只保留在至少一个组（或者至少在整个实验的三个样品中）表达的基因用于下游分析。

尽管任何合理的值都能用作表达的阈值，一般我们在分析中使用CPM值为1，因为它对于大多数数据集都能很好地分隔表达的基因与不表达的基因。在这里，CPM值为1意味着如果一个基因在测序深度最低的样品中（JMS9-P8c, 序列数量约2千万）有至少20个计数，或者在测序深度最高的样品中（JMS8-3，序列数量约7.6千万）有至少76个计数，那么它是“表达的”。如果测序的序列片段是在外显子水平而不是基因水平统计的，且／或实验的测序深度很低，可能需要考虑更低的CPM阈值。

5.3归一化基因表达分布
在样品制备或测序过程中，不具备生物学意义的外部因素会影响单个样品的表达。比如说，在实验中第一批制备的样品会总体上表达高于第二批制备的样品。假设所有样品的表达值的范围和分布都应当相似，需要进行归一化来确保整个实验中每个样本的表达分布都相似。

密度图和箱线图等展示每个样品基因表达量分布的图表可以用于判断是否有样品和其他样品分布有差异。在此数据集中，所有样品的log-CPM分布都很相似（上图B部分）。

尽管如此，我们依然需要使用edgeR中的calcNormFactors函数，用TMM(Robinson and Oshlack 2010)方法进行归一化。此处计算得到的归一化系数被用作文库大小的缩放系数。当我们使用DGEList对象时，这些归一化系数被自动存储在x$samples$norm.factors。对此数据集而言，TMM归一化的作用比较温和，这体现在所有的缩放因子都相对接近1。

简单使用DESeq2/EdgeR做差异分析
Posted: 五月 07, 2017 Under: Transcriptomics By Kai no Comments

DESeq2和EdgeR都可用于做基因差异表达分析，主要也是用于RNA-Seq数据，同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据。

这两个都属于R包，其相同点在于都是对count data数据进行处理，都是基于负二项分布模型。因此会发现，用两者处理同一组数据，最后在相同阈值下筛选出的大部分基因都是一样的，但是有一部分不同应该是由于其估计离散度的不同方法所导致的。

DESeq2的使用方法：
输入矩阵数据，行名为sample，列名为gene；DESeq2不支持无生物学重复的数据，因此我选择了2个样本，3个生物学重复的数据；并对count data取整（经大神指点，这里需要说明下，我的测试数据readcount是RSEM定量的结果，并不是常见的htseq-count的结果，所以count值会有小数点，而DESeq2包不支持count数有小数点，所以这里需要round取整）。

database_all <- read.table(file = “readcount”, sep = “\t”, header = T, row.names = 1)
database <- database_all[,1:6]
#type <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
database <- round(as.matrix(database))
设置分组信息以及构建dds对象

condition <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
coldata <- data.frame(row.names = colnames(database), condition)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
使用DESeq函数进行估计离散度，然后进行标准的差异表达分析，得到res对象结果

dds <- DESeq(dds)
res <- results(dds)
最后设定阈值，筛选差异基因，导出数据

table(res$padj<0.05)res<-res[order(res$padj),]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by=“row.names”,sort=FALSE)
write.csv(resdata,file = “LC_1_vs_LC_2.csv”)
EdgeR的使用方法：
跟DESeq2一样，EdgeR输入矩阵数据，行名为sample，列名为gene；DESeq2不支持无生物学重复的数据，因此我选择了2个样本，3个生物学重复的数据。

exprSet_all <- read.table(file = “readcount”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
exprSet <- exprSet_all[,1:6]
group_list <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
设置分组信息，去除低表达量的gene以及做TMM标准化

exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]
exprSet <- DGEList(counts = exprSet, group = group_list)
exprSet <- calcNormFactors(exprSet)
使用qCML（quantile-adjusted conditional maximum likelihood）估计离散度（只针对单因素实验设计）

exprSet <- estimateCommonDisp(exprSet)
exprSet <- estimateTagwiseDisp(exprSet)
寻找差异gene(这里的exactTest函数还是基于qCML并且只针对单因素实验设计)，然后按照阈值进行筛选即可

et <- exactTest(exprSet)
tTag <- topTags(et, n=nrow(exprSet))
tTag <- as.data.frame(tTag)
write.csv(tTag,file = “LC_1_vs_LC_2_edgeR.csv”)
Summary
以上我主要针对单因素两两比较组进行差异分析，其实DESeq2和EdgeR两个R包都可以对多因素进行差异分析。

DESeq2修改以上代码的分组信息design参数以及在差异分析results函数中添加所选定的分组因素，其他代码基本一样，具体参照DESeq2手册

EdgeR则需要用Cox-Reid profile-adjusted likelihood (CR)方法来估算离散度，y <- estimateDisp(y, design)或者分别使用三个函数（y <- estimateGLMCommonDisp(y, design)，y <- estimateGLMTrendedDisp(y, design)， ）y <- estimateGLMTagwiseDisp(y, design)；然后差异表达分析也跟单因素分析不同，主要使用generalized linear model (GLM) likelihood ratio test 或者 quasi-likelihood(QL) F-test，具体代码可以参照EdgeR手册。