繪制FeaturePlot時，遇到基因在所有細(xì)胞中表達(dá)水平相同展示效果不理想的情況唠叛，本文引入函數(shù)tryCatch()旨在解決上述問題，并將警告信息保存到日志文件中便于后續(xù)追蹤沮稚。
1 加載R包
library(easypackages)
packages <- c('ggplot2', 'cowplot', 'Seurat')
libraries(packages)
2 挑選所有細(xì)胞中表達(dá)水平相同的基因
引入內(nèi)置數(shù)據(jù)pmbc_small
pbmc_small
## An object of class Seurat
## 230 features across 80 samples within 1 assay
## Active assay: RNA (230 features, 20 variable features)
## 2 dimensional reductions calculated: pca, tsne
# 從全部基因集中挑選在所有細(xì)胞中表達(dá)量相同的基因
object_seurat <- pbmc_small
dat <- as.matrix(object_seurat@assays$RNA@counts)
dat_t <- t(dat)
search_same_value_row <- function(i){
if(all(dat_t[,i]==dat_t[,i][1] )){
print(colnames(dat_t)[i])
}
}
my_genes <- purrr::map_dfc(1:dim(dat_t)[2], search_same_value_row)
gene_same.value = as.character(my_genes)
# 結(jié)果表明：不存在所有細(xì)胞表達(dá)水平都為0的基因Ｒ照印！壮虫！
# 取子集再挑選基因
sub_cells <- subset(pbmc_small, cells = sample(Cells(pbmc_small), 20))
object_seurat <- sub_cells
dat <- as.matrix(object_seurat@assays$RNA@counts)
dat_t <- t(dat)
dim(dat_t)
## [1] 20 230
search_same_value_row <- function(i){
if(all(dat_t[,i]==dat_t[,i][1] )){
print(colnames(dat_t)[i])
}
}
my_genes <- purrr::map_dfc(1:dim(dat_t)[2], search_same_value_row)
## [1] "NCF1"
## [1] "CD180"
## [1] "IGLL5"
## [1] "DLGAP1-AS1"
## [1] "ZNF76"
## [1] "PTPN22"
## [1] "VSTM1"
## [1] "CD1C"
gene_same.value = as.character(my_genes)
# 結(jié)果表明：存在少數(shù)在所有細(xì)胞表達(dá)水平相同的基因０南帷！囚似！
# 高可變基因集
gene_highly = VariableFeatures(sub_cells)[! VariableFeatures(sub_cells) %in% as.character(my_genes)]
3 基因表達(dá)水平的可視化
seurat_object <- sub_cells
gene_set <- c(gene_highly[c(13,18)], gene_same.value[1:2])
feature_plot_fun <- function(gene_set){FeaturePlot(seurat_object , gene_set, cols = c("lightgrey", "blue"), pt.size=2, reduction="tsne")}
VlnPlot_plot_fun <- function(gene_set){VlnPlot(seurat_object, gene_set, pt.size=2)+ NoLegend() + ggtitle(label=gene_set)}
feature_plot <- purrr::map(gene_set, feature_plot_fun)
VlnPlot_plot <- purrr::map(gene_set, VlnPlot_plot_fun)
featureplot1_cluster <- CombinePlots(plots=feature_plot, nrow=1)
VlnPlot_plot_cluster <- CombinePlots(plots=VlnPlot_plot, nrow=1)
plot_grid(plotlist=list(VlnPlot_plot_cluster, featureplot1_cluster), nrow=2)
對比小提琴圖可以看出剩拢，當(dāng)基因在所有細(xì)胞中表達(dá)水平相同時，即使表達(dá)量都為零卻高亮顯示饶唤，容易對實(shí)際表達(dá)解讀造成誤解徐伐，影響了可視化效果，故引入函數(shù)tryCatch()募狂。
4 tryCatch容錯函數(shù)
try就像一個網(wǎng)办素，把try{}里面的代碼所跑出的異常都網(wǎng)住，然后把異常就給catch{}里面的代碼去執(zhí)行祸穷，最后執(zhí)行finally之中的代碼性穿。無論try中代碼有沒有異常，也無論catch是否被異常捕獲到雷滚，finally中的代碼都一定會被執(zhí)行需曾。
有時需要判斷一行命令運(yùn)行的狀態(tài)，然后再做出反應(yīng)祈远，整體來說：
1 是否出現(xiàn)warning呆万，出現(xiàn)了怎么處理？
2 是否出現(xiàn)Error车份，出現(xiàn)了怎么處理谋减？
3 沒有出現(xiàn)怎么處理？
tryCatch({
命令
}, warning = function(w){
# 這里是出現(xiàn)warning狀態(tài)時扫沼，應(yīng)該怎么做出爹，可以用print打印出來，可以執(zhí)行其它命令
}, error = function(e){
# 這里是出現(xiàn)Error狀態(tài)時充甚，應(yīng)該怎么做以政，可以用print打印出來，也可以執(zhí)行其它命令
},finally = {
# 這里是運(yùn)行正常時伴找，應(yīng)該怎么做盈蛮，可以用print打印出來，也可以執(zhí)行其它命令
})
## NULL
5 保存警告信息到日志文件中
# 創(chuàng)建空日志文件
file.create('my_log.txt')
## [1] TRUE
log.path = 'my_log.txt'
feature_plot_fun <- function(gene_set){
tryCatch({
f1 <- FeaturePlot(seurat_object, gene_set, cols = c("lightgrey", "blue"), pt.size=2, reduction="tsne")
}, warning = function(w){
# 這里是出現(xiàn)warning狀態(tài)時技矮，應(yīng)該怎么做抖誉，可以用print打印出來殊轴，可以執(zhí)行其它命令
# print(paste('warning:', w))
f2 <- FeaturePlot(seurat_object, gene_set, cols = c("lightgrey", "lightgrey"), pt.size=2, reduction="tsne")
write(toString(w), log.path, append=TRUE) # 保存警告信息到日志文件中
return(f2)
})
}
6 再次基因表達(dá)水平的可視化
feature_plot <- purrr::map(gene_set, feature_plot_fun)
VlnPlot_plot <- purrr::map(gene_set, VlnPlot_plot_fun)
featureplot1_cluster <- CombinePlots(plots=feature_plot, nrow=1)
VlnPlot_plot_cluster <- CombinePlots(plots=VlnPlot_plot, nrow=1)
plot_grid(plotlist=list(VlnPlot_plot_cluster, featureplot1_cluster), nrow=2)
References
[1] R語言tryCatch使用方法：判斷Warning和Error: http://blog.sciencenet.cn/blog-2577109-1251678.html
[2] Basic Error Handing in R with tryCatch(): https://www.r-bloggers.com/2020/10/basic-error-handing-in-r-with-trycatch/
[3] Feature Plot with 0 count data: https://github.com/satijalab/seurat/issues/1557
轉(zhuǎn)載來自：單細(xì)胞天地
