单细胞单细胞单细胞测序技术

单细胞数据处理小细节汇总

2021-05-11  本文已影响0人  Hayley笔记

1. Seurat对象查看当前的Assay

pbmc <- Read10X('./filtered_gene_bc_matrices/hg19/')
pbmc <- CreateSeuratObject(pbmc,project = 'pbmc3k',min.cells = 3,min.features = 200)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc) %>% FindVariableFeatures(nfeatures = 2000) %>% ScaleData(vars.to.regress = "percent.mt")
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt")
DefaultAssay(pbmc)
#[1] "SCT"

在进行了SCTransform操作后，矩阵默认会变成SCT矩阵，如果不加设置，后续的PCA等操作都是基于SCT矩阵。
修改DefaultAssay：

DefaultAssay(pbmc) <- 'RNA'
DefaultAssay(pbmc) 
#[1] "RNA"

2. Seurat使用FindVariables找到高变基因后增删高变基因

var.gene <- pbmc@assays$RNA@var.features #提取高变基因
pbmc@assays$RNA@var.features <- c(var.gene,c('CD4','CD8'))%>%unique() #增加高变基因
#unique()是为了保证增加的基因不和原有的重复

3. 不同运行步骤中的文件来源和储存位置⚠️

⚠️：PCA的值是可以被覆盖的，使用三步法对矩阵进行标准化后进行PCA后再使用SCT矩阵进行标准化，PCA的矩阵变成了SCT的PCA矩阵，原有的PCA矩阵不会保留。后续的TSNE和UMAP降维图也和三步法不一样。

FindVariableFeatures()就是在细胞与细胞间进行比较，选择表达量差别最大的基因，默认情况下，会返回2,000个高可变基因。高可变基因的数目用于做PCA。
Q：PCA输入的是scale.data的数据，scale.data默认做的是FindVariableFeatures()找到的高变基因的scale，但是scale.data如果scale了全部基因，做PCA也就是用的全部基因，那FindVariableFeatures()还有什么用？
A：虽然很多教程都提示内存足够的话建议scale所有的基因，但PCA的时候还是要选择高变基因，否则会引入噪声。低丰度，变化低的基因，很多都是噪音，如果用全部基因去跑PCA，这种低丰度基因的噪音是没有办法通过选择pc数去除的。仔细看Seurat官网，虽然它推荐scale全部基因，但是做PCA的时候还是限定了使用高变基因。

https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

此外如果使用了SCTransform，默认是得到3000个高变基因，Scale的也是高变基因。再下一步直接做RunPCA()，也是默认使用高变基因。

4. 关于有些细胞属于同一个cluster但是在umap或者tsne图上相聚较远的问题：UMAP和TSNE是各自的算法在PCA降维的基础上再进行非线形降维，在二维图上把其各自算法认为相近的细胞聚在一起。但是FindClusters输入的不是UMAP或TSNE降维的数据，而是FindNeighbors的数据，而FindNeighbors输入的数据是PCA降维数据，是用另外一种算法计算的细胞之间的距离。因此会出现有些细胞被认为是同一个cluster，但是在umap或者tsne图上相聚较远（尤其是一些散在的，脱离主群的细胞）

5. marker基因鉴定，查看marker基因的表达是使用RNA矩阵还是sct矩阵？
这是一个争议性问题，两个都可以，目前建议最好使用RNA矩阵。
sct的到的count并不是真实的基因表达值，而是通过scaledata倒推出来的，它是一个回归，运算之后的残差。

6. 关于FindAllMarks找到的基因
如下图，先看cluster0的Marker基因：cluster0的差异基因是cluster0的细胞和剩下的所有的cluster合在一起的细胞做对比得到的。pct.1是这个基因在cluster0中的表达比例（S100A8在cluster0的细胞中的表达比例是100%），pct.2是这个基因在除了cluster0以外的所有细胞中的表达比例（S100A8在除cluster0以外的细胞中的表达比例是51.2%）。avg_log2FC是表达差异倍数，p_val_adj是校正后的p值。

