single-cell_R_seurat(KOBIC)

00. Set-up

rm(list=ls())
library(Seurat)
library(SeuratData)
library(png)
library(ggplot2)
library(patchwork)
library(dplyr)
library(harmony)
library(openxlsx)

setwd("/Users/Shared/Data/runAnalysis/singleSell/230304_Angio_2")
dir.create("./04.seurat")
cellrangers=dir(paste0(getwd(),"/03.outs_cellranger"))

a. 여러 샘플 동시에 돌리는 경우를 위한 ID insertion

for (i in 1:length(cellrangers))
  {data.i = Read10X(data.dir=paste0(getwd(),"/03.outs_cellranger/", cellrangers[i])) ; 
  colnames(data.i) = paste0(cellrangers[i], ".", colnames(data.i)) ; 
  obj.i = CreateSeuratObject(counts= data.i, project="angio_org", min.cells=3, min.features=300) ;
  obj.i$orig.ident = cellrangers[i] ; obj.i[["percent.mt"]] = PercentageFeatureSet(obj.i, "^MT-") ;
  cat(paste0("i = ",i," | ", cellrangers[i], "\n")) ; 
  if(i==1){obj_ori = obj.i} else {obj_ori = merge(obj_ori, obj.i)}} ; 

## i = 1 | Angio_2

save(obj_ori, file = paste0(getwd(),"/04.seurat/","obj_ori.rda")) 

01. Import된 data 확인

obj_ori

## An object of class Seurat 
## 18417 features across 691 samples within 1 assay 
## Active assay: RNA (18417 features, 0 variable features)

head(obj_ori@meta.data)

##                            orig.ident nCount_RNA nFeature_RNA percent.mt
## Angio_2.AAACGCTGTAGCACGA-1    Angio_2        620          311   1.451613
## Angio_2.AAACGCTGTCTCGCGA-1    Angio_2        974          378   1.334702
## Angio_2.AAAGAACAGATGCTGG-1    Angio_2     158981         9984   3.875935
## Angio_2.AAAGGATAGTACAGAT-1    Angio_2      31468         6188   4.595144
## Angio_2.AAAGGGCTCATGGCCG-1    Angio_2       6069         2152   6.986324
## Angio_2.AAAGGTAAGAGTTGAT-1    Angio_2     138819         9771   5.783790

table(obj_ori@meta.data$orig.ident)

## 
## Angio_2 
##     691

18417 features ; gene 691 samples ; cell

a. Data back-up & Violin plot으로 data preview

obj = subset(obj_ori, subset = nFeature_RNA > 200 & percent.mt < 20)
Idents(obj) = "orig.ident"
VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"))

02. Nomalization

obj <- NormalizeData(obj, normalization.method = "LogNormalize", scale.factor = 10000)
obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 2000)

a. Nomalization 후 상위 10개 표기해서 Variable plot

variable_plot <- VariableFeaturePlot(obj)
top10 <- head(VariableFeatures(obj), 10)
variable_plot.Top10 <- LabelPoints(plot = variable_plot, points = top10, repel = TRUE)

## When using repel, set xnudge and ynudge to 0 for optimal results

variable_plot.Top10

## Warning: Transformation introduced infinite values in continuous x-axis

03. PCA plot

## in case of removing unwanted sources of variation like MT contamination or cell cycling, regress out such heterogeneity
## obj = ScaleData (obj, features=rownames(obj), vars.to.regress="percent.mt") ;
obj = ScaleData(obj, features = rownames(obj)) ; 

## Centering and scaling data matrix

obj = RunPCA(obj, features=VariableFeatures (object=obj), npcs=100) ; 

