前段時間宏蛉，一個單細胞分析同行提問阵具，發(fā)現(xiàn)AverageExpression()和通過aggregate計算得到的基因表達均值，兩者的結(jié)果是不一樣的噪叙，疑惑問題出在哪里。
這個問題我們或多或少都有過疑問霉翔，自己手動計算基因在每個cluster的表達均值睁蕾，繪制heatmap熱圖；跟AverageExpression()處理后繪制DoHeatmap會有差異债朵。網(wǎng)上網(wǎng)友也提過類似的問題子眶，也不知其解。決心看下AverageExpression()的源碼序芦。
一臭杰、測試AverageExpression()函數(shù)
library(Seurat)
table(Idents(seurat_obj_Bcell))
# B0 B1 B2 B3 B4 B5 B6 B7 B8
# 2573 2115 2093 1536 1355 1344 1005 92 74
cluster.averages <- AverageExpression(seurat_obj_Bcell, slot = 'data', return.seurat = TRUE)
# assays有三種模式：RNA, SCT, integrated
# -----------------------------------
dim(cluster.averages@assays$RNA@data)
# [1] 21972 9
head(cluster.averages@assays$RNA@data)
# B0 B1 B2 B3 B4
# AL627309.1 0.004422988 0.003237424 0.001007899 0.004426136 0.001509931
# AL627309.5 0.015784045 0.026414631 0.032074415 0.023697387 0.006019230
# AP006222.2 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
# LINC01409 0.008551480 0.009181925 0.019694557 0.019493293 0.025692960
# FAM87B 0.001375785 0.002404144 0.001298187 0.001751876 0.000000000
# LINC01128 0.138134556 0.118331411 0.149791178 0.063821427 0.072906583
# B5 B6 B7 B8
# AL627309.1 0.001075568 0.003654187 0.00000000 0.0000000
# AL627309.5 0.011408478 0.029097945 0.00000000 0.0000000
# AP006222.2 0.000000000 0.000000000 0.00000000 0.0000000
# LINC01409 0.013514985 0.003152831 0.06921692 0.0000000
# FAM87B 0.000000000 0.000000000 0.00000000 0.0000000
# LINC01128 0.054664322 0.153354677 0.02744229 0.1327957
write.csv(cluster.averages@assays$RNA@data, file = "cluster.averages_assays_RNA_data.csv")
#-----------------------------------
dim(cluster.averages@assays$SCT@data)
# [1] 16717 9
head(cluster.averages@assays$SCT@data)
# B0 B1 B2 B3 B4
# AL627309.1 0.0019413712 0.0009451797 0.000477669 0.0019512201 0.0007377352
# AL627309.5 0.0061991675 0.0117510165 0.011873804 0.0122939109 0.0029476808
# LINC01409 0.0031043874 0.0047169899 0.007141187 0.0116506172 0.0110092855
# FAM87B 0.0003885759 0.0009451797 0.000477669 0.0006508298 0.0000000000
# LINC01128 0.0551916324 0.0596507087 0.058458200 0.0351811146 0.0376583238
# LINC00115 0.0177200322 0.0150167064 0.019868203 0.0225307245 0.0204535984
# B5 B6 B7 B8
# AL627309.1 0.000743771 0.0009945302 0.00000000 0.00000000
# AL627309.5 0.006674107 0.0089153637 0.00000000 0.00000000
# LINC01409 0.005934736 0.0009945302 0.03208831 0.00000000
# FAM87B 0.000000000 0.0000000000 0.00000000 0.00000000
# LINC01128 0.023530497 0.0523375251 0.02150621 0.07796154
# LINC00115 0.011834458 0.0187289851 0.02150621 0.05264373
write.csv(cluster.averages@assays$SCT@data, file = "cluster.averages_assays_SCT_data.csv")
#-----------------------------------
dim(cluster.averages@assays$integrated@data)
# [1] 2000 9
head(cluster.averages@assays$integrated@data)
# B0 B1 B2 B3 B4
# LYZ 0.005967343 0.01146686 0.010868760 6.508555e-03 0.012815378
# S100A9 0.009264568 0.03849551 0.005806276 1.388812e-05 0.018757573
# S100A8 0.003377878 0.02276050 0.010335783 1.019810e-03 0.009345999
# IGKV3-15 1.374223462 0.96204743 1.358024467 1.397423e+00 1.096290152
# IGKV3-20 1.529314027 0.96532312 1.316039517 1.776763e+00 1.215281468
# IGKV3-11 1.290120895 0.57287190 1.234111953 9.728570e-01 1.339886382
# B5 B6 B7 B8
# LYZ 0.007843731 0.009692197 -0.005189707 0.002062152
# S100A9 0.006850602 0.008352554 -0.005198381 0.002006344
# S100A8 0.001704852 0.003578382 -0.004440559 0.015471037
# IGKV3-15 1.402736193 0.711017421 0.445023238 0.727528712
# IGKV3-20 1.456745150 0.853644113 1.736057719 1.514083481
# IGKV3-11 1.137085218 0.554071284 1.533300791 0.738172971
write.csv(cluster.averages@assays$integrated@data, file = "cluster.averages_assays_integrated_data.csv")
二、測試AverageExpression()源碼
從GitHub下載seurat源碼谚中，找到AverageExpression()函數(shù)的代碼渴杆，把該函數(shù)的調(diào)用函數(shù)代碼也找齊，運用自己的數(shù)據(jù)集宪塔，一步步執(zhí)行磁奖。我們的目的是計算基因在每個cluster(B0-B8)的平均表達值。
step1: 調(diào)用AverageExpression()
AverageExpression()調(diào)用內(nèi)部函數(shù)()蝌麸，把自身的參數(shù)傳遞給它点寥。
# step1: AverageExpression函數(shù)是傳遞變量給PseudobulkExpression
object <- seurat_obj_Bcell
pb.method = 'average'
assays = NULL
features = NULL
return.seurat = TRUE
group.by = 'ident'
add.ident = NULL
slot = 'data'
verbose = TRUE
PseudobulkExpression(
object = object,
pb.method = 'average',
assays = assays,
features = features,
return.seurat = return.seurat,
group.by = group.by,
add.ident = add.ident,
slot = slot,
verbose = verbose,
...
)
step2: PseudobulkExpression()
查看PseudobulkExpression()函數(shù)，pb.method參數(shù)只能選'average' 或者'aggregate'方法来吩；aggregate方法先不考慮敢辩，只考慮average方法蔽莱。該函數(shù)代碼較多，只提取重要代碼和參數(shù)戚长。
操作1：獲取細胞的cluster矩陣信息data盗冷；即獲取每個細胞對應的cluster信息，‘sample1-AAACGGGCACCTTGTC’細胞的cluster為B4同廉。data為12187*1的矩陣仪糖，測試數(shù)據(jù)共12187個細胞。
assays
# [1] "RNA" "integrated" "SCT"
data <- FetchData(object = object, vars = rev(x = group.by))
dim(data)
# [1] 12187 1
head(data)
# ident
# sample1-AAACGGGCACCTTGTC B4
# sample1-AAACGGGGTAACGACG B2
# sample1-AAAGATGGTACCGCTG B6
# sample1-AAAGATGGTGACGCCT B4
# sample1-AAAGTAGAGTATCGAA B2
# sample1-AAAGTAGGTCTCATCC B3
# data是每個細胞對應的cluster信息(B1-B8)迫肖，共12187個細胞锅劝；
# 去除為na的數(shù)據(jù)；該案例無NA數(shù)據(jù)蟆湖；
data <- data[which(rowSums(x = is.na(x = data)) == 0), , drop = F]
dim(data)
# [1] 12187 1
操作2：cluster矩陣數(shù)據(jù)data轉(zhuǎn)成稀疏矩陣category.matrix故爵，列為cluster的列名（B0, B1, B2...B8）,共9列。data將由12187*1的矩陣轉(zhuǎn)成12187*9的矩陣隅津。
library(Matrix)
category.matrix <- sparse.model.matrix(object = as.formula(
object = paste0(
'~0+',
paste0(
"data[,",
1:length(x = group.by),
"]",
collapse = ":"
)
)
))
dim(category.matrix)
# [1] 12187 9
head(category.matrix)
# 6 x 9 sparse Matrix of class "dgCMatrix"
# data[, 1]B0 data[, 1]B1 data[, 1]B2 data[, 1]B3 data[, 1]B4 data[, 1]B5
# 1 . . . . 1 .
# 2 . . 1 . . .
# 3 . . . . . .
# 4 . . . . 1 .
# 5 . . 1 . . .
# 6 . . . 1 . .
# data[, 1]B6 data[, 1]B7 data[, 1]B8
# 1 . . .
# 2 . . .
# 3 1 . .
# 4 . . .
# 5 . . .
# 6 . . .
# category.matrix是轉(zhuǎn)成一個稀疏矩陣诬垂；
# head(data)
# ident
# sample1-AAACGGGCACCTTGTC B4
# sample1-AAACGGGGTAACGACG B2
# sample1-AAAGATGGTACCGCTG B6
# sample1-AAAGATGGTGACGCCT B4
# sample1-AAAGTAGAGTATCGAA B2
# sample1-AAAGTAGGTCTCATCC B3
# 將12187*1的矩陣轉(zhuǎn)成12187*9的矩陣，
# 如sample1-AAACGGGCACCTTGTC為B4伦仍，即data[, 1]B4為1结窘，其余為0；
# 如sample1-AAACGGGCACCTTGTC為B2充蓝，即data[, 1]B2為1隧枫，其余為0；
# 對每列(B0-B8)求和谓苟，統(tǒng)計屬于該cluster的細胞總數(shù)
colsums <- colSums(x = category.matrix)
colsums
# data[, 1]B0 data[, 1]B1 data[, 1]B2 data[, 1]B3 data[, 1]B4 data[, 1]B5
# 2573 2115 2093 1536 1355 1344
# data[, 1]B6 data[, 1]B7 data[, 1]B8
# 1005 92 74
table(data$ident)
# B0 B1 B2 B3 B4 B5 B6 B7 B8
# 2573 2115 2093 1536 1355 1344 1005 92 74
# 數(shù)據(jù)過濾悠垛，如果某cluster沒有任一細胞，剔除該cluster
category.matrix <- category.matrix[, colsums > 0]
colsums <- colsums[colsums > 0]
操作3：稀疏矩陣category.matrix進行average操作娜谊，計算每個細胞在cluster的占比。例如斤讥，屬于B0這個cluster的細胞群有2573個細胞纱皆，那每個細胞的占比是1/2573=0.00038865。
想一想為什么要這么操作芭商？
if (pb.method == 'average') {
category.matrix <- sweep(
x = category.matrix,
MARGIN = 2,
STATS = colsums,
FUN = "/")
}
1/colsums
# data[, 1]B0 data[, 1]B1 data[, 1]B2 data[, 1]B3 data[, 1]B4 data[, 1]B5
# 0.0003886514 0.0004728132 0.0004777831 0.0006510417 0.0007380074 0.0007440476
# data[, 1]B6 data[, 1]B7 data[, 1]B8
# 0.0009950249 0.0108695652 0.0135135135
操作4：默認slot = "data"派草，獲取歸一化之后每個細胞的基因表達矩陣，矩陣大小為21972*12187铛楣。
# 如果為"RNA"近迁，
data.use <- GetAssayData(
object = object,
assay = "RNA",
slot = "data"
)
# data.use為一個21972*12187的基因表達矩陣，21972為基因個數(shù)簸州，12187為細胞個數(shù)鉴竭；slot = "data"表明是歸一化之后的基因表達矩陣歧譬；
dim(data.use)
# [1] 21972 12187
data.use[1:20,1:20]
操作5：如slot = "data"，基因表達矩陣需進行expm1處理搏存。這個很關(guān)鍵瑰步。
if (slot[i] == 'data') {
data.use <- expm1(x = data.use)
}
data.use[1:20,1:20]
問題1：為什么要進行expm1處理？
數(shù)據(jù)均一化的時候進行了log1p處理璧眠， expm1和log1p是相對的缩焦；
操作6：計算最終的每個細胞的基因在cluster的平均表達量；即基因表達矩陣data.use和細胞在cluster的占比矩陣category.matrix乘積责静。
兩個矩陣的乘積袁滥，要求第一個矩陣的行數(shù)data.use與第二個矩陣category.matrix列數(shù)相等。
data.return[[i]] <- as.matrix(x = (data.use %*% category.matrix))
問題2：為什么通過源碼推算的基因平均表達量矩陣和AverageExpression()計算的不完全一致灾螃？
猜測是category.matrix等矩陣小數(shù)點保留位數(shù)不一致题翻，導致后面幾位數(shù)不一致；
問題3：pb.method有'average' 或者'aggregate'方式睦焕，是否可以切換到aggregate方式：
否藐握，AverageExpression()沒有外置參數(shù)，調(diào)用PseudobulkExpression()函數(shù)時垃喊，已設(shè)置為pb.method = 'average'猾普，不能切換成aggregate
問題4：回歸到最初的提問，AverageExpression()函數(shù)本谜，計算出來的結(jié)果與aggregate得到的均值不一樣初家？
嘗試進行expm1()處理；
group.by <- "recluster"
mat <- as.data.frame(t(as.matrix(GetAssayData(object, assay = "RNA", slot = "data"))))
mat <- expm1(x = mat)
mat <- aggregate(mat, by=list(object@meta.data[[group.by]]), FUN="mean")
rownames(mat) <- mat$Group.1
mat <- t(mat[,-1])
head(mat)
比較expm1()處理后的aggregate得到的均值和AverageExpression()計算的均值乌助，是一樣的溜在；
通過上面的測試，AverageExpression()函數(shù)他托，計算出來的結(jié)果與aggregate得到的均值不一樣的原因是AverageExpression()對歸一化的表達矩陣進行了expm1()處理掖肋。
