本文是参考学习 CNS图表复现02—Seurat标准流程之聚类分群的学习笔记。可能根据学习情况有所改动。

今天讲解第二步：完成Seurat标准流程之聚类分群。

直接上代码：

> load(file = "main_tiss_filtered.RData") #加载   load之后右侧environment就可以看到变量名20210109
Loading required package: Seurat
Error: package or namespace load failed for ‘Seurat’ in .doLoadActions(where, attach):
 error in load action .__A__.1 for package RcppAnnoy: loadModule(module = "AnnoyAngular", what = TRUE, env = ns, loadNow = TRUE): Unable to load module "AnnoyAngular": attempt to apply non-function
Error in .requirePackage(package) : 
  unable to find required package ‘Seurat’
In addition: Warning message:
package ‘Seurat’ was built under R version 4.0.3 
Error: no more error handlers available (recursive errors?); invoking 'abort' restart

错了，重来,加上library(Seurat)

今天讲解第二步：完成Seurat标准流程之聚类分群。

直接上代码：

> library(Seurat)

Seurat v4 will be going to CRAN in the near future;
 for more details, please visit https://satijalab.org/seurat/v4_changes

Warning message:
程辑包‘Seurat’是用R版本4.0.3 来建造的 
> load(file = "main_tiss_filtered.RData") #加载   load之后右侧environment就可以看到变量名20210109
> raw_sce <- main_tiss_filtered
> raw_sce
An object of class Seurat 
26577 features across 21620 samples within 1 assay 
Active assay: RNA (26577 features, 0 variable features)
> rownames(raw_sce)[grepl('^mt-',rownames(raw_sce),ignore.case = T)]
character(0)
> rownames(raw_sce)[grepl('^Rp[sl]',rownames(raw_sce),ignore.case = T)]
  [1] "RPL10"          "RPL10A"         "RPL10L"         "RPL11"          "RPL12"         
  [6] "RPL13"          "RPL13A"         "RPL13AP17"      "RPL13AP20"      "RPL13AP3"      
 [11] "RPL13AP5"       "RPL13AP6"       "RPL13P5"        "RPL14"          "RPL15"         
 [16] "RPL17"          "RPL17-C18orf32" "RPL18"          "RPL18A"         "RPL19"         
 [21] "RPL19P12"       "RPL21"          "RPL21P28"       "RPL21P44"       "RPL22"         
 [26] "RPL22L1"        "RPL23"          "RPL23A"         "RPL23AP32"      "RPL23AP53"     
 [31] "RPL23AP64"      "RPL23AP7"       "RPL23AP82"      "RPL23AP87"      "RPL23P8"       
 [36] "RPL24"          "RPL26"          "RPL26L1"        "RPL27"          "RPL27A"        
 [41] "RPL28"          "RPL29"          "RPL29P2"        "RPL3"           "RPL30"         
 [46] "RPL31"          "RPL31P11"       "RPL32"          "RPL32P3"        "RPL34"         
 [51] "RPL34-AS1"      "RPL35"          "RPL35A"         "RPL36"          "RPL36A"        
 [56] "RPL36A-HNRNPH2" "RPL36AL"        "RPL37"          "RPL37A"         "RPL38"         
 [61] "RPL39"          "RPL39L"         "RPL3L"          "RPL4"           "RPL41"         
 [66] "RPL5"           "RPL6"           "RPL7"           "RPL7A"          "RPL7L1"        
 [71] "RPL8"           "RPL9"           "RPLP0"          "RPLP0P2"        "RPLP1"         
 [76] "RPLP2"          "RPS10"          "RPS10-NUDT3"    "RPS10P7"        "RPS11"         
 [81] "RPS12"          "RPS13"          "RPS14"          "RPS14P3"        "RPS15"         
 [86] "RPS15A"         "RPS15AP10"      "RPS16"          "RPS16P5"        "RPS17"         
 [91] "RPS18"          "RPS18P9"        "RPS19"          "RPS19BP1"       "RPS2"          
 [96] "RPS20"          "RPS21"          "RPS23"          "RPS24"          "RPS25"         