## PC_ 1 
## Positive:  FBLN5, HYOU1, COL4A2, HTRA1, VASN, CXXC5, DNAJB9, ACSL3, AEBP1, CMKLR1 
##     MYADM, SGCD, LIPA, RRAS2, ALDH1L2, CD99, CPE, OAF, GLIS2, LOXL2 
##     FMOD, UNC5B, DUSP6, COL4A1, STAT1, SLC31A1, HERPUD1, TIMP3, HSPA5, RFLNB 
## Negative:  MT2A, MT-ND6, CCL5, MT-ATP6, NMB, MT-ND3, LINC02154, CXCL8, TFPI2, G0S2 
##     MT-CO1, MYEOV, CCDC151, MT-CYB, UBD, CCL3L1, EDN1, S100Z, AC090409.1, HSPA6 
##     EPGN, NLGN4X, AL360270.1, TMEM178B, AL022068.1, IL13RA2, HLA-DRB1, AC024597.1, AL355304.1, LINC02376 
## PC_ 2 
## Positive:  SLC22A3, RAMP1, IGFBP3, RPS4Y1, RFLNB, SPOCK1, LPL, EMILIN2, LMO2, APOD 
##     MGP, UNC5B, CXCL13, POSTN, EDNRA, FAM43B, TRIB3, ADM2, ITGA11, ATRNL1 
##     SYT7, SPON1, MFAP5, EDN3, IGF1, S100A4, CADM3, AL513283.1, EGFL6, MAF 
## Negative:  ARHGAP20, DHRS3, B3GALNT1, USP53, SCG2, CPNE7, SLC16A6, MPV17L, COL8A1, TGM2 
##     PGM2L1, DKK3, HOXB5, PCLO, RNF207, LSAMP, PDE5A, NCAM1, PCOLCE2, TRPV2 
##     TCEAL7, WNT2B, PLAT, NRP2, DGKH, SEC14L1, OLFML1, MTSS1, FNDC5, SRGN 
## PC_ 3 
## Positive:  IGFBP6, TIMP1, OLFM2, TXN, S100A6, PCBD1, CTSK, MT2A, CYBA, SRGN 
##     APOE, SDC1, FXYD5, STMN1, TCEAL7, ARL4C, CD68, PLPP2, PTGES, CYGB 
##     DHRS3, CSTB, APOD, CD99, DKK3, CAV1, JPT1, MAFB, OLFML1, OMG 
## Negative:  GOLGA8B, GOLGA8A, MEG8, AC058791.1, CCNL1, AC090673.1, AL354733.3, AC016831.1, CCDC14, MIR222HG 
##     MPV17L, NOTCH2NL, DGKH, KCNQ1OT1, LINC00624, ADAMTS6, SLC4A7, AC092807.3, NABP1, ART4 
##     PRUNE2, NEAT1, AFF1, SGK494, FSIP2, MIR29B2CHG, AC089983.1, EID3, RNF207, AKAP6 
## PC_ 4 
## Positive:  SOX2, DIO2, COL11A1, CMKLR1, NEGR1, OLFML1, PLXDC2, BOC, FOSB, B3GALNT1 
##     DDIT4, COL15A1, SYNE2, COL12A1, EDIL3, TXNIP, OLFML3, OMD, SOSTDC1, ITGA10 
##     DKK3, DHRS3, SERPINI1, IL11RA, MECOM, ADGRL3, COL14A1, SCN7A, EYA2, ARL4C 
## Negative:  SFRP2, IL32, MMP1, LAMA1, BIRC3, SOD2, TDO2, SLC12A8, AC004988.1, TNFAIP8L3 
##     GFRA1, HGF, DPP4, CPXM2, CA12, SNTG1, EPS8L2, HRK, SBSN, VEGFC 
##     KYNU, IFIT3, AKR1C2, CHI3L1, CD44, MEDAG, XG, AKR1C1, POU2F2, VNN1 
## PC_ 5 
## Positive:  SCARA5, AC002546.1, GPX3, SMOC2, EBF2, CILP, CIDEC, ENPP5, MEOX2, C1QL1 
##     PDGFA, NEGR1, SHISA9, PRG4, EFHD1, EYA2, SOSTDC1, OMD, INMT, CD4 
##     STC2, FSTL5, AFF2, BGN, TAGLN, ARRDC2, CX3CR1, NEURL1B, C7, NDUFA4L2 
## Negative:  ATRNL1, EMILIN2, GALNT9, SLC22A3, COL26A1, CRMP1, KIAA1549L, HS6ST2, SIPA1L2, CADM3 
##     MMP27, FXYD6, CAPN6, RELN, PCSK1, ANK3, DSC1, RTN4R, ENTPD1, RNF165 
##     EGLN3, SYT7, NXPH3, POSTN, CNTN1, EXPH5, HOPX, MXRA5, CDH2, PLPPR4