7. 在提取Marker基因时比较好的办法：因为单细胞矩阵算出来的结果，p_val_adj==0的有很多，所以可以先把p_val_adj==0先提出来。再把p_val_adj<0.01的按差异倍数靠前的20/30/500...（按需要）个基因提出来，然后把两个矩阵合在一起（取交集）用来做细胞鉴定。（结合p值和fc来做筛选效果更好）

top = 30   #可根据需要调整
TopMarkers1 <- ClusterMarker %>% filter(p_val_adj == 0) %>% group_by(cluster) %>% 
  top_n(n = top, wt = avg_log2FC)
TopMarkers2 <- ClusterMarker %>% filter(p_val_adj < 0.01) %>% group_by(cluster) %>%
  top_n(n = top, wt = avg_log2FC)
TopMarkers <- rbind(TopMarkers1, TopMarkers2) %>% unique.matrix() %>% arrange(cluster)
write.csv(TopMarkers,'CellType/TopMarkers.csv', row.names=F)

⚠️：提取没有核糖体和线粒体的marker基因更好。（这些基因对鉴定没有帮助）

人的核糖体基因是RPS/RPL开头
人的线粒体基因是MT-开头
小鼠的核糖体基因是Rps/Rpl开头
小鼠的线粒体基因是mt-开头

##人
### 提取没有线粒体、核糖体的Markers
ClusterMarker_noRibo <- ClusterMarker[!grepl("^RP[SL]", ClusterMarker$gene, ignore.case = F),]
ClusterMarker_noRibo_noMito <- ClusterMarker_noRibo[!grepl("^MT-", ClusterMarker_noRibo$gene, ignore.case = F),]
top = 100   #可根据需要调整
TopMarkers1 <- ClusterMarker_noRibo_noMito %>% filter(p_val_adj == 0) %>% 
  group_by(cluster) %>% top_n(n = top, wt = avg_log2FC)
TopMarkers2 <- ClusterMarker_noRibo_noMito %>% filter(p_val_adj < 0.01) %>% 
  group_by(cluster) %>% top_n(n = top, wt = avg_log2FC)
ClusterMarker_noRibo_noMito <- rbind(TopMarkers1, TopMarkers2) %>% unique.matrix() %>% arrange(cluster)
write.csv(ClusterMarker_noRibo_noMito,'CellType/TopMarkers_noRibo_noMito.csv', row.names=F)

##小鼠
### 提取没有线粒体、核糖体的Markers
ClusterMarker_noRibo <- ClusterMarker[!grepl("^Rp[sl]", ClusterMarker$gene, ignore.case = F),]
ClusterMarker_noRibo_noMito <- ClusterMarker_noRibo[!grepl("^mt-", ClusterMarker_noRibo$gene, ignore.case = F),]

有些基因比如Foxp3，对细胞鉴定很重要，但常常在筛选Marker基因的时候筛选不出来。非负矩阵分解可能更好。
参考：过滤线粒体核糖体基因

8. 提取亚群

Idents(pbmc) <- "celltype" #先把celltype设成Idents（active ident)
CD4T <- subset(pbmc, idents = "CD4_T") #通过idents来提取
tmp <- c("orig.ident", "percent.mt","percent.rb","percent.HB","S.Score",
         "G2M.Score","Phase","DF.classifications")
CD4T <- CreateSeuratObject(CD4T@assays$RNA@counts, meta.data = CD4T@meta.data[,tmp])
#metadata保留了，降维聚类结果没有保留

⚠️新提取的亚群需要重新进行降维聚类（和大群相比，标记基因发生了变化），并重新寻找marker基因，重新分群，注释。❗️subset提取子集后，不同样本间批次校正的信息也被去除，需要重新进行批次矫正

参考：Seurat取子集时会用到的函数和方法⚠️⚠️⚠️

9. 取子集后如何去除空子集（还存在这个level，但里面包含的细胞为0，如何去除）as.factor(as.character())

