ggtree 绘制 scRNA-seq 细胞聚类圈图

图片参考：郭老师的nc

点击查看代码

library(Matrix)
library(Matrix.utils)
library(plyr)
library(dplyr)
library(Seurat)
library(sctransform)
library(igraph)
library(factoextra)
library(ComplexHeatmap)
library(circlize)
library(EpicTools)
require(Hmisc)
require(dplyr)
require(openxlsx)
require(ggplot2)
library(ggpubr)
require(cowplot)
library(data.table)
library(RColorBrewer)
library(rowr)
library(SingleR)
library(scater)
library(pheatmap)
#library(nichenetr)
library(tidyverse)

path = "./"
cm.list = paste0(path, list.files(pattern = "*.matrices.rds", path = path))
cm.files <- lapply(cm.list, readRDS)

names(cm.files) <- sub(path,"",
                       sub("\\_cell.counts.matrices.rds", "", cm.list))

options(repr.plot.height = 4,repr.plot.width = 4)

cm.pp <- mapply(EpicPreHS, cm.files, orig.ident = names(cm.files), SIMPLIFY = F)

save(cm.pp,file = 'cm.pp.rda')

cm.pp1 = cm.pp[['560']]

covid_combined.emat = cm.pp1[['emat']]
covid_combined.emat = as.matrix(covid_combined.emat)

covid_combined.nmat = cm.pp1[['nmat']]
covid_combined.nmat = as.matrix(covid_combined.nmat)

set.seed(1245)

covid_combined <- CreateSeuratObject(counts = covid_combined.emat, min.cells = 10, names.field = 1, names.delim = "\\.")
covid_combined <- PercentageFeatureSet(covid_combined, pattern = "^MT-", col.name = "percent.mt")
covid_combined <- PercentageFeatureSet(covid_combined, pattern = "^RPS", col.name = "percent.rps")
covid_combined <- PercentageFeatureSet(covid_combined, pattern = "^RPL", col.name = "percent.rpl")
covid_combined <- PercentageFeatureSet(covid_combined, pattern = "^RNA\\d8S5", col.name = "percent.rrna")
covid_combined <- SCTransform(covid_combined, vars.to.regress = c("percent.mt", "percent.rps", "percent.rpl", "percent.rrna", "nCount_RNA", "nFeature_RNA"), verbose = FALSE, return.only.var.genes = TRUE) #expect "iteration limit reached" warning unless suppressed per https://github.com/satijalab/seurat/issues/1426
covid_combined <- RunPCA(covid_combined, verbose = FALSE)
covid_combined <- RunUMAP(covid_combined, dims = 1:50, verbose = FALSE)
covid_combined <- FindNeighbors(covid_combined, dims = 1:50, verbose = FALSE)
covid_combined <- FindClusters(covid_combined, resolution = 1, verbose = FALSE)
DimPlot(covid_combined, label = TRUE) + NoLegend()

new_dotplot <- function(object = NULL, features = NULL, group.by = NULL, genes.on.x = TRUE, 
                        size.breaks.values = NULL, color.breaks.values = c(-3, -2, -1, 0, 1, 2, 3), shape.scale = 12, 
                        dend_x_var = "Average expression", dend_y_var = "Average expression",
                        cols.use = c("lightgrey", "blue"), scale.by = "radius", col.min = -2.5, col.max = 2.5,
                        dot.min = 0) {
  scale.func <- switch(EXPR = scale.by, size = scale_size, 
                       radius = scale_radius, stop("'scale.by' must be either 'size' or 'radius'"))
  
  data.features <- FetchData(object = object, vars = features)
  object[[group.by, drop = TRUE]]
  data.features$id <- object[[group.by, drop = TRUE]]
  if (!is.factor(x = data.features$id)) {
    data.features$id <- factor(x = data.features$id)
  }
  id.levels <- levels(x = data.features$id)
  data.features$id <- as.vector(x = data.features$id)
  data.plot <- lapply(X = unique(x = data.features$id), FUN = function(ident) {
    data.use <- data.features[data.features$id == ident, 
                              1:(ncol(x = data.features) - 1), drop = FALSE]
    avg.exp <- apply(X = data.use, MARGIN = 2, FUN = function(x) {
      return(mean(x = expm1(x = x)))
    })
    pct.exp <- apply(X = data.use, MARGIN = 2, FUN = Seurat:::PercentAbove, 
                     threshold = 0)
    return(list(avg.exp = avg.exp, pct.exp = pct.exp))
  })
  names(x = data.plot) <- unique(x = data.features$id)
  data.plot <- lapply(X = names(x = data.plot), FUN = function(x) {
    data.use <- as.data.frame(x = data.plot[[x]])
    data.use$features.plot <- rownames(x = data.use)
    data.use$id <- x
    return(data.use)
  })
  
  data.plot <- do.call(what = "rbind", args = data.plot)
  if (!is.null(x = id.levels)) {
    data.plot$id <- factor(x = data.plot$id, levels = id.levels)
  }
  avg.exp.scaled <- sapply(X = unique(x = data.plot$features.plot), 
                           FUN = function(x) {
                             data.use <- data.plot[data.plot$features.plot == 
                                                     x, "avg.exp"]
                             data.use <- scale(x = data.use)
                             data.use <- MinMax(data = data.use, min = col.min, 
                                                max = col.max)
                             return(data.use)
                           })
  avg.exp.scaled <- as.vector(x = t(x = avg.exp.scaled))
  data.plot$avg.exp.scaled <- avg.exp.scaled
  data.plot$features.plot <- factor(x = data.plot$features.plot, 
                                    levels = rev(x = features))
  data.plot$pct.exp[data.plot$pct.exp < dot.min] <- NA
  data.plot$pct.exp <- data.plot$pct.exp * 100
  
  if(genes.on.x){
    data.final <- data.frame(features.plot = data.plot$features.plot, id = data.plot$id)
    data.final <- cbind(data.final, data.plot[,c(2,5)])
  }
  else {
    data.final <- data.frame(id = data.plot$id, features.plot = data.plot$features.plot)
    data.final <- cbind(data.final, data.plot[,c(2,5)])
  }
  colnames(data.final)[3:4] <- c("Percent expressed", "Average expression")
  save(data.final,file = 'res.rda')
}

