上一次我们分享过单细胞多组差异分析-多组火山图
那么如何使用定制数据画多组火山图呢?
准备数据
#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.all
load("~/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 x
return(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 sapiens
dim(all_gene_sets_hs)
# saveRDS(all_gene_sets_hs,file="~/datasets/all_gene_sets_mm_msigdb.rds")
all_gene_sets_hs
table(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)
h2
getwd()
#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=eachslide
tmpimages=slide@images[[image_name]]
slide@images=list()
slide@images[[image_name]]=tmpimages
return(slide)
}
SiO2_7=get_one_slideimage(d.all,image_name = "SiO2_7") ;SiO2_7
table(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.8
SpatialDimPlot(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为seurat
clusters.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$Cluster1
head(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() ;mygenes
for (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_spatial
head(d_all@meta.data)
cell_signals_in_spatial=d_all@meta.data
rownames(cell_signals_in_spatial)=d_all@meta.data$col
getwd()
#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 x
return(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_celltype
combined.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$clusters
combined.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]] = correlations
print(paste0("done === ", my_desired_celltype))
}
celltype_feature_list
if (T) {
# Function to extract top 20 positive and top 20 negative correlations
extract_top_correlations <- function(cor_tibble, n = 20) {
# Sort by correlation value in descending order to get top positive correlations
top_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 correlations
top_neg <- cor_tibble %>% filter(P.value<0.05 ) %>%
arrange(Correlation) %>%
slice_head(n = n)
# Combine the top positive and negative correlations
result <- list(
top_positive = top_pos,
top_negative = top_neg
)
return(result)
}
# Apply the function to each cell type tibble in the list
top_correlations_list <- lapply(celltype_feature_list, extract_top_correlations)
# Print the results for the first cell type as an example
top_correlations_list[[1]]
}
celltype_feature_list
celltype_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$Correlation
celltype_feature_df$p_val_adj=celltype_feature_df$adj.pval
celltype_feature_df$cluster=celltype_feature_df$celltype
celltype_feature_df$gene=celltype_feature_df$features
celltype_feature_df$p_val=celltype_feature_df$P.value
celltype_feature_df$pct.2=celltype_feature_df$p_val_adj
celltype_feature_df$pct.1=0
celltype_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 conditions
Correlation_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 conditions
Correlation_type2 = case_when(
Correlation >= top1_threshold ~ "top1%",
Correlation <= bottom1_threshold ~ "bottom1%",
TRUE ~ "middle"
)
)
# View the updated dataframe
print(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_df
if (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_list
celltype_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 gene
head(celltype_feature_df)
library(scRNAtoolVis)
# test
data('pbmc.markers')
head(pbmc.markers)
# plot
jjVolcano(diffData = pbmc.markers )
celltype_feature_df$avg_log2FC=celltype_feature_df$Correlation
celltype_feature_df$p_val_adj=celltype_feature_df$adj.pval
celltype_feature_df$cluster=celltype_feature_df$celltype
celltype_feature_df$gene=celltype_feature_df$features
celltype_feature_df$p_val=celltype_feature_df$P.value
celltype_feature_df$pct.2=celltype_feature_df$p_val_adj
celltype_feature_df$pct.1=0
celltype_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 genes
mygene <- 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, use
trace(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 genes
mygene <- 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")
)