-
Notifications
You must be signed in to change notification settings - Fork 0
/
Preprocessing_HumanAKI_snRNA-seq_step2
128 lines (108 loc) · 5.63 KB
/
Preprocessing_HumanAKI_snRNA-seq_step2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#Data integration after preprocessing of individual datasets
#Preprocessing / aggregation of control datasets were similarly performed before
#https://github.com/TheHumphreysLab/Multimodal_analysis_ADPKD/blob/main/1_RNA_individual_data.R
library(Seurat)
library(ggplot2)
library(sctransform)
library(harmony)
library(Rcpp)
library(dplyr)
set.seed(1234)
######### AKI dataset ##########
#AKI_1/2/3/4
AKI_1@meta.data[["orig.ident"]] <- "AKI1"
AKI_2@meta.data[["orig.ident"]] <- "AKI2"
AKI_3@meta.data[["orig.ident"]] <- "AKI3"
AKI_4@meta.data[["orig.ident"]] <- "AKI4"
# Integration (AKI dataset)
AKI.list <- list(AKI_1,AKI_2,AKI_3,AKI_4)
AKI.list <- lapply(X = AKI.list, FUN = SCTransform)
features <- SelectIntegrationFeatures(object.list = AKI.list, nfeatures = 3000)
AKI.list <- PrepSCTIntegration(object.list = AKI.list, anchor.features = features)
AKI.anchors <- FindIntegrationAnchors(object.list = AKI.list, normalization.method = "SCT",
anchor.features = features)
AKI <- IntegrateData(anchorset = AKI.anchors, normalization.method = "SCT")
DefaultAssay(AKI) <- "integrated"
# clustering1
AKI <- ScaleData(AKI, verbose = FALSE)
AKI <- RunPCA(AKI, npcs = 30, verbose = FALSE)
AKI <- RunUMAP(AKI, reduction = "pca", dims = 1:20)
AKI <- FindNeighbors(AKI, reduction = "pca", dims = 1:20)
AKI <- FindClusters(AKI, resolution = 0.4)
DimPlot(AKI, reduction = "umap", label = TRUE)
# annotation
Idents(AKI) <- "seurat_clusters"
new.cluster.ids <- c("TAL","PT","PT","CNT_PC","DCT","PT","FIB","ENDO2",
"PT_Injured","ICA","LowQC","ICB","LEUK1","LEUK2","PODO","ENDO1","PEC","Cycling")
names(new.cluster.ids) <- levels(AKI)
AKI <- RenameIdents(AKI, new.cluster.ids)
levels(AKI) <- c("PT","PT_Injured","PEC","TAL","DCT","CNT_PC","ICA","ICB","PODO",
"ENDO1","ENDO2","FIB","LEUK1","LEUK2","Cycling","LowQC")
#LowQC cluster does not have specific markers, and they enrich mitochondrial genes
#Remove LowQC & doublets (Amulet) from AKI dataset
AKI <- subset(AKI,idents = "LowQC",invert=T)
Idents(AKI) <- "doubletdata"
AKI <- subset(AKI,idents = "Singlet")
DefaultAssay(AKI) <- "integrated"
# clustering2
AKI <- ScaleData(AKI, verbose = FALSE)
AKI <- RunPCA(AKI, npcs = 30, verbose = FALSE)
AKI <- RunUMAP(AKI, reduction = "pca", dims = 1:20)
AKI <- FindNeighbors(AKI, reduction = "pca", dims = 1:20)
AKI <- FindClusters(AKI, resolution = 0.4)
DimPlot(AKI, reduction = "umap", label = TRUE)
#annotation
Idents(AKI) <- "seurat_clusters"
new.cluster.ids <- c("PT","TAL","PT","DCT","CNT_PC","PT","FIB2","ENDO2",
"PT_Injured","ICA","ICB","LYMP","MAC","PEC","ENDO1","PODO","Cycling","FIB1")
names(new.cluster.ids) <- levels(AKI)
AKI <- RenameIdents(AKI, new.cluster.ids)
levels(AKI) <- c("PT","PT_Injured","PEC","TAL","DCT","CNT_PC","ICA","ICB","PODO",
"ENDO1","ENDO2","FIB1","FIB2","LYMP","MAC","Cycling")
AKI@meta.data[["celltype"]] <- AKI@active.ident
FigS12A <- DimPlot(AKI, reduction = "umap", label = TRUE), FigS12A
DefaultAssay(AKI) <- "RNA"
levels(AKI) <- rev(levels(AKI))
features <- c("SLC34A1","LRP2","HAVCR1","CFH","SLC12A1","SLC12A3",
"SLC8A1","AQP2","SLC26A7","SLC26A4","NPHS2",
"EMCN","HECW2","GPM6A","FLRT2","PDGFRB","PIEZO2","PDGFA","C7"
,"PTPRC","CD96","THEMIS","CD163","MRC1","MKI67")
FigS12B<- DotPlot(AKI, features = features, cols = c("lightyellow","royalblue")) +
RotatedAxis() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) #FigS12B
levels(AKI) <- rev(levels(AKI))
#control aggregated dataset (n=5) was prepared previously
#https://github.com/TheHumphreysLab/Multimodal_analysis_ADPKD/blob/main/1_RNA_individual_data.R
cont <- readRDS("control_aggregated_dataset.rds")
cont@meta.data[["celltype"]] <- paste("Cont",cont@active.ident,sep = "-")
AKI@meta.data[["celltype"]] <- paste("AKI",AKI@active.ident,sep = "-")
#merging ADPKD and control data
rnaAggr <- merge(AKI,cont,add.cell.ids = c("AKI","Cont"))
rnaAggr@meta.data <- rnaAggr@meta.data[,c(1:6,21)]
Idents(rnaAggr) <- "orig.ident"
current.ids<-c("AKI1","AKI2","AKI3","AKI4","control1","control2","control3","control4","control5")
new.ids<-c('male','female','female','male','male','male','female','male','female')
rnaAggr@meta.data$gender <-plyr::mapvalues(rnaAggr@meta.data$orig.ident,from = current.ids,to=new.ids)
current.ids<-c("AKI1","AKI2","AKI3","AKI4","control1","control2","control3","control4","control5")
new.ids<-c('AKI','AKI','AKI','AKI','control','control','control','control','control')
rnaAggr@meta.data$disease <-plyr::mapvalues(rnaAggr@meta.data$orig.ident,from = current.ids,to=new.ids)
#Integration with Harmony on assay "RNA"
DefaultAssay(rnaAggr) <- "RNA"
rnaAggr <- NormalizeData(rnaAggr)
rnaAggr <- FindVariableFeatures(rnaAggr, selection.method = "vst", nfeatures = 3000)
rnaAggr <- ScaleData(rnaAggr, verbose = FALSE)
rnaAggr <- RunPCA(rnaAggr, verbose = FALSE)
rnaAggr <- RunHarmony(rnaAggr,group.by.vars = "orig.ident")
rnaAggr <- RunUMAP(rnaAggr, reduction = "harmony", dims = 1:20)
rnaAggr <- FindNeighbors(rnaAggr, reduction = "harmony", dims = 1:20)
rnaAggr <- FindClusters(rnaAggr, resolution = 0.5)
#Annotation
Idents(rnaAggr) <- "seurat_clusters"
new.cluster.ids <- c("PT","TAL2","DCT","CNT_PC","TAL1","ENDO","FIB","iPT","ICA","ICB","PEC","LEUK","PODO","PT")
names(new.cluster.ids) <- levels(rnaAggr)
rnaAggr <- RenameIdents(rnaAggr, new.cluster.ids)
levels(rnaAggr) <- c("PT","iPT","PEC","TAL1","TAL2","DCT","CNT_PC","ICA","ICB","PODO","ENDO","FIB","LEUK")
rnaAggr@meta.data[["celltype_all"]] <- rnaAggr@active.ident
#sessionInfo()
#[1] dplyr_1.0.5 harmony_0.1.0 Rcpp_1.0.8.3 sctransform_0.3.2
#[5] ggplot2_3.3.3 SeuratObject_4.0.0 Seurat_4.0.0