Idents(CD4T) <- "seurat_clusters"
CD4T.sub <- subset(CD4T, idents = c(0,1,2,3,4,5,6))
CD4T.sub$seurat_clusters <- as.factor(as.character(CD4T.sub$seurat_clusters))
CD4T.sub$celltype <- as.factor(as.character(CD4T.sub$celltype))
saveRDS(CD4T.sub, "CD4T/CD4T.classified.rds")

10. 双细胞的预测和去除如DoubletFinder建议单样本进行，不建议双样本一起预测。除此之外，其他步骤都可以多样本一起做，质控的时候也可以多样本一起做，但是建议每个样本都单独看一看。

11. 单细胞多样本整合：merge()；多样本拆分：SplitObject()

scRNA <- merge(scRNAlist[[1]], scRNAlist[2:length(scRNAlist)])
scRNAlist <- SplitObject(scRNA, split.by = "orig.ident")

12. 在做多组数据整合，每个组又有多个样本的时候，最好把单独的每个样本当成批次，而不是把不同的组当成批次。

13. 多核运算

library(future)
options(future.globals.maxSize = 20 * 1024^3)
plan(multisession, workers = 8) #开启多核运算 (8个核)

plan('sequential') #终止多核运算

14. pbmc3k.final@commands$FindClusters 可以查看FindClusters运行时间和记录。Seurat是记录其分析过程的，也可查看command下其他操作

15. 关于质控标准：同一组织的最好用同样的标准，不同组织的可以不一样。不同组织线粒体含量等可能存在差异。

16. 可视化的方法总结
参考：https://www.jianshu.com/p/0d1e2c7d21a4

16. circos图绘制

17. 单细胞数据思维导图，有利于查看单细胞数据格式。
https://www.jianshu.com/p/7560f4fd0d77

18. 对于旧版本Seurat对象的更新
scRNA <- UpdateSeuratObject(scRNA)
UpdateSeuratObject {SeuratObject} ：Update old Seurat object to accommodate new features

19. 对有些操作需要用到python设置的情况

Sys.setenv(RETICULATE_PYTHON="/usr/bin/python3")
# ⚠️要根据自己python3的路径来设置，可以在终端使用which python3来查看

20. 单细胞数据做pooling的好处：可以尽量的降低dropout的问题。（dropout就是矩阵中的zero，这些zero实际上并不是0，而是每个液滴里面起始反应量太低了。而一般的反转录效率只能到30%左右，70%的转录本实际上在反转录那一步是被丢掉的，这是单细胞测序一个比较大的问题）。
但是一旦做了pooling，你必须要证明pooling对结果是没有影响的（下图的右面三个图）。

这个文章用了pooling

21. Seurat的VlnPlot中的combine参数，在如下画三个基因的情况下，设置成T就画一张图，设置成False，会将三个基因各画一张图。

VlnPlot(seuratObj, features = c("Il15", "Cxcl10","Cxcl16"), split.by = "aggregate", pt.size = 0, combine = T)

22. rev()这一步是将横坐标的基因反过来排序

DotPlot(scRNA, features = nichenet_output$top_ligands, cols = "RdYlBu") + RotatedAxis()
DotPlot(seuratObj, features = nichenet_output$top_ligands %>% rev(), cols = "RdYlBu") + RotatedAxis()

这两个画出来的图横坐标基因的顺序是相反的（见NicheNet）

23. 堆叠小提琴图的绘制
完成这个需求有以下几种实现方法：
1. Seurat包直接就可以实现（stack = T）
2. 通过Seuart->scanpy来实现，第一张是Seurat包VlnPlot函数画的图，第二张是scanpy中stacked_violin函数画的图，那么现在问题就变成为Seurat对象到scanpy对象的转换
3. 用R原生函数实现StackedVlnPlot的方法
4. 使用基于scanpy包衍生的scanyuan包来说实现
5. 使用R包MySeuratWrappers来实现

最简单的方法1如下：

pbmc = readRDS('pbmc.rds')
markers <- c('CD14','CD68','CD4','CD19','MS4A1','CD79A','GNLY','NKG7','GZMB','CD8A','CD8B','FCGR3A',
             'FOXP3','PPBP','FCER1A','CD34')
