单细胞小提琴图自己画

# 前面构建 sce 对象的代码，见往期教程：
markers_df <- FindMarkers(object = sce, ident.1 = 0, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features =markers_genes,log =T )
FeaturePlot(object = sce, features=markers_genes )

#violin plots represent gene expression for each subpopulation. The color of each violin represents the mean gene expression after log2 transformation.
#gene: Gene name of interest as string. DATAuse: Gene expression matrix with rownames containing gene names. tsne.popus = dbscan output, axis= if FALSE no axis is printed. legend_position= default "none" indicates where legend is plotted. gene_name = if FALSE gene name will not be included in the plot.
plot.violin2 <- function(gene, DATAuse, tsne.popus, axis=FALSE, legend_position="none", gene_name=FALSE){
  
  testframe<-data.frame(expression=as.numeric(DATAuse[paste(gene),]), Population=tsne.popus$cluster)
  testframe$Population <- as.factor(testframe$Population)
  colnames(testframe)<-c("expression", "Population")
  
  col.mean<-vector()
  for(i in levels(testframe$Population)){
    col.mean<-c(col.mean,mean(testframe$expression[which(testframe$Population ==i)]))
  }
  col.mean<-log2(col.mean+1)
  
  col.means<-vector()
  
  for(i in testframe$Population){
    col.means<-c(col.means,col.mean[as.numeric(i)])
  }
  testframe$Mean<-col.means
  testframe$expression<-log2(testframe$expression+1)
  
  
  
  p <- ggplot(testframe, aes(x=Population, y=expression, fill= Mean, color=Mean))+
    geom_violin(scale="width") +
    labs(title=paste(gene), y ="log2(expression)", x="Population")+
    theme_classic() +
    
    
    
    
    scale_color_gradientn(colors = c("#FFFF00", "#FFD000","#FF0000","#360101"), limits=c(0,14))+
    scale_fill_gradientn(colors = c("#FFFF00", "#FFD000","#FF0000","#360101"), limits=c(0,14))+
    
    theme(axis.title.y =  element_blank())+
    theme(axis.ticks.y =  element_blank())+
    theme(axis.line.y =   element_blank())+
    theme(axis.text.y =   element_blank())+
    theme(axis.title.x = element_blank())+
    
    
    theme(legend.position=legend_position )
  
  if(axis == FALSE){
    p<-p+
      theme( axis.line.x=element_blank(),
             axis.text.x = element_blank(),
             axis.ticks.x = element_blank())
    
  }
  
  if(gene_name == FALSE){
    p<-p+  theme(plot.title = element_blank())
  }else{ p<-p + theme(plot.title = element_text(size=10,face="bold"))}
  
  p
  
}

plot.violin2(gene = "Pdgfra", DATAuse = RPKM.full, tsne.popus = CAFgroups_full)

# 其中X是稀疏矩阵，拿到指定基因表达量
v=x[,match('CD8A',genes)]
# 然后meta是细胞的表型信息，其中meta$Major.cell.type列是单细胞分群结果
boxplot(v~meta$Major.cell.type,las=2)
library(ggpubr)
df=data.frame(v=v,type=meta$Most.likely.LM22.cell.type)
ggboxplot(df, "type", "v",  color = "type" ) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

v=x[,match('CD8A',genes)]
boxplot(v~meta$Major.cell.type,las=2)
library(ggpubr)
df=data.frame(v=v,type=meta$Most.likely.LM22.cell.type)
cm=unlist(lapply(split(df,df$type), function(x){
  mean(x[,1])
}))
# 主要的代码就是依据表达量均值，来确定细胞亚群的顺序。
df$type=factor(df$type,levels = names(sort(cm,decreasing = T)))
ggboxplot(df, "type", "v",  color = "type" ) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))