DimPlot(obj, reduction = "pca") ;

plot(obj@reductions$pca@stdev)

04. Jack Straw plot

#fig.width=10, fig.height=5, eval = FALSE
obj = JackStraw(obj, num.replicate=60, dims=60) ; obj = ScoreJackStraw(obj, dims=1:60)

#png(file=paste0(getwd(),"/saving_plot.png"), width=1000, height=500, unit="px", bg="white")
#pdf(file=paste0(getwd(),"/saving_plot.pdf"), 20,10)
#JackStrawPlot(obj, dims=1:100)
#dev.off()           
JackStrawPlot(obj, dims=1:60)

## Warning: Removed 96860 rows containing missing values (`geom_point()`).

05. tsne & umap

obj = RunTSNE(obj, reduction="pca", dims=1:60) ; 
obj = RunUMAP(obj, reduction="pca", dims=1:60) ; 

## 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

## 22:36:10 UMAP embedding parameters a = 0.9922 b = 1.112

## 22:36:10 Read 658 rows and found 60 numeric columns

## 22:36:10 Using Annoy for neighbor search, n_neighbors = 30

## 22:36:10 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 22:36:10 Writing NN index file to temp file /var/folders/z3/9m4zqbk57cbbr6n1rjvmh3980000gp/T//RtmpWbcJrz/file2ec84f41e499
## 22:36:10 Searching Annoy index using 1 thread, search_k = 3000
## 22:36:10 Annoy recall = 100%
## 22:36:11 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 22:36:12 Initializing from normalized Laplacian + noise (using irlba)
## 22:36:12 Commencing optimization for 500 epochs, with 30056 positive edges
## 22:36:13 Optimization finished

DimPlot(obj, reduction = "tsne") ; 

DimPlot(obj, reduction = "umap")

06. Batch effect ; sample 1개이므로 error

correction with “Harmony”

# library(harmony);
# luadobj = RunHarmony(luadobj, group.by.vars="orig.ident") ;
# ElbowPlot(luadobj, reduction="harmony", ndims=100) ;
# luadobj = RunTSNE (luadobj, reduction="harmony", dims=1:60, seed.use=1234) ;
# luadobj = RunUMAP (luadobj, reduction="harmony", dims=1:60, seed.use=1234) ;
# DimPlot(luadobj, reduction = "tsne") ;
# DimPlot(luadobj, reduction = "umap") ;
# save (luadobj, file="luadobj_batch.rda")

07. Clustering

#obj = FindNeighbors (obj, reduction= "harmony", k.param=20, dims=1:60)
obj = FindNeighbors (obj, k.param=20, dims=1:60)

## Computing nearest neighbor graph

## Computing SNN

obj = FindClusters (obj, resolution=1)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 658
## Number of edges: 41505
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.5889
## Number of communities: 4
## Elapsed time: 0 seconds

DimPlot (obj, reduction="umap", group.by="seurat_clusters", pt.size=0.001, label=T)

DimPlot (obj, reduction="tsne", group.by="seurat_clusters", pt.size=0.001, label=T)

Maximum modularity in 10 random starts: 0.5889
Number of communities: 4

08. Sub Clustering

#MAST has good FDR control and is faster than DESeq2
##BiocManager::install("MAST")
obj.markers = FindAllMarkers(obj, only.pos=TRUE, min.pct=0, logfc.threshold=0, test.use= "MAST")

## Calculating cluster 0

## 
## Done!

## Combining coefficients and standard errors

## Calculating log-fold changes

## Calculating likelihood ratio tests

## Refitting on reduced model...

## 
## Done!

## Calculating cluster 1

## 
## Done!

## Combining coefficients and standard errors

## Calculating log-fold changes

## Calculating likelihood ratio tests

## Refitting on reduced model...

## 
## Done!

## Calculating cluster 2

## 
## Done!

## Combining coefficients and standard errors

## Calculating log-fold changes

