前言 前面我们演示了生信工程师Roman Hillje分享的scRNA-seq workflow的部分（你还缺scRNA-seq的workflow吗？ ）, 今天继续后面的学习。
Roman Hillje的github主页：
https://github.com/romanhaa workflow：
https://romanhaa.github.io/projects/scrnaseq_workflow/ 12 Perform analysis for Cerebro 为了进一步研究Cerebro https://github.com/romanhaa/Cerebro 中的数据集，我们现在添加了更多的meta data，计算标记基因(marker genes)，标记基因的通路富集，并导出数据集。这个过程在很大程度上遵循了cerebroApp网站 https://romanhaa.github.io/cerebroApp/articles/cerebroApp_workflow_Seurat.html 上描述的工作流程。
12.1 Add meta data load("labels_assinged_seurat.Rdata")12.2 Calculate relationship trees 这里，我们为所有三个分组变量:样本、聚类和细胞类型计算关系树——就像我们之前为聚类所做的那样。它们中的每一个都被添加到Seurat对象的@misc槽中，这样它就可以被提取到Cerebro中。
#Samples:12.3 cerebroApp steps 在这里，我们计算线粒体和核糖体转录本的百分比，得到最多表达的基因和标记基因，使用富集过的通路进行测试，并进行基因集富集分析GSEA。因为cereBroAPP作者提桶跑路，已经不支持seurat5了，所以没办法只能自己改函数，改过的函数存在cereBroFunction.R脚本里面
#感觉是cerebroApp不能兼容Seurat5，找不到原来位置的counts(因为v5在layers下面了)cereBroFunction.R:
#函数 addPercentMtRibo function (object, assay = "RNA" , organism ="hg" ,gene_nomenclature = 'name' ){function (object, assay = "RNA" , genes ,message) {function (x) {if (length(genes_here) == 1 ) {"Only one " , message,"gene was found" ))else {#配置文件中记录了线粒体和核糖体基因 "extdata/genes_mt_" , "_" , gene_nomenclature, ".tsv.gz" ), package = "cerebroApp" ), FALSE ) %>% dplyr::select(1 ) %>% "extdata/genes_ribo_" , "_" , gene_nomenclature, ".tsv.gz" ), package = "cerebroApp" ), FALSE ) %>% dplyr::select(1 ) %>% "mitochondrial" , list(genes_mt = genes_mt_here))"ribosomal" , list(genes_ribo = genes_ribo_here))[email protected] $mt_percent <- as.vector(values_mt[[1 ]])[email protected] $ribo_percent <- as.vector(values_ribo[[1 ]])return (object)#函数 getMostExpressedGenes function (object, assay = "RNA" , groups = NULL ) {for (i in seq_along(groups)) {if (is.factor([email protected] [[current_group]])) {[email protected] [[current_group]])else if (is.character([email protected] [[current_group]])) {[email protected] [[current_group]])if (any(is.na(group_levels))) {[email protected] [[current_group]] %>% TRUE ) %>% length()warning (paste0("Found " , number_of_cells_without_group_assignment, " cell(s) without group assignment (NA) for `" , "`. These cells will be ignored during the analysis." ), FALSE )if (length(group_levels) == 1 ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] Group `" , current_group, "` contains only one subgroup. Will calculate most expressed genes " , "across all cells of this group..." ))TRUE )1 :100 ], 1 :100 ])"most_expressed_genes" ]][[current_group]] <- most_expressed_geneselse if (length(group_levels) > 1 ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] Get most expressed genes for " , length(group_levels), " groups in `" , current_group, "`..." ))function (x) {[email protected] )[which([email protected] [[current_group]] == if (is.vector(transcript_count_matrix) == TRUE ) {else {TRUE )1 :100 ], 1 :100 ])"most_expressed_genes" ]][[current_group]] <- most_expressed_genesreturn (object)#函数 getMarkerGenes function (object, organism = NULL , groups = NULL , "cerebro_seurat" , only_pos = TRUE , min_pct = 0.7 , 0.25 , thresh_p_val = 0.01 , test = "wilcox" , TRUE , ... ) {if (is.null(object@misc$marker_genes)) {if (organism %in % c("hg" , "mm" ) == FALSE ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] No information about genes on cell surface because organism is " , "either not specified or not human/mouse." ))else {if (organism == "hg" || organism == "human" ) {"hgnc_symbol" "hsapiens_gene_ensembl" else if (organism == "mm" || organism == "mouse" ) {"external_gene_name" "mmusculus_gene_ensembl" 1 while (!exists("genes_on_cell_surface" ) && attempt <= 3 ) {try (genes_on_cell_surface <- biomaRt::getBM(attributes = temp_attributes, "go" , values = "GO:0009986" , mart = biomaRt::useMart("ensembl" , 1 ])if (!exists("genes_on_cell_surface" )) {"[" , format(Sys.time(), "%H:%M:%S" ), "] Genes in GO term \"cell surface\" (GO:0009986) could not be " , "retrieved, possibly due to the server not being reachable at the " , "moment." ))for (i in seq_along(groups)) {if (is.factor([email protected] [[current_group]])) {[email protected] [[current_group]])else if (is.character([email protected] [[current_group]])) {[email protected] [[current_group]])if (any(is.na(group_levels))) {[email protected] [[current_group]] %>% TRUE ) %>% length()warning (paste0("Found " , number_of_cells_without_group_assignment, " cell(s) without group assignment (NA) for `" , "`. These cells will be ignored during the analysis." ), FALSE )if (length(group_levels) == 1 ) {warning (paste0("Only one group level found for group `" , "`. Will skip this group and proceed to next." ), FALSE )else if (length(group_levels) > 1 ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] Get marker genes for " , length(group_levels), " groups in `" , current_group, "`..." ))... )if (nrow(results) == 0 ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] No marker genes found for any of the level of `" , "`." ))"no_markers_found" else if (nrow(results) > 0 ) {if (exists("genes_on_cell_surface" ) && "gene" %in % in % if ("cluster" %in % colnames(results)) {"gene" ), dplyr::everything())"marker_genes" ]][[name]][[current_group]] <- resultsreturn (object)#函数 performGeneSetEnrichmentAnalysis function (object, assay = "RNA" , GMT_file, NULL , name = "cerebro_GSVA" , 0.05 , thresh_q_val = 0.1 , ... ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] Loading gene sets..." ))function (file) ";" , col_names = c("X1" ), for (i in seq_len(nrow(gmt))) {"\t" )[[1 ]] %>% 3 :length(temp_genes)]"\t" ), "[" , 1 ) %>% unlist(), geneset.description = lapply(strsplit(gmt$X1, "\t" ), "[" , 2 ) %>% unlist())return (gene_set_loaded)NA , NA )for (i in seq_along(gene_sets$genesets)) {"," )"[" , format(Sys.time(), "%H:%M:%S" ), "] Loaded " , " gene sets from GMT file." ))0 )"[" , format(Sys.time(), "%H:%M:%S" ), "] Extracting transcript counts " , "from `data` slot of `" , assay, "` assay..." ))for (i in seq_along(groups)) {if (is.factor([email protected] [[current_group]])) {[email protected] [[current_group]])else if (is.character([email protected] [[current_group]])) {[email protected] [[current_group]])if (any(is.na(group_levels))) {[email protected] [[current_group]] %>% TRUE ) %>% length()warning (paste0("Found " , number_of_cells_without_group_assignment, " cell(s) " , "without group assignment (NA) for `" , "`. These cells will be ignored during the analysis." ), FALSE )"[" , format(Sys.time(), "%H:%M:%S" ), "] Performing analysis for " , "{length(group_levels)} subgroups of group `{current_group}`..." ))if (length(group_levels) == 1 ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] Only one group level found " , "for group `" , "`. Will skip this group and proceed " , "to next." ), call. = FALSE )NULL else if (length(group_levels) > 1 ) {TRUE , simplify = TRUE , future.globals = FALSE , function (x) {[email protected] [[current_group]] == if (length(cells_in_this_group) == 1 ) {else {... )for (i in seq_len(ncol(enrichment_scores))) {if (nrow(temp_results) < 2 ) {"[" , format(Sys.time(), "%H:%M:%S" ), "] Only 1 gene set " , "available, therefore p- and q-values cannot be calculated." ))NA , NA )else if (nrow(temp_results) >= 2 ) {1 ]))1 , ties = min)$lfdr)"name" ) %>% dplyr::select(c("group" , "name" , "description" , "length" , "genes" , "enrichment_score" , "p_value" , "q_value" )) %>% dplyr::rename(`:=`(!!current_group, in % "[" , format(Sys.time(), "%H:%M:%S" ), "] {nrow(results)} gene sets " , "passed the thresholds across all subgroups of group `{current_group}`." ))if (nrow(results) == 0 ) {"no_gene_sets_enriched" "enriched_pathways" ]][[name]][[current_group]] <- resultsreturn (object)13 单基因的表达 13.1单基因的小提琴图 首先，我们发现用宽度而不是面积来衡量小提琴是很有用的。 其次，我们建议对小提琴进行修剪，以避免小提琴延伸到负尺度，这是没有意义的，因为没有任何负的表达值。 第三，在数据中添加一些噪声可能是有意义的，以避免突出显示刻度底部的不同值，例如log(1)， log(2)等。在使用Seurat对表达量进行可视化时，这是默认情况下完成的(甚至没有选项)。 在这种情况下，差异非常细微，如下图所示(红细胞和单核细胞)。
p1 <- [email protected] %>%[email protected] %>%13.2 多基因点图 点图可以展示多个基因的表达，通过一个特定的分组变量。接下来的例子中，我们将使用检测到的细胞类型和聚类的标记基因。
13.2.1 根据细胞类型划分 # cells will be grouped by samples that they have been assigned to[email protected] $cell_type_singler_blueprintencode_main)[email protected] $cell_type_singler_blueprintencode_main == x)[email protected] $cell_type_singler_blueprintencode_main == x)[email protected] %>%13.2.2 根据聚类结果进行划分 group_ids <- levels([email protected] $seurat_clusters)[email protected] $seurat_clusters == x)[email protected] $seurat_cluster == x)[email protected] %>%14 基因集中的表达量 我们根据细胞类型划分检查GO通路B_CELL_ACTIVATION中基因的表达量
library(msigdbr)14.1 Mean log-expression temp_labels <- [email protected] %>%[email protected] $cell_type_singler_blueprintencode_main,14.2 Module score seurat <- AddModuleScore([email protected] %>%[email protected] %>%15 GSEA富集分析 15.1 GSVA 为了说明原因，我们将使用来自MSigDB的50个Hallmark基因集，并将单核细胞与造血干细胞(HSC;控制)进行比较。感谢Alessandro Raveane的代码。
# function to read GMT file[email protected] %>%[email protected] %>% filter(cell_type_singler_blueprintencode_main == 'HSC') %>% nrow()),[email protected] %>% filter(cell_type_singler_blueprintencode_main == 'Monocytes') %>% nrow()) 16 保存结果和Session info 最后我们将seurat对象保存，可以用于后续分析
saveRDS(seurat, 'data/seurat.rds')保存我们使用的R和R包版本信息，以便后续进行可重复性分析。
devtools::session_info()