[101] "RPS26"          "RPS26P11"       "RPS27"          "RPS27A"         "RPS27L"        
[106] "RPS28"          "RPS29"          "RPS2P32"        "RPS3"           "RPS3A"         
[111] "RPS4X"          "RPS4Y1"         "RPS4Y2"         "RPS5"           "RPS6"          
[116] "RPS6KA1"        "RPS6KA2"        "RPS6KA2-AS1"    "RPS6KA2-IT1"    "RPS6KA3"       
[121] "RPS6KA4"        "RPS6KA5"        "RPS6KA6"        "RPS6KB1"        "RPS6KB2"       
[126] "RPS6KC1"        "RPS6KL1"        "RPS7"           "RPS7P5"         "RPS8"          
[131] "RPS9"           "RPSA"           "RPSAP52"        "RPSAP58"        "RPSAP9"        
> rownames(raw_sce)[grepl('^MT-',rownames(raw_sce),ignore.case = T)]
character(0)
> rownames(raw_sce)[grepl('^Rp[sl]',rownames(raw_sce),ignore.case = T)]
  [1] "RPL10"          "RPL10A"         "RPL10L"         "RPL11"          "RPL12"         
  [6] "RPL13"          "RPL13A"         "RPL13AP17"      "RPL13AP20"      "RPL13AP3"      
 [11] "RPL13AP5"       "RPL13AP6"       "RPL13P5"        "RPL14"          "RPL15"         
 [16] "RPL17"          "RPL17-C18orf32" "RPL18"          "RPL18A"         "RPL19"         
 [21] "RPL19P12"       "RPL21"          "RPL21P28"       "RPL21P44"       "RPL22"         
 [26] "RPL22L1"        "RPL23"          "RPL23A"         "RPL23AP32"      "RPL23AP53"     
 [31] "RPL23AP64"      "RPL23AP7"       "RPL23AP82"      "RPL23AP87"      "RPL23P8"       
 [36] "RPL24"          "RPL26"          "RPL26L1"        "RPL27"          "RPL27A"        
 [41] "RPL28"          "RPL29"          "RPL29P2"        "RPL3"           "RPL30"         
 [46] "RPL31"          "RPL31P11"       "RPL32"          "RPL32P3"        "RPL34"         
 [51] "RPL34-AS1"      "RPL35"          "RPL35A"         "RPL36"          "RPL36A"        
 [56] "RPL36A-HNRNPH2" "RPL36AL"        "RPL37"          "RPL37A"         "RPL38"         
 [61] "RPL39"          "RPL39L"         "RPL3L"          "RPL4"           "RPL41"         
 [66] "RPL5"           "RPL6"           "RPL7"           "RPL7A"          "RPL7L1"        
 [71] "RPL8"           "RPL9"           "RPLP0"          "RPLP0P2"        "RPLP1"         
 [76] "RPLP2"          "RPS10"          "RPS10-NUDT3"    "RPS10P7"        "RPS11"         
 [81] "RPS12"          "RPS13"          "RPS14"          "RPS14P3"        "RPS15"         
 [86] "RPS15A"         "RPS15AP10"      "RPS16"          "RPS16P5"        "RPS17"         
 [91] "RPS18"          "RPS18P9"        "RPS19"          "RPS19BP1"       "RPS2"          
 [96] "RPS20"          "RPS21"          "RPS23"          "RPS24"          "RPS25"         
[101] "RPS26"          "RPS26P11"       "RPS27"          "RPS27A"         "RPS27L"        
[106] "RPS28"          "RPS29"          "RPS2P32"        "RPS3"           "RPS3A"         
[111] "RPS4X"          "RPS4Y1"         "RPS4Y2"         "RPS5"           "RPS6"          
[116] "RPS6KA1"        "RPS6KA2"        "RPS6KA2-AS1"    "RPS6KA2-IT1"    "RPS6KA3"       
[121] "RPS6KA4"        "RPS6KA5"        "RPS6KA6"        "RPS6KB1"        "RPS6KB2"       
[126] "RPS6KC1"        "RPS6KL1"        "RPS7"           "RPS7P5"         "RPS8"          
[131] "RPS9"           "RPSA"           "RPSAP52"        "RPSAP58"        "RPSAP9"        
> raw_sce[["percent.mt"]] <- PercentageFeatureSet(raw_sce, pattern = "^MT-")
> fivenum(raw_sce[["percent.mt"]][,1])
[1] 0 0 0 0 0
> rb.genes <- rownames(raw_sce)[grep("^RP[SL]",rownames(raw_sce),ignore.case = T)]
> C<-GetAssayData(object = raw_sce, slot = "counts")
> percent.ribo <- Matrix::colSums(C[rb.genes,])/Matrix::colSums(C)*100
> fivenum(percent.ribo)
A12_B001464 L19_B003105  M3_B001543 I21_B003528  E5_B003659 
   0.000000    2.196870    3.409555    5.444660   49.341911 