## Calculating likelihood ratio tests

## Refitting on reduced model...

## 
## Done!

## Calculating cluster 3

## 
## Done!

## Combining coefficients and standard errors

## Calculating log-fold changes

## Calculating likelihood ratio tests

## Refitting on reduced model...

## 
## Done!

#obj.markers
table(obj.markers$cluster)

## 
##    0    1    2    3 
## 9513 4757 4832 8043

sum(table(obj.markers$cluster))

## [1] 27145

obj.markers.top20 = obj.markers %>% dplyr::group_by(cluster) %>% dplyr:: top_n(n = 20, wt=avg_log2FC) ;
mySeuratClusters=unique(obj.markers.top20$cluster) ;
for(c in 1:length(mySeuratClusters)){
  obj.markers.top20.c = data.frame(
    cluster=obj.markers.top20[obj.markers.top20$cluster %in% mySeuratClusters[c], "gene"]) ;
  colnames(obj.markers.top20.c) = mySeuratClusters[c] ;
  if (c == 1) {obj.markers.top20s = obj.markers.top20.c} else {
    obj.markers.top20s = cbind(obj.markers.top20s, obj.markers.top20.c)}} ;
obj.markers.top20s

##            0         1         2      3
## 1   ARHGAP20     RBM43  C15orf61 RPS4Y1
## 2      ABCA6      CKS2     SESN3  UNC5B
## 3       TGM2    PTPN11  MTRNR2L8 IGFBP3
## 4   B3GALNT1     WDR41     PSMB6  TRIB3
## 5      DHRS3   LRRC75A      NQO1 S100A4
## 6       CLTC   MMP24OS MTRNR2L12 PDGFRB
## 7     SEZ6L2      NBDY   MT-ND4L  PLAC9
## 8       KLF6     DDIT3    MT-ND6  MFAP5
## 9       SCG2     HSPE1   MT-ATP6    LPL
## 10      MRC2     PLCG2    MT-ND5  HTRA1
## 11       JUN     COX7B    MT-ND1   CD99
## 12      CTGF       TKT    MT-ND2  TIMP3
## 13    PLXDC2     IFIT3     TFPI2 IGFBP5
## 14 POLR2J3.1    PHLDA2    MT-CO3    GSN
## 15    CYP1B1     RAB13    MT-CYB    MGP
## 16     ARL4C     CCND1    MT-ND3  MXRA8
## 17     NEAT1      SOD2    MT-CO1   IFI6
## 18      DDX5      G0S2    MT-CO2  POSTN
## 19     SOCS3      CCL5    MT-ND4  MFAP4
## 20    MALAT1 LINC02154     CXCL8   APOD

xlsx <- createWorkbook("xlsx") ;
addWorksheet(xlsx, cellrangers) ;
writeDataTable(xlsx, cellrangers, obj.markers.top20s, colNames = TRUE, sep = ", ") ;
saveWorkbook(xlsx, file=paste0(getwd(),"/04.seurat/",cellrangers,"_clustering.xlsx"))

09. Discovery of cluster idenity (20 -> 50 -> 100)

nonImm.markers = c("EPCAM", "KRT18", "KRT19", "VWF", "CLDN4", "COL1A2") ; 
FeaturePlot(obj, features=nonImm.markers, reduction="umap")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: KRT19

FeaturePlot(obj, features=nonImm.markers, reduction="tsne")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: KRT19

Imm.markers = c("CD79A", "CD3E", "NKG7" , "CD68" , "LYZ", "LGALS2", "TPSAB1", "JCHAIN") ; 
FeaturePlot(obj, features=Imm.markers, reduction="umap")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: CD3E, NKG7, LYZ, LGALS2,
## TPSAB1

FeaturePlot(obj, features=Imm.markers, reduction="tsne")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: CD3E, NKG7, LYZ, LGALS2,
## TPSAB1

Proliferating.markers = c("MKI67", "STMN1", "TOP2A") ; 
FeaturePlot(obj, features=Proliferating.markers, reduction="umap")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: MKI67

FeaturePlot(obj, features=Proliferating.markers, reduction="tsne")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: MKI67

