当所有细胞基因表达量相同时如何更好的可视化？

分享是一种态度

绘制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

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注