> raw_sce <- AddMetaData(raw_sce, percent.ribo, col.name = "percent.ribo")
> plot1 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "percent.ercc")
Error: Feature 2 (percent.ercc) not found.
In addition: Warning message:
In FetchData(object = object, vars = c(feature1, feature2, group.by),  :
  The following requested variables were not found: percent.ercc
> plot2 <- FeatureScatter(raw_sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
> CombinePlots(plots = list(plot1, plot2))
Error in CombinePlots(plots = list(plot1, plot2)) : 
  object 'plot1' not found
In addition: Warning message:
CombinePlots is being deprecated. Plots should now be combined using the patchwork system. 
> VlnPlot(raw_sce, features = c("percent.ribo", "percent.ercc"), ncol = 2)
Warning message:
In FetchData(object = object, vars = features, slot = slot) :
  The following requested variables were not found: percent.ercc
> VlnPlot(raw_sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
> VlnPlot(raw_sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)
> raw_sce
An object of class Seurat 
26577 features across 21620 samples within 1 assay 
Active assay: RNA (26577 features, 0 variable features)
> sce=raw_sce
> sce
An object of class Seurat 
26577 features across 21620 samples within 1 assay 
Active assay: RNA (26577 features, 0 variable features)
> sce <- NormalizeData(sce, normalization.method =  "LogNormalize", 
+                      scale.factor = 10000)
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> GetAssay(sce,assay = "RNA")
Assay data with 26577 features for 21620 cells
First 10 features:
 A1BG, A1BG-AS1, A1CF, A2M, A2M-AS1, A2ML1, A2MP1, A3GALT2, A4GALT, A4GNT 
> sce <- FindVariableFeatures(sce, 
+                             selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> # 步骤 ScaleData 的耗时取决于电脑系统配置（保守估计大于一分钟）
> sce <- ScaleData(sce)
Centering and scaling data matrix
  |=======================================================================================| 100%
> sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
PC_ 1 
Positive:  CORO1A, CXCR4, IL2RG, CD52, RHOH, ALOX5AP, GPR183, CD69, ISG20, LTB 
       CST7, CCL5, LCK, UCP2, FERMT3, SERPINB9, ARHGAP30, TUBA4A, CD247, AMICA1 
       CD27, CCR7, NKG7, TBC1D10C, SASH3, S1PR4, SELL, INPP5D, CTSW, TRAF3IP3 
Negative:  SPARC, CALD1, DCN, COL1A2, IGFBP7, LUM, MGP, COL3A1, THY1, RARRES2 
       FBLN1, MFAP4, SPARCL1, COL1A1, TIMP3, TPM2, CNN3, CYR61, TAGLN, SERPINF1 
       LHFP, CTSK, PDGFRB, CTGF, CD248, CRISPLD2, PRELP, COL5A1, ACTA2, OLFML3 
PC_ 2 
Positive:  TSPAN1, SLPI, KRT18, PIGR, RSPH1, SLC34A2, PSCA, PIFO, SNTN, AGR2 
       C20orf85, FAM183A, CAPS, C9orf24, TMEM190, LDLRAD1, CAPSL, C1orf194, ZMYND10, CCDC78 
       C11orf88, TEKT1, WDR38, ROPN1L, RSPH9, FAM92B, TEKT2, DEGS2, TUBA4B, LCN2 
Negative:  CORO1A, CXCR4, DCN, COL1A2, A2M, LUM, SPARC, ZEB2, COL3A1, FN1 
       THY1, IL2RG, CD52, SERPINB9, MFAP4, PDGFRB, ALOX5AP, CTSK, TAGLN, CALD1 
       SERPINF1, ERCC-00171, GPR183, OLFML3, SPARCL1, CD248, PRELP, RHOH, SFRP2, CCND2 
PC_ 3 
Positive:  NAPSA, SFTPB, RNASE1, SERPINA1, SFTPA1, SFTPA2, SFTPD, SFTA3, C4BPA, CEACAM6 
       C16orf89, SLC22A31, PON3, SCGB3A2, EFNA1, SLC34A2, LGMN, AQP4, ABCA3, PEBP4 
       SCGB3A1, S100A9, SFTPC, SCD, PGC, HSD17B6, CTSE, FTL, SUSD2, PON2 
Negative:  C20orf85, C9orf24, TMEM190, SNTN, CAPSL, C1orf194, FAM183A, RSPH1, C11orf88, ZMYND10 
       LDLRAD1, TEKT1, WDR38, TUBA4B, FAM92B, ROPN1L, RSPH9, CCDC78, PIFO, C2orf40 
       TEKT2, CAPS, C22orf15, SPAG8, PSCA, TPPP3, WDR54, GSTA1, CRIP1, SCGB2A1 
PC_ 4 
Positive:  FCGR2A, MS4A7, FTL, IL1B, TREM2, MS4A4A, OLR1, CLEC7A, APOE, APOC1 
       MARCKS, GPNMB, TGFBI, FOLR2, CPVL, IL1RN, CXCL3, SPP1, HMOX1, HLA-DMB 
       ZEB2, FCN1, CFD, NLRP3, CCL2, S100A8, PLA2G7, SLC8A1, C20orf85, TMEM190 
Negative:  IL32, TUBA4A, LCK, CD247, PRF1, IL2RG, CCL5, NKG7, OCIAD2, CD96 
       CTSW, RHOH, GZMM, HOPX, ZAP70, GZMA, SH2D1A, CD8A, CST7, CCND3 
       CXCR6, FAM46C, CXCR3, TBCC, CD27, TBC1D10C, EFNA1, LEPROTL1, CD69, GZMH 
PC_ 5 
Positive:  FBLN1, DCN, SFRP2, SERPINF1, LUM, CTSK, COL1A2, COL3A1, COL1A1, RARRES2 
       SFRP4, OLFML3, MFAP4, CXCL14, EFEMP1, IGF1, ADH1B, MFAP5, COL5A1, DPT 
       HTRA3, FGF7, SULF1, CRISPLD2, FNDC1, SLIT2, COL12A1, FBLN5, WISP2, PRELP 
Negative:  CLEC14A, RAMP3, CLDN5, RAMP2, VWF, CDH5, ESAM, ECSCR, PLVAP, SOX18 
       EGFL7, HYAL2, GNG11, FAM107A, AQP1, CXorf36, CD34, FCN3, SDPR, ACKR1 
       CLEC3B, PODXL, KANK3, TEK, COL4A1, NES, AFAP1L1, TINAGL1, ARHGEF15, C2CD4B 
> DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)
> ElbowPlot(sce)
> sce <- FindNeighbors(sce, dims = 1:15)
Computing nearest neighbor graph
Computing SNN
> sce <- FindClusters(sce, resolution = 0.2)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 21620
Number of edges: 759616

Running Louvain algorithm...
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9684
Number of communities: 17
Elapsed time: 4 seconds
> table(sce@meta.data$RNA_snn_res.0.2)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
4447 4001 2661 2381 1881 1406  997  957  740  554  434  407  219  184  173  130   48 
> sce <- FindClusters(sce, resolution = 0.8)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 21620
Number of edges: 759616

Running Louvain algorithm...
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9276
Number of communities: 30
Elapsed time: 5 seconds
> table(sce@meta.data$RNA_snn_res.0.8)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18 
2148 2075 1652 1607 1208 1045 1020  999  997  971  957  891  739  725  623  554  511  434  405 
  19   20   21   22   23   24   25   26   27   28   29 
 380  333  231  219  189  173  136  130  118  102   48 
> library(gplots)

载入程辑包：‘gplots’

The following object is masked from ‘package:stats’:

    lowess

Warning message:
程辑包‘gplots’是用R版本4.0.3 来建造的 
> tab.1=table(sce@meta.data$RNA_snn_res.0.2,sce@meta.data$RNA_snn_res.0.8)
> balloonplot(tab.1)
> set.seed(123)
> sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
> DimPlot(sce,reduction = "tsne",label=T)
Warning message:
Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
Please use `as_label()` or `as_name()` instead.
This warning is displayed once per session. 
> phe=data.frame(cell=rownames(sce@meta.data),
+                cluster =sce@meta.data$seurat_clusters)
> head(phe)
            cell cluster
1 A10_1001000329       2
2 A10_1001000407      10
3 A10_1001000408      10
4 A10_1001000410       1
5 A10_1001000412      10
6    A10_B000420      16
> table(phe$cluster)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18 
2148 2075 1652 1607 1208 1045 1020  999  997  971  957  891  739  725  623  554  511  434  405 
  19   20   21   22   23   24   25   26   27   28   29 
 380  333  231  219  189  173  136  130  118  102   48 
> tsne_pos=Embeddings(sce,'tsne')
> DimPlot(sce,reduction = "tsne",group.by  ='orig.ident')
> DimPlot(sce,reduction = "tsne",label=T,split.by ='orig.ident')
> head(phe)
            cell cluster
1 A10_1001000329       2
2 A10_1001000407      10
3 A10_1001000408      10
4 A10_1001000410       1
5 A10_1001000412      10
6    A10_B000420      16
> table(phe$cluster)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18 
2148 2075 1652 1607 1208 1045 1020  999  997  971  957  891  739  725  623  554  511  434  405 
  19   20   21   22   23   24   25   26   27   28   29 
 380  333  231  219  189  173  136  130  118  102   48 
> head(tsne_pos)
                   tSNE_1     tSNE_2
A10_1001000329 -19.916177 -20.300083
A10_1001000407 -26.883484  -8.166466
A10_1001000408 -37.016645 -14.658582
A10_1001000410  23.665744 -15.812134
A10_1001000412 -37.017060 -11.048993
A10_B000420     -6.313442 -18.903210
> dat=cbind(tsne_pos,phe)
> pro='first'
> save(dat,file=paste0(pro,'_for_tSNE.pos.Rdata'))
> load(file=paste0(pro,'_for_tSNE.pos.Rdata'))
> library(ggplot2)
> p=ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
> p=p+stat_ellipse(data=dat,aes(x=tSNE_1,y=tSNE_2,fill=cluster,color=cluster),
+                  geom = "polygon",alpha=0.2,level=0.9,type="t",linetype = 2,show.legend = F)+coord_fixed()
> print(p)
Warning message:
In MASS::cov.trob(data[, vars]) : Probable convergence failure
> theme= theme(panel.grid =element_blank()) +   ## 删去网格
+   theme(panel.border = element_blank(),panel.background = element_blank()) +   ## 删去外层边框
+   theme(axis.line = element_line(size=1, colour = "black"))
> p=p+theme+guides(colour = guide_legend(override.aes = list(size=5)))
> print(p)
Warning message:
In MASS::cov.trob(data[, vars]) : Probable convergence failure
> ggplot2::ggsave(filename = paste0(pro,'_tsne_res0.8.pdf'))
Saving 6.4 x 3.77 in image
Warning message:
In MASS::cov.trob(data[, vars]) : Probable convergence failure
> sce <- RunUMAP(object = sce, dims = 1:15, do.fast = TRUE)
Warning: The following arguments are not used: do.fast
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
20:53:01 UMAP embedding parameters a = 0.9922 b = 1.112
20:53:01 Read 21620 rows and found 15 numeric columns
20:53:01 Using Annoy for neighbor search, n_neighbors = 30
20:53:01 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:53:07 Writing NN index file to temp file C:\Users\Nano\AppData\Local\Temp\RtmpGkHwMb\file2032223c6e
20:53:07 Searching Annoy index using 1 thread, search_k = 3000
20:53:18 Annoy recall = 100%
20:53:18 Commencing smooth kNN distance calibration using 1 thread
20:53:21 Initializing from normalized Laplacian + noise
20:53:27 Commencing optimization for 200 epochs, with 909916 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:54:05 Optimization finished
> DimPlot(sce,reduction = "umap",label=T)
> DimPlot(sce,reduction = "umap",group.by = 'orig.ident')
> plot1 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "percent.mt")
Warning message:
In cor(x = data[, 1], y = data[, 2]) : 标准差为零
> plot2 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
> CombinePlots(plots = list(plot1, plot2))
Warning message:
CombinePlots is being deprecated. Plots should now be combined using the patchwork system. 
> ggplot2::ggsave(filename = paste0(pro,'_CombinePlots.pdf'))
Saving 6.4 x 3.77 in image
> VlnPlot(sce, features = c("percent.ribo", "percent.mt"), ncol = 2)
Warning message:
In SingleExIPlot(type = type, data = data[, x, drop = FALSE], idents = idents,  :
  All cells have the same value of percent.mt.
> ggplot2::ggsave(filename = paste0(pro,'_mt-and-ribo.pdf'))
Saving 6.4 x 3.77 in image
> VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
> ggplot2::ggsave(filename = paste0(pro,'_counts-and-feature.pdf'))
Saving 6.4 x 3.77 in image
> VlnPlot(sce, features = c("percent.ribo", "nCount_RNA"), ncol = 2)
> table(sce@meta.data$seurat_clusters)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18 
2148 2075 1652 1607 1208 1045 1020  999  997  971  957  891  739  725  623  554  511  434  405 
  19   20   21   22   23   24   25   26   27   28   29 
 380  333  231  219  189  173  136  130  118  102   48 

下面这一步时间较长

16G内存电脑跑了2个小时

> for( i in unique(sce@meta.data$seurat_clusters) ){
+   markers_df <- FindMarkers(object = sce, ident.1 = i, 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 )
+   ggsave(filename=paste0(pro,'_VlnPlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
+   FeaturePlot(object = sce, features=markers_genes )
+   ggsave(filename=paste0(pro,'_FeaturePlot_subcluster_',i,'_sce.markers_heatmap.pdf'))
+ }
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 24s
      p_val avg_logFC pct.1 pct.2 p_val_adj
LYZ       0  2.371829 0.973 0.240         0
FCN1      0  2.295688 0.651 0.063         0
IL1B      0  2.164884 0.809 0.163         0
EREG      0  1.880129 0.539 0.081         0
OLR1      0  1.755714 0.697 0.119         0
CXCL3     0  1.593757 0.604 0.210         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07m 48s
       p_val avg_logFC pct.1 pct.2 p_val_adj
LCN2       0  3.340067     1 0.093         0
MUC20      0  2.946946     1 0.133         0
CD24       0  2.855278     1 0.153         0
KRT7       0  2.848640     1 0.168         0
SCCPDH     0  2.700061     1 0.220         0
WFDC2      0  2.677898     1 0.202         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 44s
       p_val avg_logFC pct.1 pct.2 p_val_adj
IL7R       0  2.109133 0.813 0.144         0
CCR7       0  1.834167 0.598 0.105         0
LCK        0  1.677304 0.708 0.093         0
CXCR4      0  1.622574 0.893 0.395         0
SARAF      0  1.604799 0.971 0.756         0
SPOCK2     0  1.527563 0.731 0.107         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 54s
         p_val avg_logFC pct.1 pct.2 p_val_adj
CD1C         0  2.012028 0.499 0.022         0
NAPSB        0  1.974966 0.822 0.134         0
HLA-DQA1     0  1.857886 0.975 0.353         0
CD1E         0  1.829083 0.472 0.014         0
FCER1A       0  1.728109 0.419 0.028         0
CLEC10A      0  1.233114 0.499 0.058         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 13s
     p_val avg_logFC pct.1 pct.2 p_val_adj
SPP1     0  3.180647 0.731 0.120         0
C1QB     0  2.705537 0.924 0.088         0
APOE     0  2.656615 0.878 0.229         0
C1QA     0  2.558506 0.942 0.090         0
CD14     0  2.246694 0.975 0.175         0
C1QC     0  2.182844 0.932 0.074         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 56s
        p_val avg_logFC pct.1 pct.2 p_val_adj
NAPSA       0  2.023432 0.621 0.107         0
AZGP1       0  1.741359 0.336 0.027         0
EPCAM       0  1.716374 0.769 0.196         0
KRT19       0  1.706181 0.842 0.238         0
CEACAM6     0  1.678523 0.709 0.158         0
KRT18       0  1.625329 0.812 0.260         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 53s
     p_val avg_logFC pct.1 pct.2 p_val_adj
GNLY     0  3.114463 0.322 0.019         0
CCL5     0  2.708346 0.949 0.133         0
NKG7     0  2.681242 0.801 0.058         0
PRF1     0  2.647103 0.706 0.052         0
CTSW     0  2.243186 0.660 0.058         0
GZMB     0  2.220752 0.441 0.036         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 09s
        p_val avg_logFC pct.1 pct.2 p_val_adj
SCGB3A2     0  2.781470 0.608 0.042         0
SCGB3A1     0  2.650964 0.829 0.067         0
SFTPB       0  2.519461 0.978 0.104         0
AQP4        0  2.478347 0.825 0.043         0
SFTPD       0  2.159180 0.832 0.052         0
C4BPA       0  2.063080 0.819 0.052         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 10s
       p_val avg_logFC pct.1 pct.2 p_val_adj
RGS5       0  3.199030 0.850 0.046         0
ACTA2      0  3.025543 0.845 0.150         0
COL4A1     0  2.565525 0.863 0.108         0
THY1       0  2.555382 0.692 0.091         0
IGFBP7     0  2.550436 1.000 0.315         0
TAGLN      0  2.445983 0.824 0.159         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 57s
      p_val avg_logFC pct.1 pct.2 p_val_adj
CLDN5     0  3.735707 0.854 0.010         0
ACKR1     0  3.552193 0.356 0.007         0
RAMP3     0  2.997716 0.799 0.009         0
HYAL2     0  2.955281 0.877 0.196         0
VWF       0  2.936765 0.788 0.018         0
AQP1      0  2.864674 0.773 0.096         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 49s
        p_val avg_logFC pct.1 pct.2 p_val_adj
IGLL5       0  5.361883 0.962 0.033         0
JCHAIN      0  5.272503 0.866 0.056         0
MZB1        0  4.222333 0.994 0.047         0
DERL3       0  3.083926 0.981 0.089         0
HERPUD1     0  2.982221 0.984 0.645         0
SSR4        0  2.940790 0.995 0.776         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 44s
       p_val avg_logFC pct.1 pct.2 p_val_adj
TPSAB1     0  6.019722 1.000 0.010         0
TPSB2      0  5.329733 0.998 0.008         0
CPA3       0  3.712487 0.998 0.005         0
MS4A2      0  3.448981 0.991 0.006         0
CTSG       0  3.240622 0.502 0.003         0
TPSD1      0  2.874191 0.713 0.003         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 07s
       p_val avg_logFC pct.1 pct.2 p_val_adj
MMP11      0  3.735834 0.684 0.034         0
COL3A1     0  3.715876 0.987 0.083         0
COL1A2     0  3.642927 0.996 0.107         0
SPARC      0  3.066119 0.994 0.209         0
LUM        0  2.895565 0.870 0.077         0
SFRP2      0  2.722759 0.724 0.037         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 28s
      p_val avg_logFC pct.1 pct.2 p_val_adj
MFAP4     0  3.935752 0.954 0.054         0
INMT      0  3.328297 0.742 0.047         0
MGP       0  3.267342 0.970 0.135         0
CFD       0  3.164279 0.738 0.255         0
DCN       0  3.098411 0.959 0.083         0
FBLN1     0  3.043089 0.834 0.137         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=08m 37s
       p_val avg_logFC pct.1 pct.2 p_val_adj
CSF3R      0  3.018346 0.780 0.157         0
G0S2       0  2.801950 0.729 0.244         0
ADGRG3     0  2.746309 0.500 0.014         0
S100A8     0  2.636231 0.785 0.176         0
PROK2      0  2.427165 0.411 0.018         0
FCGR3B     0  2.365168 0.611 0.037         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~14s          Error in UseMethod("depth") : 
  no applicable method for 'depth' applied to an object of class "NULL"
In addition: Warning messages:
1: In grid.Call.graphics(C_setviewport, vp, TRUE) :
  'layout.pos.row'的值不对
2: In doTryCatch(return(expr), name, parentenv, handler) :
  无法弹到最上层的視窗('grid'和'graphics'输出有混合？)
3: In UseMethod("depth") :
  no applicable method for 'depth' applied to an object of class "NULL"
Error: no more error handlers available (recursive errors?); invoking 'abort' restart
Graphics error: Plot rendering error
Error in UseMethod("depth") : 
  no applicable method for 'depth' applied to an object of class "NULL"
Error: no more error handlers available (recursive errors?); invoking 'abort' restart
Graphics error: Plot rendering error
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06m 41s
        p_val avg_logFC pct.1 pct.2 p_val_adj
KRT13       0  4.230822 0.880 0.016         0
KRT6A       0  3.914537 0.943 0.022         0
ALDH3A1     0  3.790334 0.886 0.049         0
AKR1B10     0  3.135880 0.825 0.021         0
AKR1C2      0  3.017500 0.945 0.062         0
AKR1C3      0  2.991269 0.893 0.107         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 59s
       p_val avg_logFC pct.1 pct.2 p_val_adj
MS4A1      0  2.154018 0.583 0.021         0
TCL1A      0  2.086960 0.274 0.007         0
NAPSB      0  1.864919 0.753 0.133         0
SPIB       0  1.775952 0.636 0.046         0
LILRA4     0  1.755557 0.311 0.027         0
BCL11A     0  1.603087 0.772 0.106         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 08s
                p_val  avg_logFC pct.1 pct.2     p_val_adj
PSAP    2.338026e-220 -1.3965945 0.331 0.850 6.213770e-216
CTSD    5.511924e-209 -1.5798582 0.200 0.743 1.464904e-204
IFITM3  3.069978e-187 -1.8009139 0.350 0.770 8.159080e-183
CD63    1.130407e-183 -0.7787599 0.352 0.822 3.004283e-179
ITM2B   5.809383e-179 -0.7746119 0.474 0.910 1.543960e-174
LAPTM4A 5.368902e-176 -1.0974041 0.217 0.709 1.426893e-171
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04m 37s
      p_val avg_logFC pct.1 pct.2 p_val_adj
ALB       0  5.481862 0.985 0.013         0
FGB       0  4.388274 0.646 0.010         0
AMBP      0  4.179769 0.954 0.010         0
APOA2     0  4.032604 0.631 0.009         0
APOA1     0  4.021442 0.600 0.017         0
VTN       0  3.921150 0.946 0.017         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02m 13s
      p_val avg_logFC pct.1 pct.2 p_val_adj
MKI67     0 1.6254500 0.905 0.053         0
TOP2A     0 1.0759263 0.762 0.059         0
BIRC5     0 0.9886379 0.788 0.060         0
RRM2      0 0.9704843 0.619 0.039         0
TPX2      0 0.9488848 0.709 0.050         0
AURKB     0 0.9137951 0.582 0.031         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06m 06s
       p_val avg_logFC pct.1 pct.2 p_val_adj
SFTPC      0  6.187119 0.988 0.039         0
SFTPA2     0  4.745303 0.988 0.063         0
SFTPA1     0  4.371754 0.991 0.059         0
SFTPB      0  3.764590 1.000 0.115         0
SFTPD      0  3.395693 0.991 0.060         0
PGC        0  3.343723 0.933 0.026         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03m 46s
         p_val avg_logFC pct.1 pct.2 p_val_adj
SCGB3A2      0  3.150003 0.970 0.052         0
SFTPA1       0  2.937482 0.996 0.068         0
NAPSA        0  2.693990 0.991 0.122         0
C16orf89     0  2.628919 0.970 0.112         0
C4BPA        0  2.539893 0.974 0.068         0
HPGD         0  2.522995 0.987 0.143         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07m 41s
         p_val avg_logFC pct.1 pct.2 p_val_adj
TPPP3        0  4.134702 1.000 0.110         0
TSPAN1       0  3.701205 1.000 0.157         0
C20orf85     0  3.642036 1.000 0.006         0
CAPS         0  3.583035 1.000 0.102         0
TMEM190      0  3.490252 0.993 0.007         0
RSPH1        0  3.293992 0.993 0.041         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05m 19s
       p_val avg_logFC pct.1 pct.2 p_val_adj
CXCL13     0  2.965052 0.347 0.013         0
MKI67      0  2.328364 0.973 0.051         0
RRM2       0  2.035326 0.817 0.036         0
AURKB      0  1.594833 0.744 0.028         0
CDC20      0  1.581626 0.626 0.042         0
ASF1B      0  1.557680 0.831 0.058         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05m 17s
        p_val avg_logFC pct.1 pct.2 p_val_adj
MUC5B       0  3.877584 0.985 0.064         0
DMBT1       0  2.888433 0.688 0.037         0
CEACAM6     0  2.716520 0.973 0.172         0
AGR2        0  2.702805 0.979 0.170         0
LRIG3       0  2.116126 0.955 0.085         0
TNC         0  1.950363 0.889 0.082         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05m 39s
       p_val avg_logFC pct.1 pct.2 p_val_adj
AGER       0  5.424912 1.000 0.058         0
CYP4B1     0  3.777285 0.960 0.081         0
CLDN18     0  3.359197 0.983 0.047         0
UPK3B      0  2.764430 0.896 0.045         0
SUSD2      0  2.528722 0.908 0.092         0
RTKN2      0  2.460933 0.896 0.063         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06m 04s
                p_val avg_logFC pct.1 pct.2     p_val_adj
TCL1A    0.000000e+00 3.1446182 0.792 0.013  0.000000e+00
AICDA    0.000000e+00 1.8957092 0.729 0.006  0.000000e+00
PAX5     0.000000e+00 0.7329505 0.917 0.024  0.000000e+00
SNX29P2 1.053131e-300 0.8000521 0.750 0.018 2.798906e-296
AURKB   7.249962e-298 2.0674165 1.000 0.034 1.926822e-293
DTX1    5.860053e-277 0.5553459 0.750 0.019 1.557426e-272
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=12m 12s
               p_val avg_logFC pct.1 pct.2     p_val_adj
COL1A1 2.893780e-110  5.474136 0.941 0.204 7.690800e-106
COL3A1 1.432534e-108  1.854252 0.882 0.120 3.807246e-104
TWIST1 1.613721e-108  2.407123 0.490 0.044 4.288785e-104
COL1A2  6.068805e-82  1.138744 0.873 0.143  1.612906e-77
B2M     2.067944e-53 -2.357415 0.843 0.989  5.495975e-49
CFL1    7.346010e-52 -1.995883 0.588 0.951  1.952349e-47
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04m 13s
                p_val avg_logFC pct.1 pct.2     p_val_adj
DLK1     0.000000e+00 2.9797640 0.309 0.003  0.000000e+00
ASCL1    0.000000e+00 0.8161364 0.309 0.002  0.000000e+00
INSM1    0.000000e+00 0.6706563 0.265 0.003  0.000000e+00
SIX3     0.000000e+00 0.4142789 0.257 0.004  0.000000e+00
ZIC2    4.602548e-295 0.8194995 0.279 0.006 1.223219e-290
ADCYAP1 4.918342e-257 0.6012040 0.272 0.007 1.307148e-252
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05m 09s
      p_val avg_logFC pct.1 pct.2 p_val_adj
PMEL      0  5.042566 0.932 0.069         0
MLANA     0  4.037025 0.932 0.050         0
TYRP1     0  4.023393 0.932 0.042         0
DCT       0  3.715292 0.932 0.033         0
TYR       0  3.243960 0.932 0.029         0
BCAN      0  2.605913 0.932 0.033         0
Saving 6.4 x 3.77 in image
Saving 6.4 x 3.77 in image

找marker也耗时近30min

sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25, 
                              thresh.use = 0.25)

DT::datatable(sce.markers)
write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
library(dplyr) 
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)
DoHeatmap(sce,top10$gene,size=3)
ggsave(filename=paste0(pro,'_sce.markers_heatmap.pdf'))

save(sce,sce.markers,file = 'first_sce.Rdata')