angio.markers = c("CD31", "CD34", "VEGF", "VEGFR", "PDGF", "PDGFR") ; 
FeaturePlot(obj, features=angio.markers, reduction="umap")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: CD31, VEGF, VEGFR, PDGF,
## PDGFR

FeaturePlot(obj, features=angio.markers, reduction="tsne")

## Warning in FetchData.Seurat(object = object, vars = c(dims, "ident", features),
## : The following requested variables were not found: CD31, VEGF, VEGFR, PDGF,
## PDGFR

  obj$annot_1st = "" ;
obj@meta.data[obj@meta.data$seurat_clusters %in% c(,), ]$annot_1st <- "nonImm" ; obj@meta.data[obj@meta.data$seurat_clusters %in% c(0:1,), ]$annot_1st <- "Imm" ; obj@meta.data[obj@meta.data$seurat_clusters %in% c(),
]$annot_1st <- "prolif" ;

Idents(obj) = "annot_1st" ;
DimPlot(obj, group.by = "annot_1st", reduction-"umap", label-T) ;
DimPlot(obj, group.by = "annot_1st", reduction-"tsne", label-T)

obj.nonImm = subset(obj, subset = seurat_clusters %in% c (,)) ;
obj.nonImm = NormalizeData(obj.nonImm, normalization.method = "LogNormalize", scale.factor = 10000) ;
obj.nonImm = FindVariableFeatures(obj.nonImm, selection.method = "vst", nfeatures=2000) ;
obj.nonImm = ScaleData(obj.nonImm, features = rownames(obj.nonImm)) ;
obj.nonImm = RunPCA(obj.nonImm, features = VariableFeatures(object=obj.nonImm), npcs=100) ;
obj.nonImm = RunHarmony(obj.nonImm, "orig.ident") ;
plot(obj.nonImm@reductions$harmony@stdev) ;

obj.nonImm = RunUMAP (obj.nonImm, reduction = "harmony", dims=1:40, seed.use=1234) ;
obj.nonImm = RunTSNE (obj.nonImm, reduction = "harmony", dims=1:40, seed.use=1234) ;
obj.nonImm = FindNeighbors (obj.nonImm, reduction = "harmony", dims=1:40) ;
obj.nonImm = FindClusters (obj.nonImm, resolution = 0.5) ;
DimPlot(obj.nonImm, reduction="map", group.by="seurat_clusters", pt.size=0.001, label=T) ;

Epithelial.markers = c("EPCAM", "CLDN4", "ELF3") ;
FeaturePlot(obj.nonImm, features=Epithelial.markers, reduction= "umap", ncol=3) ;
FeaturePlot(obj.nonImm, features=Epithelial.markers, reduction= "tsne", ncol=3) ;
Endothelial.markers = c("CDH5", "CLDN5", "RAMP2") ;
FeaturePlot(obj.nonImm, features-Endothelial.markers, reduction= "umap", ncol=3) ;
FeaturePlot(obj.nonImm, features-Endothelial.markers, reduction= "umap", ncol=3) ;
Mesenchymal.markers = c("ITLNI", "COL1A2", "IGFBP6")   ;
FeaturePlot(obj.nonImm, features = Mesenchymal.markers, reduction-"umap", ncol=3) ;
FeaturePlot(obj.nonImm, features = Mesenchymal.markers, reduction-"umap", ncol=3) 

obj.nonImm$annot_2nd = "" ;
obj.nonImm@meta.data[obj.nonImm@meta.data$seurat_clusters %in% c(,),]$annot_2nd <- "Epithelial" ;
obj.nonImm@meta.data[obj.nonImm@meta.data$seurat_clusters %in% c(,),]$annot_2nd <- "Endothelial" ;
obj.nonImm@meta.data[obj.nonImm@meta.data$seurat_clusters %in% c(,),]$annot_2nd <- "Mesenchymal" ;
Idents(obj.nonImm) = "annot_2nd" ;
DimPlot(obj.nonImm, group.by = "annot_2nd", reduction = "umap", label=T) ;

