ixxmu
commented
4 days ago 单细胞空间转录组中细胞类型与通路关联分析（单细胞丰度与细胞通路相关性分析)-多组差异分析、相关性分析 热图展示 by 生信小博士
上一次我们分享过单细胞多组差异分析-多组火山图
那么如何使用定制数据画多组火山图呢？
准备数据
#request 2.libPaths(c( "/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2","/home/data/t040413/R/yll/usr/local/lib/R/site-library","/refdir/Rlib/", "/usr/local/lib/R/library"))dir.create("~/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types")setwd("~/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types")library(Seurat)library(dplyr)library(scRNAtoolVis)load("~/silicosis/spatial/new_spatial/d.all_with_metadata.rds")d_all=d.allload("~/silicosis/spatial_c2l_results/silicosis_mapped_after_kaiti/barcodes_cellltypes_proportion.Rdata")load("~/silicosis/spatial_c2l_results/silicosis_mapped_after_kaiti/sc_silicosis_analysis/cell2location_map/st_meta_rotate.rds")load("~/silicosis/spatial_c2l_results/silicosis_mapped_after_kaiti/major_abundance.rds")head(major_abundance)head(d.all@meta.data)head(st_meta)head(wide_df)boxplot(wide_df$AT2.cell.1)boxplot(wide_df$AT2.cell.2)#1 获取细胞丰度信息------major_abundance=major_abundance[d.all$col,] ;head(major_abundance)identical( unname(d.all$col) ,rownames(major_abundance))dim(d.all)dim(major_abundance)head(d.all@meta.data)head(major_abundance)delete_dot=function(x ) {x <- gsub(pattern = "\\.+", replacement = " ", x, perl = TRUE) # Apply gsub directly to xreturn(x)}delete_dot("asfdafd.afd")colnames(major_abundance)=delete_dot(colnames(major_abundance)) ;head(major_abundance)head(d.all@meta.data)head(major_abundance)#2获取细胞信号信息----head(d.all@meta.data)#2.1提取并制备人的mouse列表---------if(F){msigdbr::msigdbr_collections()all_gene_sets_hs = msigdbr::msigdbr(species = "Homo sapiens") #Mus musculus Homo sapiens ,category = "C2",subcategory = "CP"dim(all_gene_sets_hs)# saveRDS(all_gene_sets_hs,file="~/datasets/all_gene_sets_hs_msigdb.rds")all_gene_sets_hs = msigdbr::msigdbr(species = "Mus musculus" ) #Mus musculus Homo sapiensdim(all_gene_sets_hs)# saveRDS(all_gene_sets_hs,file="~/datasets/all_gene_sets_mm_msigdb.rds")all_gene_sets_hstable(all_gene_sets_hs$gs_subcat)table(all_gene_sets_hs$gs_cat)#假设我们这里想要寻找的是APOPTOSIS相关通路#pattern参数内输入想要寻找的关键词，这里用的是"APOPTOSIS"h2 <- all_gene_sets_hs #all_gene_sets_hs[grep(pattern = "APOPTOSIS",x = all_gene_sets_hs$gs_name,ignore.case = TRUE),]length(unique(h2$gs_name))#查看唯一通路length(unique(h2$human_gene_symbol))#查看所有通路中的唯一基因length(unique(h2$gene_symbol))#查看所有通路中的唯一基因table(h2$gs_subcat)#table(h2$gs_name)h2getwd()#openxlsx::write.xlsx(h2,"allPathway.xlsx")#保存结果#save(h2,file = "allPathway.rds")# 将数据转换为嵌套列表------------nested_list <- h2 %>%group_by(gs_cat, gs_name) %>%summarise(gene_symbol = list(gene_symbol), .groups = 'drop') %>%group_by(gs_cat) %>%summarise(gs_name = list(set_names(gene_symbol, gs_name)), .groups = 'drop') %>%deframe()# 查看结果# nested_list$C1[[1]]# head(names(nested_list$C1))mygenes=nested_list}#load("~/silicosis/spatial/new_spatial/d.all_with_metadata.rds")DimPlot(d.all,label = TRUE)SpatialDimPlot(d.all,ncol = 2)if (F) {get_one_slideimage=function(slide,image_name="NS_7"){# slide=eachslidetmpimages=slide@images[[image_name]]slide@images=list()slide@images[[image_name]]=tmpimagesreturn(slide)}SiO2_7=get_one_slideimage(d.all,image_name = "SiO2_7") ;SiO2_7table(Idents(SiO2_7)=='AM3⁺ Inmt⁻')head(d.all@meta.data)SiO2_7$leiden_res0.8= as.factor(SiO2_7$leiden_res0.8)is.na(SiO2_7$leiden_res0.8) %>%table()Idents(SiO2_7)=SiO2_7$leiden_res0.8SpatialDimPlot(d.all, #group.by = 'leiden_res0.8',cells.highlight = colnames(SiO2_7)[SiO2_7$leiden_res0.8=='1'],ncol = 2 #, images = "SiO2_7")Idents(sample_seurat) = sample_seurat$opt_clust # opt_clust为seuratclusters.colors=mycol[15:25]names(colors) <- Idents(sample_seurat) %>% levels() ## 命名非常重要spatial_embedding <- SpatialDimPlot( sample_seurat, images = "slice1", group.by = c("opt_clust"), pt.size.factor = 1.15, label = TRUE, label.size = 6, repel = TRUE, combine = FALSE, cols = colors)spatial_embedding ## 成功！}mylist=SplitObject(d.all,split.by = "stim")lapply( mylist, dim)# 含每行非零元素的个数# d.all@assays$Spatial@counts %>% apply(., 2, function(x) sum(x != 0)) #sapply(., function(x) sum(x != 0)) #colSums()d.all$genes_per_spot=d.all@assays$SCT@counts %>% apply(., 2, function(x) sum(x != 0)) #sapply(., function(x) sum(x != 0)) #colSums()d.all@assays$SCT@meta.features$spot_per_gene=d.all@assays$SCT@counts %>% apply(., 1, function(x) sum(x != 0))
计算细胞丰度与细胞通路相关性
#2.2cc包装成函数 -----if (T) {# 导入必要的包library(Seurat)library(dplyr)# 定义函数visualize_pathway <- function(h2_path, d_all, pathway_name,min.cutoff = "q0.8",max.cutoff = "q99.99") {# 读取h2数据# 检查 h2 对象是否存在if (!exists("h2", envir = .GlobalEnv)) {# 如果 h2 对象不存在，则加载它h2 <<- readRDS(h2_path)} else {# 如果 h2 对象已经存在，则打印提示信息print("h2 对象已存在，不需要加载。")} #在这个版本中，h2 对象使用 <<- 操作符在全局环境中定义。如果你不想使用全局变量，可以在函数外部加载 h2 对象，并将其作为参数传递给函数# 打印h2对象# print(h2)# pathway_name="GOCC_ALVEOLAR_LAMELLAR_BODY" #"gocc-ALVEOLAR-LAMELLAR-BODY"# 检查含有特定字符串的通路matched_pathways <- h2[grepl(h2$gs_name, pattern = pathway_name, ignore.case = TRUE),]# 打印匹配的通路print(matched_pathways)# 筛选特定通路desired_list <- h2[h2$gs_name == gsub(pattern = "-", replacement = "_", x = pathway_name),] %>%split(x = .$gene_symbol, f = .$gs_name);print(desired_list)# 添加模块得分d_all <- AddModuleScore(d_all, features = desired_list)# 更新元数据d_all@meta.data[[names(desired_list)]] <- d_all@meta.data$Cluster1head(d_all@meta.data)return(d_all)# 可视化特定通路# SpatialFeaturePlot(d_all,crop = FALSE, keep.scale = 'all', pt.size.factor=1,# features = names(desired_list),# ncol = 2,# # stroke = 0,# max.cutoff =max.cutoff,# min.cutoff =min.cutoff# )## jpeg( paste0(names(desired_list),".jpeg"), res = 300, width = 22, height = 22, units = "in")# p=SpatialFeaturePlot(d_all,crop = FALSE, keep.scale = 'all', pt.size.factor=1,# features = names(desired_list),# ncol = 2 ,# # stroke = 0,# max.cutoff =max.cutoff,# min.cutoff =min.cutoff# )# plot(p)## dev.off()#}}#2.3示例调用----h2_path <- "~/datasets/all_gene_sets_mm_msigdb.rds"d_all <- d.all # 替换为你的d.all对象head(d_all@meta.data)pathway_name <- "ALTEMEIER-RESPONSE-TO-LPS-WITH-MECHANICAL-VENTILATION" #"GOCC-LAMELLAR-BODY"#"REACTOME-RELEASE-OF-APOPTOTIC-FACTORS-FROM-THE-MITOCHONDRIA" #"GOBP-POSITIVE-REGULATION-OF-APOPTOTIC-CELL-CLEARANCE"# "GOBP-NEGATIVE-REGULATION-OF-MITOCHONDRIAL-MEMBRANE-PERMEABILITY-INVOLVED-IN-APOPTOTIC-PROCESS" #"GOBP-POSITIVE-REGULATION-OF-APOPTOTIC-CELL-CLEARANCE" # "GOBP-INTRINSIC-APOPTOTIC-SIGNALING-PATHWAY-IN-RESPONSE-TO-HYPOXIA" #"REACTOME-APOPTOTIC-CLEAVAGE-OF-CELL-ADHESION-PROTEINS"#"GOBP-NEGATIVE-REGULATION-OF-MITOCHONDRIAL-MEMBRANE-PERMEABILITY-INVOLVED-IN-APOPTOTIC-PROCESS" #"GOCC_ALVEOLAR_LAMELLAR_BODY"visualize_pathway(h2_path = h2_path, d_all = d_all,pathway_name = pathway_name, min.cutoff = 'q0.1',max.cutoff = 'q99.9')getwd()head(d_all@meta.data)#2.4 批量产生信号通路信息----markers=read.csv("~/silicosis_spatial/features_in spatial/gsva_markers_Deseq2_difpct0.5.csv")head(markers);dim(markers)mygenes= markers %>% group_by(cluster) %>% slice_max(order_by = avg_log2FC,n = 50) %>%dplyr:: select(gene) %>%.$gene %>%unique() ;mygenesfor (each in mygenes) {d_all=visualize_pathway(h2_path = h2_path, d_all = d_all,pathway_name = each,min.cutoff = "q1",max.cutoff = "q99.9")}load("/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/cell_signal_in_spatial.rds")cell_signals_in_spatialhead(d_all@meta.data)cell_signals_in_spatial=d_all@meta.datarownames(cell_signals_in_spatial)=d_all@meta.data$colgetwd()#save(cell_signals_in_spatial,file = "/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/cell_signal_in_spatial.rds")load("/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/cell_signal_in_spatial.rds")#3.1 celltypes_info----load("~/silicosis/spatial_c2l_results/silicosis_mapped_after_kaiti/major_abundance.rds")head(major_abundance)cell_abundance=major_abundance[d.all$col,] ;head(major_abundance)identical( unname(d.all$col) ,rownames(cell_abundance))delete_dot=function(x ) {x <- gsub(pattern = "\\.+", replacement = " ", x, perl = TRUE) # Apply gsub directly to xreturn(x)}delete_dot("asfdafd.afd")colnames(cell_abundance)=delete_dot(colnames(cell_abundance)) ;head(cell_abundance)head(d_all@meta.data)head(cell_abundance);dim(cell_abundance)#3.2 features_info----features_abundance= cell_signals_in_spatial[ ,which( colnames(cell_signals_in_spatial)=="NAKAMURA_BRONCHIAL_AND_BRONCHIOLAR_EPITHELIA" ) :which( colnames(cell_signals_in_spatial)=="GOBP_REGULATION_OF_CARDIAC_MUSCLE_CELL_MEMBRANE_POTENTIAL")]head(features_abundance);dim(features_abundance)colnames(features_abundance)features_abundance= features_abundance [ unname(d.all$col),] ;head(features_abundance )identical( unname(d.all$col) ,rownames(features_abundance ))head(cell_abundance)[,1:6]head(features_abundance)[,1:6]##验证是否正确-----d.all$'Tissue structure'=ifelse(grepl(pattern = "Bronchial",x = d.all$clusters),"Airway",ifelse( grepl(pattern = "ascula",x = d.all$clusters) ,"Vessel","Alveolar") )print( identical(rownames(features_abundance),rownames(cell_abundance)))#j计算相关性，并添加adj.pval-------celltype_feature_list = list()for (my_desired_celltype in colnames(cell_abundance)) {combined.data = cbind(cell_abundance[, my_desired_celltype], features_abundance)colnames(combined.data)[1] = my_desired_celltypecombined.data$group = stringr::str_split(string = rownames(combined.data),pattern = "_", simplify = TRUE)[, c(1, 2)] %>%apply(., 1, function(x) paste(x[1], x[2], sep = "_"))combined.data$clusters = d.all$clusterscombined.data$`Tissue structure` = d.all$`Tissue structure`# 数据筛选data_filtered <- combined.data %>%filter(grepl(pattern = "_", x = .$group))# 将数据从宽格式转换为长格式data_long <- data_filtered %>%tidyr::pivot_longer(cols = 2:237, names_to = "features", values_to = "Value")# 计算相关系数和p值correlations <- data_long %>%group_by(features) %>%summarise(Correlation = cor.test(!!sym(my_desired_celltype), Value)$estimate,P.value = cor.test(!!sym(my_desired_celltype), Value)$p.value) %>%mutate(Significance = case_when(P.value < 0.001 ~ "***",P.value < 0.01 ~ "**",P.value < 0.05 ~ "*",TRUE ~ "ns"),Label = paste(features, "\nr =", round(Correlation, 2), Significance, sep = " "))# 添加调整后的p值correlations <- correlations %>%mutate(adj.pval = p.adjust(P.value, method = "BH")) # 使用Benjamini-Hochberg校正# 添加celltype列correlations$celltype = my_desired_celltype# 存储结果celltype_feature_list[[my_desired_celltype]] = correlationsprint(paste0("done === ", my_desired_celltype))}celltype_feature_listif (T) {# Function to extract top 20 positive and top 20 negative correlationsextract_top_correlations <- function(cor_tibble, n = 20) {# Sort by correlation value in descending order to get top positive correlationstop_pos <- cor_tibble %>% filter(P.value<0.05 ) %>%arrange(desc(Correlation)) %>%slice_head(n = n)# Sort by correlation value in ascending order to get top negative correlationstop_neg <- cor_tibble %>% filter(P.value<0.05 ) %>%arrange(Correlation) %>%slice_head(n = n)# Combine the top positive and negative correlationsresult <- list(top_positive = top_pos,top_negative = top_neg)return(result)}# Apply the function to each cell type tibble in the listtop_correlations_list <- lapply(celltype_feature_list, extract_top_correlations)# Print the results for the first cell type as an exampletop_correlations_list[[1]]}celltype_feature_listcelltype_feature_df= Reduce( rbind,celltype_feature_list) #%>%head()View(celltype_feature_df)library(dplyr)#给这个数据框添加新的一列Correlation_type-----if (F) {celltype_feature_df$avg_log2FC=celltype_feature_df$Correlationcelltype_feature_df$p_val_adj=celltype_feature_df$adj.pvalcelltype_feature_df$cluster=celltype_feature_df$celltypecelltype_feature_df$gene=celltype_feature_df$featurescelltype_feature_df$p_val=celltype_feature_df$P.valuecelltype_feature_df$pct.2=celltype_feature_df$p_val_adjcelltype_feature_df$pct.1=0celltype_feature_df$log_adj.pval=log ( celltype_feature_df$p_val_adj +0.00000000001)celltype_feature_df$pct.2=celltype_feature_df$log_adj.pval}celltype_feature_df <- celltype_feature_df %>%group_by(cluster) %>%mutate(# Define quantiles for top and bottom 5%top5_threshold = quantile(Correlation, 0.95),bottom5_threshold = quantile(Correlation, 0.05),# Create Correlation_type based on conditionsCorrelation_type = case_when(Correlation >= top5_threshold ~ "top5%",Correlation <= bottom5_threshold ~ "bottom5%",TRUE ~ "middle"))# Remove temporary columns used for threshold calculations# %>% select(-top5_threshold, -bottom5_threshold)celltype_feature_df <- celltype_feature_df %>%group_by(cluster) %>%mutate(# Define quantiles for top and bottom 1%top1_threshold = quantile(Correlation, 0.99),bottom1_threshold = quantile(Correlation, 0.01),# Create Correlation_type based on conditionsCorrelation_type2 = case_when(Correlation >= top1_threshold ~ "top1%",Correlation <= bottom1_threshold ~ "bottom1%",TRUE ~ "middle"))# View the updated dataframeprint(celltype_feature_df)# 定义函数#save(celltype_feature_list,celltype_feature_df,top_correlations_list,file = "/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/celltype_features_correlations.rdata")load("/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/celltype_features_correlations.rdata")
选择想要的细胞通路与细胞类型
h2_path="~/datasets/all_gene_sets_mm_msigdb.rds"# 检查 h2 对象是否存在if (!exists("h2", envir = .GlobalEnv)) {# 如果 h2 对象不存在，则加载它h2 <<- readRDS(h2_path)} else {# 如果 h2 对象已经存在，则打印提示信息print("h2 对象已存在，不需要加载。")} #在这个版本中，h2 对象使用 <<- 操作符在全局环境中定义。如果你不想使用全局变量，可以在函数外部加载 h2 对象，并将其作为参数传递给函数h2_list= split( x = h2[,c('gene_symbol' ,'gs_name')] ,f = h2$gs_name)grep (pattern = "BP",x = h2$gs_subcat,ignore.case = T,value = TRUE)intersect(h2$gs_name,celltype_feature_df$features)intersect(h2$gs_name,names(h2_list))celltype_feature_df$genes# for (each in celltype_feature_df$features ) {# # each="AUJLA_IL22_AND_IL17A_SIGNALING"## if ( each %in% names(h2_list) ) {# celltype_feature_df [ celltype_feature_df$features==each,]$genes = list(h2_list[[each]]$gene_symbol)## }## }# 确保 genes 列存在，如果没有，可以初始化一个空列if (!"genes" %in% colnames(celltype_feature_df)) {celltype_feature_df$genes <- NA}for (each in celltype_feature_df$features) {# 如果 'features' 在 h2_list 中if (each %in% names(h2_list)) {# 使用 paste 函数将基因符号合并成一个字符串gene_symbols <- paste(h2_list[[each]]$gene_symbol, collapse = ", ")# 找到相应的行，并将基因符号字符串赋值到 'genes' 列celltype_feature_df$genes[celltype_feature_df$features == each] <- gene_symbols}}celltype_feature_dfif (T) {library(dplyr)library(tidyr)# 将每一行的基因名称按逗号分隔开，创建一个新的行来表示每个基因df_genes <- celltype_feature_df %>% ungroup( ) %>% filter(adj.pval<0.01 &grepl(x = .$features,pattern="blood|vess|endotheli",ignore.case=T)) %>%dplyr::select(features,genes) %>%# 使用 separate_rows 来分割 genes 列separate_rows(genes, sep = ", ") %>%# 去除前后空格（确保一致性）mutate(genes = trimws(genes))# 统计每个基因的出现频率gene_frequency <- df_genes %>%group_by(genes) %>%summarise(frequency = n()) %>%arrange(desc(frequency)) ;gene_frequency # 按频率降序排列# 输出结果print(gene_frequency)gene_frequency$genes}
计算基因出现频率
画出所有细胞通路与细胞之间的相关性分析图
celltype_feature_listcelltype_feature_selected= celltype_feature_df %>%group_by(celltype) %>%filter( P.value<0.01 & abs(Correlation)>0.2 )table(celltype_feature_selected$celltype)length(unique(celltype_feature_selected$celltype))## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster genehead(celltype_feature_df)library(scRNAtoolVis)# testdata('pbmc.markers')head(pbmc.markers)# plotjjVolcano(diffData = pbmc.markers )celltype_feature_df$avg_log2FC=celltype_feature_df$Correlationcelltype_feature_df$p_val_adj=celltype_feature_df$adj.pvalcelltype_feature_df$cluster=celltype_feature_df$celltypecelltype_feature_df$gene=celltype_feature_df$featurescelltype_feature_df$p_val=celltype_feature_df$P.valuecelltype_feature_df$pct.2=celltype_feature_df$p_val_adjcelltype_feature_df$pct.1=0celltype_feature_df$log_adj.pval=log ( celltype_feature_df$p_val_adj +0.000000000000000001)celltype_feature_df$pct.2=celltype_feature_df$log_adj.pval#write.csv(celltype_feature_df,file = "/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/celltype_feature_df.csv")getwd()head(celltype_feature_df)# supply own genesmygene <- c('NAKAMURA_BRONCHIAL_AND_BRONCHIOLAR_EPITHELIA',"GOBP_SEQUESTERING_OF_IRON_ION","GOBP_LYMPHATIC_ENDOTHELIAL_CELL_DIFFERENTIATION")trace(markerVocalno, edit = T)#修改函数------markerVocalno <- function (markers = NULL, ownGene = NULL, topn = 5, log2FC = 0.25, toppos=NULL, topnegtive=NULL, lableSize=3,labelCol = NULL, hlineSize = 1, hlineColor = "grey50", random_number=0.001,pforce = 5, nforce = 2.5, nudge_x = 0.8, pnudge_y = 0.25,nnudge_y = 0, base_size = 14, facetColor = NA, facetFill = "white",ylab = "Log2-Fold Change", nrow = 1) {# 筛选顶层正负基因if (is.null(ownGene)) {#此处可以使用自己筛选好想要添加label的基因----toppos <- toppos %>% dplyr::group_by(cluster) %>% dplyr::slice_max(avg_log2FC, n = topn) #markers %>% dplyr::group_by(cluster) %>% dplyr::slice_max(avg_log2FC, n = topn)topnegtive <- topnegtive %>% dplyr::group_by(cluster) %>% dplyr::slice_min(avg_log2FC, n = topn) #markers %>% dplyr::group_by(cluster) %>% dplyr::slice_min(avg_log2FC, n = topn)topgene <- rbind(toppos, topnegtive)} else {topgene <- markers %>% dplyr::filter(gene %in% ownGene)toppos <- topgene %>% dplyr::filter(avg_log2FC > 0)topnegtive <- topgene %>% dplyr::filter(avg_log2FC < 0)}# 生成基础散点图p1 <- ggplot2::ggplot(markers, ggplot2::aes(x = pct.1 - pct.2, y = avg_log2FC)) +ggplot2::geom_point(color = "grey80") +ggplot2::geom_hline(yintercept = c(-log2FC, log2FC), lty = "dashed",size = hlineSize, color = hlineColor)# 加入文本标签和颜色set.seed(1)p2 <- p1 +ggrepel::geom_text_repel(data = toppos, ggplot2::aes(x = pct.1 - pct.2, y = avg_log2FC,label = gene, color = cluster),show.legend = F, direction = "y", hjust = 1,nudge_y = pnudge_y,#调整label大小 增加随机波动size = lableSize, force = pforce,nudge_x = -nudge_x - (toppos$pct.1 - toppos$pct.2) +runif(1,min = 0,max=random_number) ) +ggrepel::geom_text_repel(data = topnegtive, ggplot2::aes(x = pct.1 - pct.2,y = avg_log2FC, label = gene,color = cluster),show.legend = F, direction = "y", hjust = 0,nudge_y = nnudge_y,#调整label大小 增加随机波动size = lableSize, force = nforce,nudge_x = nudge_x - (topnegtive$pct.1 - topnegtive$pct.2) +runif(1,min = 0,max=random_number) ) +ggplot2::geom_point(data = topgene, ggplot2::aes(x = pct.1 - pct.2, y = avg_log2FC,color = cluster), show.legend = F) +ggplot2::scale_color_manual(name = "", values = labelCol) +ggplot2::theme_bw(base_size = base_size) +ggplot2::theme(panel.grid = ggplot2::element_blank(),axis.text.x = ggplot2::element_text(angle = 45, hjust = 1),strip.background = ggplot2::element_rect(color = facetColor, fill = facetFill)) +ggplot2::xlab(expression(Delta ~ "Percentage Difference")) +ggplot2::ylab(ylab) +ggplot2::facet_wrap(~cluster, nrow = nrow, scales = "fixed")return(p2)}#top5-----p= markerVocalno(markers = celltype_feature_df, lableSize = 1.62, #ownGene = mygene,# topn = 5,toppos = celltype_feature_df[ celltype_feature_df$Correlation_type =='top5%',],topnegtive = celltype_feature_df[ celltype_feature_df$Correlation_type =='bottom5%',],pforce = 2,nforce = 0.9,random_number =0.001,log2FC = 0.2, nrow = 5,labelCol = c("#EE3B3B", "#FF7F00", "#CD6600", "#8B2323", "darkred","darkgreen", "#008B8B", "#FFB90F", "#1F1F1F", "#66CD00","#0000FF", "darkorange", "#528B8B", "#9400D3", "#EE1289","#00BFFF", "#00FF00", "#191970", "darkblue", "#4A708B", "#00FF7F", "#8B8B00","#FF1493", "#FFA500", "#8B4513","red","purple","lightblue","#F5F5DC", "#F0FFFF"))+ylab("Correlation")+xlab("-Log adj.Pval")ggsave(plot=p,filename = "/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/top1——top5-pnfore-labelsize1-random.6234.jpeg",width = 20,height = 20,limitsize = FALSE,dpi = 300)getwd()ggsave(plot=p,filename = "/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/top1——top5-pnfore-labelsize1-random.6234.pdf",width = 20,height = 20,limitsize = FALSE,dpi = 600)getwd()runif(1,0,0).libPaths(c( # "/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2",#"/home/data/t040413/R/yll/usr/local/lib/R/site-library","/refdir/Rlib/", "/usr/local/lib/R/library"))library(export)graph2ppt(x=p,file="celltypes_signals_in spatial.pptx",width = 18,height = 18)#top1-----markerVocalno(markers = celltype_feature_df,lableSize = 2, #ownGene = mygene,topn = 5,toppos = celltype_feature_df[ celltype_feature_df$Correlation_type2=='top1%',],topnegtive = celltype_feature_df[ celltype_feature_df$Correlation_type2=='bottom1%',],log2FC = 0.2, nrow = 5,labelCol = c("#EE3B3B", "#FF7F00", "#CD6600", "#8B2323", "darkred","darkgreen", "#008B8B", "#FFB90F", "#1F1F1F", "#66CD00","#0000FF", "darkorange", "#528B8B", "#9400D3", "#EE1289","#00BFFF", "#00FF00", "#191970", "darkblue", "#4A708B", "#00FF7F", "#8B8B00","#FF1493", "#FFA500", "#8B4513","red","purple","lightblue","#F5F5DC", "#F0FFFF"))+ylab("Correlation")+xlab("-Log adj.Pval")ggsave(filename = "/home/data/t040413/silicosis/spatial/new_spatial/20241011-cell_signals_in spatial cell types/top1——top1-=.png",width = 18,height = 18,limitsize = FALSE,dpi = 900)getwd()#trance the library is loaded, usetrace(jjVolcano, edit = T)getAnywhere(jjVolcano)edit(jjVolcano)pkg_path <- find.package("scRNAtoolVis")print(pkg_path)list.files(file.path(pkg_path, "R"), full.names = TRUE)file.edit("path/to/your/file.R")# supply own genesmygene <- c('NAKAMURA_BRONCHIAL_AND_BRONCHIOLAR_EPITHELIA', "GOBP_SEQUESTERING_OF_IRON_ION","GOBP_LYMPHATIC_ENDOTHELIAL_CELL_DIFFERENTIATION")markerVocalno(markers = celltype_feature_df,labelCol = seq(1,25,1))jjVolcano(diffData = celltype_feature_df,order.by = "avg_log2FC",aesCol = c('blue','red'),col.type ='adjustP', # "updown",# log2FC.cutoff=1,# topGeneN = 3,myMarkers = mygene,# tile.col = corrplot::COL2('RdBu', 15)[4:12] ,fontface = 'italic',cluster.order = rev(unique(celltype_feature_df$cluster)),tile.col = c("#EE3B3B", "#FF7F00", "#CD6600", "#8B2323", "#DEB887", "#76EEC6","#F0FFFF", "#008B8B", "#FFB90F", "#F5F5DC", "#1F1F1F", "#66CD00","#0000FF", "#97FFFF", "#528B8B", "#9400D3", "#EE1289", "#00BFFF","#00FF00", "#191970", "#FFFF00", "#4A708B", "#00FF7F", "#8B8B00","#FF1493", "#FFA500", "#8B4513","red","grey","lightblue"),celltypeSize=3,base_size = 9) +ggplot2::xlab("Clusters") + ggplot2::ylab("Correlation") +ggplot2::scale_y_continuous(n.breaks = 6)jjVolcano(diffData = celltype_feature_df,polar = T,order.by = "avg_log2FC",aesCol = c('blue','red'),# log2FC.cutoff=1,# topGeneN = 3,myMarkers = mygene,# tile.col = corrplot::COL2('RdBu', 15)[4:12] ,fontface = 'italic', fontsize=0.1,cluster.order = rev(unique(celltype_feature_df$cluster)),tile.col = c("#EE3B3B", "#FF7F00", "#CD6600", "#8B2323", "#DEB887", "#76EEC6","#F0FFFF", "#008B8B", "#FFB90F", "#F5F5DC", "#1F1F1F", "#66CD00","#0000FF", "#97FFFF", "#528B8B", "#9400D3", "#EE1289", "#00BFFF","#00FF00", "#191970", "#FFFF00", "#4A708B", "#00FF7F", "#8B8B00","#FF1493", "#FFA500", "#8B4513","red","grey","lightblue"))
https://mp.weixin.qq.com/s/P3U2krL8D52IQZjHqev9-A