在單細(xì)胞測(cè)序下游分析中，當(dāng)重點(diǎn)關(guān)注哪些基因在所有細(xì)胞平均表達(dá)顯著時(shí)，可選取所選取的top基因進(jìn)行可視化。
0、scater包方法
在scaterR包里，有個(gè)函數(shù)plotHighestExprs，可自動(dòng)根據(jù)SingleCellExperiment對(duì)象數(shù)據(jù)進(jìn)行繪制。但有個(gè)缺點(diǎn)就是太慢了，具體操作如下
library(SingleCellExperiment)
sce <- SingleCellExperiment(
assays = list(counts = expr)
)
class(sce)
library(scater)
sce <- logNormCounts(sce)
p1 <- plotHighestExprs(sce, exprs_values = "logcounts",n=20)
#因?yàn)槌鰣D比較慢，故保存為本地后再觀察可提高效率
ggsave("./p1.pdf", plot = p1)
- 下面介紹一種手動(dòng)繪制箱圖的方法，不僅出圖快，而且箱圖反映的信息也更多一些
1、獲取單細(xì)胞表達(dá)矩陣
具體可分為三種情況
- 直接獲得原始表達(dá)矩陣；
dim(expr)
expr[1:4,1:4]
- 從Seurat對(duì)象獲??；
scRNA = CreateSeuratObject(counts=expr)
class(scRNA)
#[1] "Seurat"
#attr(,"package")
#[1] "Seurat"
scRNA@assays$RNA@counts[1:4,1:4]
counts <- scRNA@assays$RNA@counts[1:4,1:4]
注意的是Seurat是以稀疏矩陣
dgCMatrix格式儲(chǔ)存的
- 從SingleCellExpriment對(duì)象獲??；
sce <- SingleCellExperiment(
assays = list(counts = expr)
)
class(sce)
#[1] "SingleCellExperiment"
#attr(,"package")
#[1] "SingleCellExperiment"
assay(sce, "counts")[1:4,1:4]
counts <- assay(sce, "counts")[1:4,1:4]
2、原始表達(dá)矩陣標(biāo)準(zhǔn)化
最主要的原始是，不同細(xì)胞的文庫(kù)大小可能因測(cè)序過(guò)程存在差異。通過(guò)標(biāo)準(zhǔn)化，可使基因在不同細(xì)胞的表達(dá)情況具有可比性。
hist(colSums(expr))
- 手動(dòng)對(duì)原始矩陣標(biāo)準(zhǔn)化
cell_lib <- colSums(expr) #計(jì)算細(xì)胞文庫(kù)
normalized <- t(t(expr)/cell_lib)*100
上述兩步也可合并成一步，但可能跳的有點(diǎn)多，關(guān)鍵理解如下示例的第三步的除法：
（1）R語(yǔ)言的計(jì)算是向量化的；例如可以進(jìn)行向量間加減乘除運(yùn)算，具體規(guī)則，自己嘗試下理解更深刻；
（2）應(yīng)用到矩陣時(shí)，可以理解為一行代表向量的一個(gè)元素。
# 最后乘100是將小數(shù)更友好得表示（百分比）
normalized[1:4,1:4]
- 如果是Seurat對(duì)象、SingleCellExperiment對(duì)象，則可以運(yùn)用相應(yīng)的函數(shù)計(jì)算，再導(dǎo)出標(biāo)準(zhǔn)化矩陣。
scRNA <- NormalizeData(scRNA)
scRNA@assays$RNA@data[1:4,1:4]
如下圖，結(jié)果與我們上面計(jì)算的有點(diǎn)不同；
查看NormalizeData幫助文檔可知，其默認(rèn)方法是計(jì)算每個(gè)細(xì)胞中基因表達(dá)量與文庫(kù)的比值，然后乘一個(gè)size.factor(一般是10000)，最后進(jìn)行l(wèi)og轉(zhuǎn)換（加1，避免0以及零點(diǎn)幾的無(wú)意義結(jié)果log1p）
library("scater")
sce <- logNormCounts(sce)
assay(sce, "logcounts")[1:4,1:4]
如上圖結(jié)果，與之前兩種方法結(jié)果還不一樣。
scater包的
logNormCounts的主要思路是首先平衡每個(gè)細(xì)胞的文庫(kù)大小，再計(jì)算相對(duì)于所在細(xì)胞文庫(kù)的比例。具體流程可見(jiàn)下圖的示例。
綜上不同的方法，可得到不同的標(biāo)準(zhǔn)化結(jié)果。具體采用哪種方法，可自己根據(jù)情況選擇。
3、選取top基因，進(jìn)行箱圖可視化
- 這里我選擇了第一種，自己簡(jiǎn)單進(jìn)行的標(biāo)準(zhǔn)化的結(jié)果。
#選擇平均在每個(gè)細(xì)胞表達(dá)最顯著的基因
most_exp <- sort(rowSums(normalized),T)[20:1] / ncol(expr)
# 可以想一下為什么設(shè)置[20:1]，與下面繪圖有關(guān)
- boxplot
boxplot(t(normalized[names(most_exp),]),
cex=.01, las=1, cex.axis = .5,
#分別設(shè)置離群點(diǎn)大小，軸標(biāo)簽的大小及方向
xlab="% total count per cell",
col=scales::hue_pal()(20)[20:1], #配色
horizontal=TRUE)
- ggplot2
library(ggplot2)
library(reshape2)
newdata <- melt(t(normalized[names(most_exp),]))
head(newdata)
ggplot(data=newdata, aes(x=Var2, y=value, fill=Var2)) +
geom_boxplot(outlier.size = 0.5) +
coord_flip() +
theme(legend.position="none") +
labs(title="High expression gene", x="",
y="% total count per cell")
最后總結(jié)下：主要根據(jù)單細(xì)胞表達(dá)矩陣，繪制箱圖來(lái)可視化高表達(dá)基因。中間介紹了下如何讓同一基因在不同細(xì)胞表達(dá)具有可比性的幾種方法。
參考文章
https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_01_qc.html