markers <- CaseMatch(markers, rownames(pbmc))
markers <- as.character(markers)
VlnPlot(pbmc, features = markers, pt.size = 0, group.by = 'cell_type',stack = T)+NoLegend()
  1. 当绘图需要设置横坐标顺序时，先把Ident设置为需要绘图的变量，再使用factor进行设置。
Idents(scRNA) <- 'orig.ident'
My_levels <- c( "con1","con2", "con3", "case1","case2", "case3")
Idents(scRNA) <- factor(Idents(scRNA), levels= My_levels)

如果不设置level，会按字母顺序排列，case会自动排在con前面。

  1. 在绘制多个质控参数的小提琴图时，可以先生成一个空list，将每张图都存成list的一个对象，再使用patchwork包的wrap_plots函数画成一张图。
plot.featrures = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.HB")
plots = list()
for(i in seq_along(plot.featrures)){
    plots[[i]] = VlnPlot(scRNA, pt.size = 0,
                         features = plot.featrures[i],log = T) + theme.set2 + NoLegend()}
wrap_plots(plots = plots, nrow=2) 
ggsave(filename = 'beforeQC9log.pdf', width = 8, height = 6)
  1. 在进行细胞注释时，可以使用recode函数
scRNA$celltype <- scRNA$SCT_snn_res.0.3
scRNA$celltype <- recode(scRNA$celltype,
                         "0" = "CD14_Monocyte",
                         "1" = "CD14_Monocyte",
                         "2" = "NK_cell",
                         "3" = "B_cell",
                         "4" = "T_cell",
                         "5" = "T_cell",
                         "6" = "CD16_Monocyte",
                         "7" = "Megakaryocyte",
                         "8" = "NK_cell",
                         "9" = "T_cell",
                         "10" = "T_cell",
                         "11" = "T_cell",
                         "12" = "T_cell",
                         "13" = "T_cell",
                         "14" = "B_cell",
                         "15" = "DC",
                         "16" = "Neutrophil",
                         "17" = "T_cell",
                         "18" = "Proliferating_lymphocyte",
                         "19" = "CD14_Monocyte",
                         "20" = "B_cell")
  1. 对某一metadata进行重命名
#查看meta.data中有哪些内容
colnames(scRNA@meta.data) 
#meta.data中第13个进行重命名
colnames(scRNA@meta.data)[13] <- 'Cytokine_Score'

使用Seurat的RenameIdents函数也可以

scRNA <- RenameIdents(scRNA, 'oldIdent'='newIdent')
  1. 进行marker基因绘图时，使用CaseMatch函数去除不存在的基因。
markers <- c('CD79A','CD79B','MS4A1','CD19','MZB1','GNLY','NKG7','NCAM1','CD3D','CD3E')
markers <- CaseMatch(markers, rownames(scRNA))
markers <- as.character(markers)
  1. prop.table函数，生成细胞比例表格，用于绘制细胞比例条形图
prop.table(x, margin = NULL)

x: table
margin: a vector giving the margins to split by. E.g., for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns. When x has named dimnames, it can be a character vector selecting dimension names.

m <- matrix(1:4, 2)
m
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4
proportions(m, 1)
#           [,1]      [,2]
# [1,] 0.2500000 0.7500000
# [2,] 0.3333333 0.6666667

#根据orig.ident来计算细胞百分比
cell.prop<-as.data.frame(prop.table(table(Idents(scRNA), scRNA$orig.ident)))
colnames(cell.prop)<-c("cluster","Group","proportion")
write.csv(cell.prop,row.names = FALSE,file = 'cell_prop_all.csv')
cell.prop$Group <- factor(cell.prop$Group,ordered=TRUE,levels=c("HC1","HC2","HC3","Case1","Case2","Case3"))
ggplot(cell.prop,aes(Group,proportion,fill=cluster))+
geom_bar(stat="identity",position="fill")+
  ggtitle("")+
  theme_bw()+
  theme(axis.ticks.length=unit(0.5,'cm'))+
  guides(fill=guide_legend(title=NULL))