#source('func.R')

covid_combined

new_dotplot(covid_combined, features = grep("^HLA-", rownames(covid_combined.emat), value = T), group.by = "seurat_clusters", shape.scale = 6, color.breaks.values = c(-2, -1, 0, 1, 2), size.breaks.values = c(0, 25, 50, 75, 100))

save(covid_combined,file = 'covid_combined.rda')

load('covid_combined.rda')

load('res.rda')

data.final1 = data.final %>% dplyr::filter(id != 12)

dim(data.final1)

library(ggplot2)
ggplot(data.final1,aes(x=features.plot,y=id))+
  geom_point(aes(size=`Percent expressed`,
                 color=`Average expression`))+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text.x=element_text(angle=90,hjust = 1,vjust=0.5))+
  scale_color_gradient(low="lightgrey",high="blue")+
  labs(x=NULL,y=NULL)+
  guides(size=guide_legend(order=3))

df<-data.final1[,c(1,2,4)]

df1<-reshape2::dcast(df,id~features.plot,value.var = "Average expression")
rownames(df1)<-df1$id

df1.1<-df1[,2:21]

df1.1.clust<-hclust(dist(df1.1))
library(ggtree)
p2<-ggtree(df1.1.clust)
p2+
  geom_tiplab()+
  xlim(NA,7)

ggtree(df1.1.clust,layout = "circular")+xlim(0,7)+ geom_tiplab2(offset=0.1,size =5)
+ geom_highlight(node = 5,fill="red")

data = fortify(df1.1.clust)
dim(data)

node = c(0,1,2,3,4,5,6,7,8,9,10,11)
group = c('a','a','b','c','a','a','b','c','a','a','b','b')
group.info = data.frame(node,group)
#rownames(group.info) = group.info[,1];group.info = group.info[,-1];group.info = as.data.frame(group.info)

library(ape)

tree = as.phylo(df1.1.clust)# 聚类结果转成系统发育格式
write.tree(tree) # 输出newick格式文件

groupInfo <- split(group.info$node, group.info$group)

# 将分组信息添加到树中
tree <- groupOTU(tree, groupInfo)

options(repr.plot.height = 5, repr.plot.width = 5)

ggtree(tree,branch.length = "none",layout = "circular") + xlim(0,7) + 
    geom_tiplab2(offset=0.1,size =5)+ 
    geom_highlight(node = c(7,9,11),fill="red")  +
 geom_highlight(node = c(5,6,10,3),fill="blue" )+
 geom_highlight(node = c(2,12,1,8,4),fill   ="green") +
    theme(legend.position = "right") + 
    geom_nodepoint(aes(color=group), size=3) +
  geom_strip(taxa1 = "6",
             taxa2 = "8",
             offset = 2,
             barsize = 20,
             extend = 0.5,
             color="#7c90c8",
             label = "YN",
             offset.text = 3)
# 使用geom_point 会有0 

p = ggtree(tree,
       branch.length = "none",
       layout = "circular")+
  geom_tiplab(size=3)+
  geom_nodepoint(aes(color=group))+
  geom_strip(taxa1 = "6",
             taxa2 = "5",
             offset = 1.5,
             barsize = 12,
             extend = 0.5,
             color="#7c90c8")+
  geom_tiplab(size=3,align = TRUE,offset = 1)+
  geom_strip(taxa1 = "0",
             taxa2 = "10",
             offset = 1.5,
             barsize = 12,
             extend = 0.5,
             color="green")+
  geom_tiplab(size=3,align = TRUE,offset = 1)+
  geom_strip(taxa1 = "9",
             taxa2 = "11",
             offset = 1.5,
             barsize = 12,
             extend = 0.5,
             color="red",
             label = "linanlad",
             offset.text = 7)+
  geom_tiplab(size=3,align = TRUE,offset = 1)

library(ggimage)
library(ggtree)