-
Notifications
You must be signed in to change notification settings - Fork 0
/
corr_spearman.R
66 lines (48 loc) · 2.87 KB
/
corr_spearman.R
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
##
## Calculate the Spearman's rank correlation coefficient for every PBMC cell type
##
## Dataset: 33K PBMC 10xGenomics
## Loaded data: data/seurat_33k_merged (created @ pbmc33k.R)
## Saved data: data/33k_correlation_spearman_celltype
##
library(Matrix)
library(scran)
library(Seurat)
## Data loading
# Load the clustered data.
load("data/seurat_33k_merged")
# Rename cluster IDs.
new_ids=c("Megakaryocytes","Monocytes",3,"NK_1","NK_2","B_Cells_Plasma","Dendritic_Cells","B_Cells_1","B_Cells_2","B_Cells_3","T_Cells_CD4+_1","T_Cells_CD8+_1","T_Cells_CD4+_2","T_Cells_CD8+_2","T_Cells_CD8+_3","T_Cells_CD8+_4")
current.cluster.ids=1:16
pbmc33k.merged@ident <- plyr::mapvalues(pbmc33k.merged@ident, from = current.cluster.ids, to = new_ids)
## Data preperation
# Filter genes with less than 10 percent non-zero values.
fraction.non.zero <- 0.05
pbmc33k.non.zero.counts <- rowSums(as.matrix(pbmc33k.merged@data) != 0)
genes.filtered <- pbmc33k.non.zero.counts > fraction.non.zero * ncol(pbmc33k.merged@data)
pbmc33k.filtered.counts <- pbmc33k.merged@data[genes.filtered,]
dim(pbmc33k.merged@data)
dim(pbmc33k.filtered.counts)
#TSNEPlot(pbmc33k.merged)
#fraction.cell.type = 1/2
# Subset the celltypes from filtered data
nk.cells <- pbmc33k.filtered.counts[,WhichCells(pbmc33k.merged, c("NK_1", "NK_2"))]
b.cells <- pbmc33k.filtered.counts[,WhichCells(pbmc33k.merged, c("B_Cells_1", "B_Cells_2", "B_Cells_3", "B_Cells_Plasma"))]
t.h.cells <- pbmc33k.filtered.counts[,WhichCells(pbmc33k.merged, c("T_Cells_CD4+_1", "T_Cells_CD4+_2"))]
t.c.cells <- pbmc33k.filtered.counts[,WhichCells(pbmc33k.merged, c("T_Cells_CD8+_1", "T_Cells_CD8+_2", "T_Cells_CD8+_3","T_Cells_CD8+_4"))]
monocytes <- pbmc33k.filtered.counts[,WhichCells(pbmc33k.merged, c("Monocytes"))]
## Calculate correlations
# When all sample are used algorithm fails
sample.count <- 1500
nk.cells.cor.pairs <- correlatePairs(as.matrix(nk.cells)[,1:sample.count])
nk.cells.cor.sign <- nk.cells.cor.pairs[nk.cells.cor.pairs$FDR < 0.05,]
b.cells.cor.pairs <- correlatePairs(as.matrix(b.cells)[,1:sample.count])
b.cells.cor.pairs.sign <- b.cells.cor.pairs[b.cells.cor.pairs$FDR < 0.05,]
t.h.cells.cor.pairs <- correlatePairs(as.matrix(t.h.cells)[,1:sample.count])
t.h.cells.cor.pairs.sign <- t.h.cells.cor.pairs[t.h.cells.cor.pairs$FDR < 0.05,]
t.c.cells.cor.pairs <- correlatePairs(as.matrix(t.c.cells)[,1:sample.count])
t.c.cells.cor.pairs.sign <- t.c.cells.cor.pairs[t.c.cells.cor.pairs$FDR < 0.05,]
monocytes.cor.pairs <- correlatePairs(as.matrix(monocytes)[,1:sample.count])
monocytes.cor.pairs.sign <- monocytes.cor.pairs[monocytes.cor.pairs$FDR < 0.05,]
## Save the correlations
save(nk.cells.cor.sign, nk.cells.cor.pairs, b.cells.cor.pairs.sign, b.cells.cor.pairs, t.c.cells.cor.pairs.sign, t.c.cells.cor.pairs, t.h.cells.cor.pairs.sign, t.h.cells.cor.pairs, monocytes.cor.pairs.sign, monocytes.cor.pairs, file = "data/33k_correlation_spearman_celltype_frac_0.05")