ggsave(
  'NK/satck_prop_NK.pdf',
  plot = last_plot(),
  device = NULL,
  path = NULL,
  scale = 1,
  width = 10,
  height = 8,
  dpi = 300,
  limitsize = TRUE,
)
  1. 在从基因表达矩阵（行为基因，列为样本）直接生成Seurat对象时，如下：
HC_1 <- CreateSeuratObject(counts = HC_1)

得到的HC_1样本的orig.ident默认是样本名中第一个_号的前一部分。所以要保证矩阵的列名是样本名_细胞barcode这样的格式。
如果有多个分组，例如两个样本矩阵中细胞分别命名为HC_1_barcode，HC_2_barcode，在直接通过如下方法得到两个Seurat对象，再对其进行merge之后，两个样本会被合并成一个。也就是样本信息只保留了第一个_号之前的HC，没有保留_号之后的1和2。

HC_1 <- CreateSeuratObject(counts = HC_1)
HC_2 <- CreateSeuratObject(counts = HC_2)
scRNA <- merge(HC_1,HC_2)
table(scRNA$orig.ident)
# HC
#3000

为了避免这种情况，可以在构建Seurat对象时通过参数进行设置

HC_1 <- CreateSeuratObject(counts = HC_1,names.field = 1:2)
HC_2 <- CreateSeuratObject(counts = HC_2,names.field = 1:2)
scRNA <- merge(HC_1,HC_2)
table(scRNA$orig.ident)
# HC_1  HC_2 
# 1500  1500

  1. 在做umap和tsne的时候有很多聚不在一起的小群，可以考虑下面几种情况
    (1) 双细胞的问题
    双细胞倾向于独立成群，或位于大群的边缘或线状连接部分。
    (2) 高变基因数目选择问题
    FindVariableFeatures()时nfeatures选的越多，图越分散
    (3) 降维时pc数目的选择
    做tsne和umap降维时选择的pc数越多，图越分散。
    参考：如何确定细胞聚类的PC数

  2. 亚群注释注意事项

一般先选默认分辨率(0.8)，大概可能会分出十几个群。因为最终都是要注释到每一个barcode，所以首先可以看大类marker的分布(不受分辨率影响)，可以根据marker基因的分布来调整分辨率。是否需要精细的分群得看精细的分群对研究有没有决定作用，还有很重要的一点是看分出的各个cluster在Findallmarkers给出的结果中marker的热图是不是能明显分开。精细划分的细胞本来就很类似，如果有部分小群的热图明显分不开或者非常类似，就可以考虑把分辨率调小。

比如像这种情况，3和4可能就是同一类细胞
  1. 跑umap/tsne跟跑FindNeighbors输入的都是PCA的矩阵，但是用的pc数可以是不一样的。
    在一些教程里面，我们常常看到这样的代码，在看了一下ElbowPlot之后选定了PC数就直接RunTSNE，RunUMAP和FindNeighbors都用一样的pc数。
pbmc <- RunPCA(pbmc, verbose = F)
ElbowPlot(pbmc, ndims = 50)
pc.num=1:20
pbmc <- pbmc %>% RunTSNE(dims=pc.num) %>% RunUMAP(dims=pc.num) %>%
  FindNeighbors(dims = pc.num) %>% FindClusters()

这实际上是没有必要的必须保持一致的。下游的都是用pca之后的，pca是为了压缩数据。
umap和tsne是为了可视化(仅仅是可视化)，但是FindNeighbor是计算细胞间距离矩阵。找类群数目和可视化可以说没有关系。

  1. UMAP图的分群可视化
library(dplyr)
library(Seurat)
library(patchwork) 
library(tidyverse)
library(cowplot)
library(clustree)

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc,dims=1:20,resolution = seq(from=0,by=.2,length=10))
pbmc <- RunUMAP(pbmc, dims = 1:10)
clustree(pbmc)

