绘制FeaturePlot时，遇到基因在所有细胞中表达水平相同展示效果不理想的情况，本文引入函数tryCatch()旨在解决上述问题，并将警告信息保存到日志文件中便于后续追踪。
1 加载R包
library(easypackages)
packages <- c('ggplot2', 'cowplot', 'Seurat')
libraries(packages)
2 挑选所有细胞中表达水平相同的基因
引入内置数据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
# 从全部基因集中挑选在所有细胞中表达量相同的基因
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)
# 结果表明：不存在所有细胞表达水平都为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)
# 结果表明：存在少数在所有细胞表达水平相同的基因！！！
# 高可变基因集
gene_highly = VariableFeatures(sub_cells)[! VariableFeatures(sub_cells) %in% as.character(my_genes)]
3 基因表达水平的可视化
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)
图片
对比小提琴图可以看出，当基因在所有细胞中表达水平相同时，即使表达量都为零却高亮显示，容易对实际表达解读造成误解，影响了可视化效果，故引入函数tryCatch()。
4 tryCatch容错函数
try就像一个网，把try{}里面的代码所跑出的异常都网住，然后把异常就给catch{}里面的代码去执行，最后执行finally之中的代码。无论try中代码有没有异常，也无论catch是否被异常捕获到，finally中的代码都一定会被执行。
有时需要判断一行命令运行的状态，然后再做出反应，整体来说：
1 是否出现warning，出现了怎么处理？
2 是否出现Error，出现了怎么处理？
3 没有出现怎么处理？
tryCatch({
命令
}, warning = function(w){
# 这里是出现warning状态时，应该怎么做，可以用print打印出来，可以执行其它命令
}, error = function(e){
# 这里是出现Error状态时，应该怎么做，可以用print打印出来，也可以执行其它命令
},finally = {
# 这里是运行正常时，应该怎么做，可以用print打印出来，也可以执行其它命令
})
## NULL
5 保存警告信息到日志文件中
# 创建空日志文件
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){
# 这里是出现warning状态时，应该怎么做，可以用print打印出来，可以执行其它命令
# 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 再次基因表达水平的可视化
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
转载来自：单细胞天地
