現(xiàn)在回想單細(xì)胞項(xiàng)目的一些常規(guī)分析疏橄，clusters/groups間的差異表達(dá)分析（differential expression analysis）幾乎是分析必備項(xiàng)焕数，F(xiàn)indAllMarkers()和FindMarkers()的使用頻率非常高瀑凝。做差異表達(dá)分析的工具很多彪杉，我們?nèi)绾胃鶕?jù)自己的數(shù)據(jù)和實(shí)驗(yàn)處理方案來選擇合適的分析工具亩码，還是有些不清晰瞬场。
這次想通過梳理FindAllMarkers源碼，回答以下幾個(gè)問題：
1.scRNA-seq數(shù)據(jù)的哪些方面使差異基因分析復(fù)雜化
2.FindAllMarkers算法原理
3.FindAllMarkers使用的注意事項(xiàng)
一句喜、cluster差異表達(dá)分析
對(duì)cluster進(jìn)行差異表達(dá)分析预愤，使我們能夠找到特定于每個(gè)cluster的標(biāo)記基因，也有助于cluster的細(xì)胞類型鑒定咳胃。
根據(jù)cluster的定義植康，它受cluster本身定義方式的影響，因此在找到標(biāo)記之前確定正確的cluster分辨率很重要展懈。 如果某些cluster群缺少顯著性標(biāo)記销睁，可以嘗試調(diào)整聚類數(shù)目。
1.1 差異表達(dá)分析方法
- FindAllMarkers：比較一個(gè)cluster與所有其他cluster之間的基因表達(dá)
-
FindMarkers：比較兩個(gè)特定cluster之間的基因表達(dá)
運(yùn)行上面的函數(shù)存崖，會(huì)為每個(gè)cluster生成marker基因列表冻记，從而獲得一個(gè)cluster相對(duì)于其他cluster的表達(dá)顯著上調(diào)基因（up-regulated）和下調(diào)基因（down-regulated）。在分析中来惧，我們經(jīng)常設(shè)置FindAllMarkers的參數(shù)only.pos=True冗栗，只顯示當(dāng)前cluster陽性表達(dá)的基因。高表達(dá)的marker gene，有助于我們識(shí)別cluster細(xì)胞類型贞瞒，和后續(xù)的差異基因富集通路分析等偶房。
1.2. seurat執(zhí)行差異表達(dá)分析
根據(jù)單細(xì)胞數(shù)據(jù)預(yù)處理的方式不同（lognormalize和sctransform），執(zhí)行差異表達(dá)分析代碼有所不同军浆。
經(jīng)sctransform預(yù)處理的單細(xì)胞數(shù)據(jù)棕洋，進(jìn)行差異表達(dá)分析：
- 1）作者推薦使用在RNA assay上做差異表達(dá)分析，而不是 SCTransform乒融。
切換到RNA assay掰盘，對(duì)原始的基因counts矩陣進(jìn)行NormalizeData和ScaleData，獲取歸一化的seurat_obj[['RNA']]@data后赞季，進(jìn)行差異基因分析愧捕。 - 2）但也可在SCT assay上做差異表達(dá)分析
# lognormalize預(yù)處理
DefaultAssay(pbmc) <- "RNA"
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
# sctransform預(yù)處理
# method1
DefaultAssay(pbmc) <- "RNA"
pbmc <- NormalizeData(pbmc)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
all.markers.RNA <- FindAllMarkers(pbmc, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5)
# 使用SCT assay進(jìn)行差異表達(dá)分析
# method2
DefaultAssay(srat) <- "SCT"
all.markers.SCT <- FindAllMarkers(pbmc, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5)
1.3 cluster差異基因列表
Differentially expressed gene (DEG)：差異表達(dá)基因
生成的差異基因列表列名說明：
| 名稱 | 解釋 |
|---|---|
| gene | 基因名稱 |
| cluster | 該基因?qū)?yīng)的cluster |
| pct.1 | 在當(dāng)前cluster細(xì)胞中檢測(cè)到該基因表達(dá)的細(xì)胞比例 |
| pct.2 | 在其它c(diǎn)luster細(xì)胞中檢測(cè)到該基因表達(dá)的細(xì)胞比例 |
| avg_logFC | 兩組間平均logFC，Seurat v4默認(rèn)log2申钩。正值表示特征在第一組中表達(dá)得更高 |
| p_val | 未調(diào)整P-value次绘，數(shù)值越小越顯著 |
| p_val_adj | 基于使用數(shù)據(jù)集中所有特征的bonferroni校正，校正后的p值 |
1.4 為什么scRNA-seq的差異表達(dá)分析與bulk RNA-seq不同
由于單細(xì)胞測(cè)序技術(shù)的局限性撒遣，單細(xì)胞測(cè)序數(shù)據(jù)通常具有高噪音邮偎，有較高的dropout問題，即很多低表達(dá)或中度表達(dá)的基因無法有效檢測(cè)到义黎，造成基因表達(dá)矩陣大量的零值禾进。所以，以前針對(duì)傳統(tǒng)bulk RNA-seq測(cè)序數(shù)據(jù)開發(fā)的差異表達(dá)檢測(cè)方法或軟件不一定完全適用于單細(xì)胞測(cè)序數(shù)據(jù)廉涕。
- scRNA-seq 受較高噪聲影響（技術(shù)和生物因素）
- 少量可用的mRNA由于擴(kuò)增偏好性和“dropout 事件”（技術(shù)層面）
mRNA含量低 ->低counts泻云，高dropout率，擴(kuò)增偏好性差異 - 3' 端偏好狐蜕、部分覆蓋和測(cè)序深度不均（技術(shù)層面）轉(zhuǎn)錄的隨機(jī)性（生物層面）
- 基因表達(dá)的多模態(tài)宠纯； 細(xì)胞群中存在多種可能的細(xì)胞狀態(tài)（生物）
1.5 常用的統(tǒng)計(jì)檢驗(yàn)
- 非參數(shù)檢驗(yàn)，如Wilcoxon rank sum test (Mann-Whitney U test)
例如scRNA-seq 中存在dropouts(零值)的情況层释，可能會(huì)失敗 - 特定于scRNA-seq的方法征椒，如MAST
充分利用每組的大量樣本（細(xì)胞）
MAST考慮了隨機(jī)dropouts 丟失和雙細(xì)胞的情況 - 特定于bulk RNA-seq，如DESeq2
基于負(fù)二項(xiàng)分布湃累，適用于UMI數(shù)據(jù)
注意：不應(yīng)該過濾基因，因?yàn)?DESeq2 通過借用其他具有相似表達(dá)水平的基因的信息來模擬離散度
運(yùn)行非常慢碍讨！ 僅用于比較2個(gè)cluster
| 算法 | 工具 | |
|---|---|---|
| Hurdle model combining tests for difference in mean and detection rate | MAST | Appropriate for read counts or molecule/UMI counts |
| Non-parametric, based on ranking expression | Wilcoxon | Robust and very easy to use |
| Negative binomial model GLM, parametric | DESeq2/edgeR |
1.6 統(tǒng)計(jì)學(xué)中的Wilcoxon rank-sum test
統(tǒng)計(jì)學(xué)中的假設(shè)檢驗(yàn)：用來判斷樣本與樣本治力，樣本與總體的差異是由抽樣誤差引起還是本質(zhì)差別造成的統(tǒng)計(jì)推斷方法〔颍基本是思路類似數(shù)學(xué)中的“反證法”宵统，是先對(duì)總體做出某種假設(shè)，然后通過抽樣研究的統(tǒng)計(jì)推理，對(duì)此假設(shè)應(yīng)該是拒絕或者接受马澈。
在假設(shè)檢驗(yàn)中瓢省，我們?cè)O(shè)定原假設(shè)H0和備選假設(shè)H1：
- 原假設(shè)H0：兩組數(shù)據(jù)沒有顯著性差異，兩組相似
- 原假設(shè)H1：兩組數(shù)據(jù)存在顯著性差異痊班，兩組不同
在假設(shè)檢驗(yàn)中勤婚，先設(shè)定原假設(shè)(H0)，再設(shè)定與其相反的備擇假設(shè)(H1)涤伐。接下來隨機(jī)抽取樣本馒胆，若在原假設(shè)成立的情況下，樣本發(fā)生的概率(P)非常小凝果，說明原假設(shè)不成立祝迂，備擇假設(shè)成立，則拒絕原假設(shè)器净。否則型雳，接受原假設(shè)。
判斷依據(jù)：依據(jù)小概率事件的原理山害，統(tǒng)計(jì)學(xué)上纠俭，把小概率事件在一次實(shí)驗(yàn)中看成是實(shí)際不可能發(fā)生的事件，一般認(rèn)為等于或小于0.05或0.01的概率為小概率粗恢。當(dāng)P值小于或等于顯著性水平時(shí)柑晒，則拒絕原假設(shè)。
Wilcoxon檢驗(yàn)（也稱為 Mann-Withney-Wilcoxon 檢驗(yàn)）：
- 一種非參的假設(shè)檢驗(yàn)方法眷射，這意味著它不依賴于任何特定概率分布的數(shù)據(jù)匙赞，它并不需要正態(tài)性分布的假設(shè)，一般用來檢測(cè) 2 個(gè)數(shù)據(jù)集是否來自于相同分布的總體妖碉。例如涌庭，Student's t 檢驗(yàn)僅適用于數(shù)據(jù)為高斯分布或樣本量足夠大（通常 n≥30）的情況。
- 將數(shù)據(jù)按從小到大欧宜，或等級(jí)從弱到強(qiáng)轉(zhuǎn)換成秩后夭坪，再求秩和戈二，計(jì)算檢驗(yàn)統(tǒng)計(jì)量–秩和統(tǒng)計(jì)量，做出統(tǒng)計(jì)推斷。
- 包含：Wilcoxon rank-sum test（秩和檢驗(yàn)）和 Wilcoxon signed-rank test（符號(hào)秩檢驗(yàn)）
Wilcoxon秩和檢驗(yàn)：推斷連續(xù)型變量資料或有序變量資料的兩個(gè)獨(dú)立樣本代表的兩個(gè)總體分布是否有差別仁热。
- 通常在樣本的正態(tài)性假設(shè)不成立且樣本量較小時(shí)使用
- 秩和檢驗(yàn)更適用于非成對(duì)數(shù)據(jù)之間的差異性檢測(cè)，不要求被檢驗(yàn)的兩組數(shù)據(jù)包含相同個(gè)數(shù)的元素
- 非參數(shù)統(tǒng)計(jì)假設(shè)檢驗(yàn)淑掌，用于評(píng)估兩組獨(dú)立樣本是否存在顯著差異
image.png
Wilcoxon秩和檢驗(yàn)步驟：
從每個(gè)總體中選擇隨機(jī)樣本吊圾，設(shè) n1 和 n2 分別是較小樣本和較大樣本中的觀察數(shù)；
編秩：將兩組數(shù)據(jù)n1+n2個(gè)觀測(cè)值混合挂绰，由小到大統(tǒng)一編秩屎篱；不同組遇到相同數(shù)據(jù)取平均秩次；
計(jì)算各組秩和；
-
建立雙尾檢驗(yàn)假設(shè)交播，確定檢驗(yàn)水準(zhǔn)
- H0: 兩組樣本的總體中位數(shù)等于0
- H1: 兩組樣本的樣本的中位數(shù)不等于0
- α =0.05
-
根據(jù)分組樣本數(shù)重虑、各組秩和查詢臨界值表格，得出統(tǒng)計(jì)學(xué)結(jié)論秦士；
- 利用統(tǒng)計(jì)軟件或查Wilcoxon符號(hào)秩檢驗(yàn)的分布表缺厉，
-
如果P值較小（比如小于或等于給定的顯著性水平伍宦，如α=0.05）芽死，則可以拒絕原假設(shè)；如果P值較大次洼，則沒有充分的證據(jù)來拒絕原假設(shè)关贵。
image.png
1.7 統(tǒng)計(jì)學(xué)中的Bonferroni校驗(yàn)
問題：為什么要進(jìn)行P值校正？
當(dāng)同一研究問題下進(jìn)行多次假設(shè)檢驗(yàn)時(shí)卖毁，不再符合小概率原理所說的“一次試驗(yàn)”揖曾，需要進(jìn)行多重檢驗(yàn)校正。如果在該研究問題下只要有檢驗(yàn)是陽性的亥啦，就對(duì)該問題下陽性結(jié)論的話炭剪，對(duì)該問題的檢驗(yàn)的犯一類錯(cuò)誤的概率就會(huì)增大。如果同一問題下進(jìn)行n次檢驗(yàn)翔脱，每次的檢驗(yàn)水準(zhǔn)為α（每次假陽性概率為α）奴拦，則n次檢驗(yàn)至少出現(xiàn)一次假陽性的概率會(huì)比α大。換句話說届吁，在同時(shí)對(duì)多組數(shù)據(jù)進(jìn)行處理和比較的時(shí)候错妖，很可能其中部分?jǐn)?shù)據(jù)因?yàn)殡S機(jī)效應(yīng)而超過閾值，造成假陽性結(jié)果疚沐。而檢驗(yàn)的次數(shù)越多暂氯，出現(xiàn)假陽性的概率就越大。
Bonferroni校正：最嚴(yán)格的多重檢驗(yàn)校正方法亮蛔，該方法是由Carlo Emilio Bonferroni發(fā)展的痴施，故稱Bonferroni校正。
其校正原理為：在同一數(shù)據(jù)集上同時(shí)檢驗(yàn)n個(gè)相互獨(dú)立的假設(shè)究流，那么用于每一假設(shè)的統(tǒng)計(jì)顯著水平辣吃，應(yīng)為僅檢驗(yàn)一個(gè)假設(shè)時(shí)的顯著水平的1/n。舉個(gè)例子：如要在同一數(shù)據(jù)集上檢驗(yàn)兩個(gè)獨(dú)立的假設(shè)芬探，顯著水平設(shè)為常見的0.05齿尽。此時(shí)用于檢驗(yàn)該兩個(gè)假設(shè)應(yīng)使用更嚴(yán)格的0.025。即0.05* (1/2)灯节。假如有10個(gè)檢驗(yàn)，經(jīng)過Bonferroni 校正以后，新的p值可以這樣計(jì)算1-(1-p)^(1/10)=p'炎疆，p'就是 新的顯著性水平卡骂。
我們要對(duì)數(shù)千個(gè)基因進(jìn)行統(tǒng)計(jì)測(cè)驗(yàn)，所以我們可能只是偶然獲得p-values很小的基因（假陽性基因）形入，但是實(shí)際全跨，在樣本間毫無差異。每個(gè)基因的p值只是單一檢驗(yàn)(單基因)的結(jié)果亿遂。檢驗(yàn)的基因越多浓若，假陽性率就越高。所以蛇数，我們需要對(duì)基因的p值進(jìn)行多重假設(shè)檢驗(yàn)校正挪钓，從而減低假陽性結(jié)果在我們的檢驗(yàn)中出現(xiàn)的次數(shù)。
FindAllMarkers函數(shù)默認(rèn)使用Bonferroni校驗(yàn)耳舅。校驗(yàn)的次數(shù)取決于測(cè)試（基因）的數(shù)目碌上，如果我們測(cè)試較少的基因，校正就沒有那么苛刻浦徊。
Bonferroni correction: adjusted p-value = raw p-value * number of genes tested
二馏予、FindMarkers函數(shù)解析
FindMarkers函數(shù)就是一個(gè)S3泛型函數(shù)；包含Seurat盔性，DimReduc霞丧，Assay等數(shù)據(jù)類型。
# NAMESPACE
S3method(FindMarkers,Assay)
S3method(FindMarkers,DimReduc)
S3method(FindMarkers,Seurat)
S3method(FindMarkers,default)
FindMarkers <- function(object, ...) {
UseMethod(generic = 'FindMarkers', object = object)
}
# differential_expression.R
FindMarkers.default <- function(){...}
FindMarkers.Assay <- function(){...}
FindMarkers.DimReduc <- function(){...}
FindMarkers.Seurat <- function(){...}
案例1：找出seurat對(duì)象seurat_obj中cluster0和cluster1的差異基因冕香。
DefaultAssay(seurat_obj) <- "RNA"
DEG_wilcox <- FindMarkers(seurat_obj, ident.1 = 0, ident.2= 1, test.use = "wilcox", min.pct = 0.5, logfc.threshold = 0.5, only.pos = T)
head(DEG_wilcox)
# p_val avg_log2FC pct.1 pct.2 p_val_adj
# RBP7 0 1.698551 0.726 0.004 0
# PGD 0 1.670035 0.823 0.046 0
# AGTRAP 0 1.488062 0.833 0.136 0
# TNFRSF1B 0 1.696910 0.858 0.119 0
# EFHD2 0 2.097743 0.917 0.079 0
# CDA 0 1.577383 0.760 0.001 0
dim(DEG_wilcox)
# [1] 789 5
2.1 在RNA assay上做差異表達(dá)分析
LogNormalize預(yù)處理后進(jìn)行差異基因分析蛹尝，F(xiàn)indMarkers函數(shù)的代碼很多，我把主要步驟記錄如下：
step1：將以上FindMarkers命令中的參數(shù)傳至FindMarkers.Seurat()函數(shù)暂筝；
object <- seurat_obj
ident.1 = 0
ident.2= 1
test.use = "wilcox"
min.pct = 0.5
logfc.threshold = 0.5
only.pos = T
step2：獲取歸一化的基因表達(dá)矩陣seurat_obj[['RNA']]@data箩言；在FindMarkers.Assay()函數(shù)會(huì)調(diào)用FoldChange()函數(shù)。
data.use <- GetAssayData(object = object, slot = "data")
fc.results <- FoldChange(
object = object,
slot = data.slot,
cells.1 = cells.1,
cells.2 = cells.2,
features = features,
pseudocount.use = pseudocount.use,
mean.fxn = mean.fxn,
fc.name = fc.name,
base = base
)
step3：FoldChange.default()的計(jì)算非常關(guān)鍵焕襟。fold change是組間基因表達(dá)量的差異倍數(shù)陨收，avg_log2FC是基因在組間的平均表達(dá)量取log2，Seurat v4默認(rèn)log2鸵赖。正值表示特征在第一組中表達(dá)得更高务漩。
FoldChange計(jì)算方式如下：
# 基因表達(dá)均值計(jì)算方式
mean.fxn <- mean.fxn %||% switch(
EXPR = slot,
'data' = function(x) {
return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base))
},
'scale.data' = rowMeans,
function(x) {
return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
}
)
step4：每個(gè)基因計(jì)算在兩個(gè)細(xì)胞群的占比 pct.1和 pct.2，計(jì)算方式是：基因陽性表達(dá)細(xì)胞數(shù)/該群細(xì)胞總數(shù)它褪；
pct.1 <- round(
x = rowSums(x = object[features, cells.1, drop = FALSE] > thresh.min) /
length(x = cells.1),
digits = 3
)
pct.2 <- round(
x = rowSums(x = object[features, cells.2, drop = FALSE] > thresh.min) /
length(x = cells.2),
digits = 3
)
每個(gè)基因的avg_log2FC的計(jì)算：
cells.1群細(xì)胞的每個(gè)基因表達(dá)量取expm1饵骨，求平均+1，再取log2茫打；
cells.2群細(xì)胞的每個(gè)基因表達(dá)量取expm1居触，求平均+1妖混，再取log2；
每個(gè)基因的avg_log2FC計(jì)算：fc <- (data.1 - data.2)
data.1 <- mean.fxn(object[features, cells.1, drop = FALSE])
data.2 <- mean.fxn(object[features, cells.2, drop = FALSE])
mean.fxn <- mean.fxn %||% switch(
EXPR = slot,
'data' = function(x) {
return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base))
},
'scale.data' = rowMeans,
function(x) {
return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
}
)
fc <- (data.1 - data.2)
step5：我們進(jìn)入到FindMarkers.default()這個(gè)函數(shù)轮洋，先判斷cluster0和cluster1的細(xì)胞數(shù)是不是滿足最小細(xì)胞數(shù)3制市，min.cells.group默認(rèn)是3個(gè)，cluster0和cluster1最少要3個(gè)以上細(xì)胞弊予；顯然祥楣，符合要求。根據(jù)設(shè)置的最小百分比min.pct 和 log2FC閾值和僅保留陽性基因汉柒，對(duì)所有的36601基因進(jìn)行過濾误褪，僅保留789個(gè)基因；
step6：調(diào)用PerformDE()函數(shù)碾褂，執(zhí)行wilcox檢驗(yàn)和P值校正兽间；
p_val的計(jì)算：test.use = "wilcox"，采用WilcoxDETest檢測(cè)計(jì)算出p_val斋扰；
de.results <- switch(
EXPR = test.use,
'wilcox' = WilcoxDETest(
data.use = data.use,
cells.1 = cells.1,
cells.2 = cells.2,
verbose = verbose
)
)
p_val_adj的計(jì)算：對(duì)p_val進(jìn)行Bonferroni校正
de.results$p_val_adj = p.adjust(
p = de.results$p_val,
method = "bonferroni",
n = nrow(x = object)
)
匯總FindMarkers函數(shù)的主要步驟如下：
2.2 在SCT assay上做差異表達(dá)分析
經(jīng)過SCTransform預(yù)處理后進(jìn)行差異基因分析渡八，看下哪些步驟和上面的不同：
DefaultAssay(seurat_obj) <- "SCT"
DEG_wilcox <- FindMarkers(seurat_obj, ident.1 = 0, ident.2= 1, test.use = "wilcox", min.pct = 0.5, logfc.threshold = 0.5, only.pos = T)
head(DEG_wilcox)
# p_val avg_log2FC pct.1 pct.2 p_val_adj
# RBP7 0 1.475970 0.709 0.004 0
# PGD 0 1.464552 0.804 0.045 0
# AGTRAP 0 1.296850 0.803 0.135 0
# TNFRSF1B 0 1.473898 0.830 0.120 0
# EFHD2 0 1.821583 0.915 0.086 0
# CDA 0 1.373028 0.745 0.001 0
dim(DEG_wilcox)
# [1] 686 5
閱讀完源碼，RNA/SCT assay做差異表達(dá)分析只有輸入的歸一化矩陣不同传货，fold change計(jì)算方法完全一致屎鳍。
小心得：
- 求基因在細(xì)胞群的平均表達(dá)量
通過閱讀源碼，seurat在計(jì)算基因在細(xì)胞群的平均表達(dá)量不是對(duì)細(xì)胞群的每個(gè)細(xì)胞歸一化的表達(dá)值求均值问裕；而是對(duì)將細(xì)胞群的細(xì)胞歸一化的基因表達(dá)量進(jìn)行expm1還原逮壁，再求均值。這也是合理的粮宛，AverageExpression()函數(shù)在求基因的平均表達(dá)也是如此處理的窥淆。
這個(gè)怎么理解呢？
可以把細(xì)胞群看成一個(gè)大細(xì)胞巍杈，包含每個(gè)細(xì)胞校正后的counts忧饭，我們求單個(gè)細(xì)胞的基因表達(dá)量=基因的總counts/細(xì)胞內(nèi)總counts，另外筷畦，每個(gè)細(xì)胞考慮到測(cè)序深度的影響會(huì)乘以scale_factor進(jìn)行校正词裤；
計(jì)算步驟：
- 每個(gè)基因歸一化基因表達(dá)量：log1p(counts/colSums[cell-idx] *scale_factor)
- 進(jìn)行expm1還原：
expm1(log1p(counts/colSums[cell-idx] *scale_factor)) =counts/colSums[cell-idx]*scale_factor
cell1的基因表達(dá)值(expm1還原)=count1/colSums[cell1]*scale_factor
cell2的基因表達(dá)值(expm1還原)=count2/colSums[cell2]*scale_factor
...
cellN的基因表達(dá)值(expm1還原)=countN/colSums[cellN]*scale_factor - =>基因在細(xì)胞群的表達(dá)量=(count1/colSums[cell1]+count2/colSums[cell2]+...countN/colSums[cellN])*scale_factor;
- =>求mean=>基因的平均表達(dá)量
問題：為什么求基因的平均表達(dá)量不需要log1p處理？
基因的平均表達(dá)量在上面的步驟求mean后不需要log1p處理鳖宾，單個(gè)細(xì)胞的基因表達(dá)為了使其滿足正態(tài)分布會(huì)進(jìn)行l(wèi)og1p處理吼砂，但是細(xì)胞群的平均表達(dá)就不需要考慮其正態(tài)性；故不需要進(jìn)行l(wèi)og1p處理鼎文；
如果渔肩，我們對(duì)每個(gè)細(xì)胞歸一化的變量量求均值，就會(huì)差之千里拇惋。
- avg_log2FC的計(jì)算梳理
前面我們說了基因的平均表達(dá)量計(jì)算：
cells.1群細(xì)胞的每個(gè)基因表達(dá)量取expm1周偎，求平均抹剩；
cells.2群細(xì)胞的每個(gè)基因表達(dá)量取expm1，求平均蓉坎；
但是我們是計(jì)算基因的avg_log2FC吧兔；所以要取log2；
那為什么會(huì)出現(xiàn)“均值+1”袍嬉，這個(gè)怎么理解，計(jì)算方式如下：
每個(gè)基因的avg_log2FC的計(jì)算：
cells.1群細(xì)胞的每個(gè)基因表達(dá)量取expm1灶平，求平均+1伺通，再取log2；
cells.2群細(xì)胞的每個(gè)基因表達(dá)量取expm1逢享，求平均+1罐监，再取log2；
每個(gè)基因的avg_log2FC計(jì)算：fc <- (data.1 - data.2)
答：這是避免log2(基因的平均表達(dá)值)為負(fù)瞒爬，
如果基因的平均表達(dá)值=0弓柱，log2(基因的平均表達(dá)值)會(huì)負(fù)；
但是log2(基因的平均表達(dá)值+1)=0侧但，避免計(jì)算結(jié)果為負(fù)矢空。
細(xì)胞群的平均表達(dá)量的方式：total sum scaling (TSS) normalization
The expm1 does un-log the data, but the normalization persists (this would be lost in the counts slot)
expm1() transformed in order to recover normalized values not in log scale
在網(wǎng)上也請(qǐng)教過一個(gè)大佬，以下是他的答復(fù)禀横，非常感謝迷茫時(shí)有人點(diǎn)燈
你要取平均屁药，很自然的應(yīng)該用原始數(shù)據(jù)啊（總共三個(gè)蘋果柏锄，五個(gè)小孩酿箭，每個(gè)小孩0.6個(gè)蘋果）。log以后的數(shù)據(jù)取平均趾娃，你說這是個(gè)什么東西表達(dá)了什么含義呢缭嫡，就很難讓人理解的。（不過你非要log以后取平均抬闷，我猜測(cè)也能得到相似的結(jié)果妇蛀，但是很難讓人理解這東西到底是個(gè)啥……）