obj.Imm = subset (obj, subset = seurat_clusters %in% c(,)) ;
obj.Imm = NormalizeData(obj.Imm, normalization.method = "LogNormalize", scale.factor = 10000) ;
obj.Imm = FindVariableFeatures(obj.Imm, selection.method = "vst", nfeatures=2000)
obj.Imm = ScaleData(obj.Imm, features = rownames(obj.Imm)) ;
obj.Imm = RunPCA(obj.Imm, features = VariableFeatures(object=dobj.Imm), npcs=100) ;
obj.Imm = RunHarmony(obj.Imm, "orig.ident") ;
plot(obj.Imm@reductions$harmony@stdev)

obj.Imm = RunUMAP(obj.Imm, reduction = "harmony", dims=1:50, seed.use=1234) ;
obj.Imm = RunTSNE(obj.Imm, reduction = "harmony", dims=1:50, seed.use=1234) ;
obj.Imm = FindNeighbors(obj.Imm, reduction= "harmony", dims=1:50) ;
obj.Imm = FindClusters(obj.Imm, resolution=0.8) ;
DimPlot(obj.Imm, reduction="umap", group.by-"seurat_clusters", pt.size=0.001, label=T) ;
DimPlot(obj.Imm, reduction-"tsne", group.by="seurat_clusters", pt.size=0.001, label-T) ;

NK.markers = c("GNLY", "KLRD1", "KLRF1") ;
FeaturePlot (obj.Imm, features=NK.markers, reduction= "umap", ncol-3) ;
T_common.markers = c("CD2", "CD3D") ;
FeaturePlot (obj.Imm, features=T_common.markers, reduction="umap", ncol=3) ;
CD4.markers = c("CD4", "CD40LG") ;
FeaturePlot (obj.Imm, features=CD4.markers, reduction= "umap", nco1=3) ;
CD8.markers = c("CD8A", "CD8B") ;
FeaturePlot (obj.Imm, features=CD8.markers, reduction= "umap", ncol-3) ;
gdT.markers = c("TRDV2" , "TRGV9") ;
FeaturePlot (obj.Imm, features=gdT.markers, reduction="umap", ncol=3)
B.markers = c("CD79A", "MS4A1", "IGKC") ;
FeaturePlot (obj.Imm, features=B.markers, reduction-"umap", nco1=3) ;
DC.markers = c("LGALS2", "CPVL", "CD1C") ;
FeaturePlot (obj.Imm, features=DC.markers, reduction="umap", ncol-3) ;
MQ.markers = c("MARCO", "C1QA", "FABP4") ;
FeaturePlot (obj.Imm, features=MQ.markers, reduction="umap", ncol=3) ;
Mono.markers = c("G0S2", "S100A8", "FCN1");
FeaturePlot (obj.Imm, features=Mono.markers, reduction-"umap", nco1=3)
MAST.markers = c("TPSAB1", "CPA3", "GATA2") ;
FeaturePlot (obj.Imm, features=MAST.markers, reduction= "umap", ncol=3)
pDC.markers = c("JCHAIN", "IRF4", "TCL1A") ;
FeaturePlot (obj.Imm, features=pD.markers, reduction="umap", ncol=3)  

obj.Imm$annot_2nd = "" ;
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), ]$annot_2nd <- "NK" ; 
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "CD4" ;
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "CD8" ;
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "B" ; 
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "DC" ;
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "MQ" ;
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "Mono" ;
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "MAST" ;
obj.Imm@meta.data[obj.Imm@meta.data$seurat_clusters %in% c(,), J$annot_2nd <- "pDC" ;

Idents(obi.Imm) = "annot_2nd" ;
DimPlot(obj.Imm, group.by = "annot_2nd", reduction="umap", label=T) ;

obj$annot_2nd = "" ;
obj@meta.data[obj@meta.data$annot_1st %in% "prolif",]$annot_2nd <- "prolif" ;
obj@meta.data[rownames(obj@meta.data) %in% rownames(obj.nonImm@meta.data),]$ annot_2nd <- obj.nonImm@meta.data$annot_2nd ;
obj@meta.data[rownames(obj@meta.data) %in% rownames(obj.Imm@meta.data),]$ annot_2nd <- obj.Imm@meta.data$annot_2nd ;