Idents(pbmc) <- "RNA_snn_res.1"
p=list()
p<- map(c(levels(Idents(pbmc))),function(x){DimPlot(pbmc, cells.highlight = CellsByIdentities(object = pbmc, idents = x))})
plot_grid(plotlist = p)

map函数：
https://blog.sciencenet.cn/blog-2577109-1290630.html
https://blog.csdn.net/qq_18055167/article/details/1044372363

  1. RunUMAP()不同参数对分群的影响

  2. 批量读入矩阵并创建seurat对象
    ⚠️assign()的用法

# Data loading and QC

### 读入sample消息
samples <- read_excel("../data/metadata/patients_metadata.xlsx", range = cell_cols("A:A")) %>% .$sample_id
samples
#  [1] "p018t" "p018n" "p019t" "p019n" "p023t" "p024t" "p027t" "p027n" "p028n"
# [10] "p029n" "p030t" "p030n" "p031t" "p031n" "p032t" "p032n" "p033t" "p033n"
# [19] "p034t" "p034n"

### import cellranger files from different data sets
for (i in seq_along(samples)){
  assign(paste0("scs_data", i), Read10X(data.dir = paste0("../data/cellranger/", samples[i], "/filtered_feature_bc_matrix")))
}
# 读入的每一个文件都是一个对象，命名为scs_data1-20

### create seurat objects from cellranger files
for (i in seq_along(samples)){
  assign(paste0("seu_obj", i), CreateSeuratObject(counts = eval(parse(text = paste0("scs_data", i))), project = samples[i], min.cells = 3))
}
# 对每一个scs_data创建一个Seurat对象，命名为seu_obj1-20

### merge data sets
seu_obj <- merge(seu_obj1, y = c(seu_obj2, seu_obj3, seu_obj4, seu_obj5, seu_obj6, seu_obj7, seu_obj8, seu_obj9, seu_obj10, seu_obj11, seu_obj12, seu_obj13, seu_obj14, seu_obj15, seu_obj16, seu_obj17, seu_obj18, seu_obj19, seu_obj20), add.cell.ids = samples, project = "lung")

### calculate mitochondrial, hemoglobin and ribosomal gene counts
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^MT-", col.name = "pMT")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^HBA|^HBB", col.name = "pHB")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^RPS|^RPL", col.name = "pRP")

qcparams <- c("nFeature_RNA", "nCount_RNA", "pMT", "pHB", "pRP")
for (i in seq_along(qcparams)){
  print(VlnPlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident", pt.size = 0))
}
for (i in seq_along(qcparams)){
  print(RidgePlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident"))
}

VlnPlot(seu_obj, features = c("nFeature_RNA", "nCount_RNA", "pMT"), pt.size = 0, group.by = "orig.ident", ncol = 1, log = T)
ggsave2("SuppFig1B.pdf", path = "../results", width = 30, height = 20, units = "cm")
  1. 绘图时颜色的传参
ClusterName_color_panel <- c(
  "Naive CD4 T" = "#DC143C", "Memory CD4 T" = "#0000FF", "CD14+ Mono" = "#20B2AA",
  "B" = "#FFA500", "CD8 T" = "#9370DB", "FCGR3A+ Mono" = "#98FB98",
  "NK" = "#F08080", "DC" = "#0000FF", "Platelet" = "#20B2AA"
)
ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +
  geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()+ scale_fill_manual(name = "", values = ClusterName_color_panel)+scale_color_manual(name = "", values = ClusterName_color_panel)

参考：monocle2

  1. 绘图气泡图/小提琴图的时候翻转纵坐标
    使用scale_y_discrete+rev 
p1=DotPlot(pbmc,features = c('CCR2','CD3D'),group.by = 'cell_type')+theme_bw()
p2=DotPlot(pbmc,features = c('CCR2','CD3D'),group.by = 'cell_type')+theme_bw()+scale_y_discrete(limits = rev(levels(pbmc$cell_type)))
p1|p2
上一篇 下一篇

猜你喜欢

热点阅读