-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSection_5_step_1_run_lochness_calculation.R
93 lines (82 loc) · 3.52 KB
/
Section_5_step_1_run_lochness_calculation.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
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
library(Seurat)
library(Matrix)
library(irlba)
library(ggplot2)
library(dplyr)
library(plyr)
library(monocle3)
library(argparse)
library(FNN)
set.seed(2020)
safe_column_multiply = function(bmat, scale_vector) {
bmat@x <- bmat@x * rep.int(scale_vector, diff(bmat@p))
return(bmat)
}
wd = './'
sd = './'
# Load data
print('Loading data...')
pca_object_list = readRDS('../data/sub_trajectory_aligned_pca_coor.rds')
umap_object = readRDS('../data/sub_trajectory_umap_coor.rds')
pd = readRDS('../data/pd.rds')
mtx = readRDS('../data/count.rds')
pd$Background[pd$Background=="C57BL/6"]<-"C57BL6"
pd$Background[pd$Background=="BALB/C"]<-"BALBC"
pd$Background_Mutant = paste0(pd$Background,':',pd$Mutant)
full_seurat_obj = CreateSeuratObject(mtx,meta.data=pd)
full_seurat_obj[['umap']] = Seurat::CreateDimReducObject(embeddings=as.matrix(umap_object), key='UMAP_', assay='RNA')
##pool
for (i in unique(full_seurat_obj$Background_Mutant)) {
mt = sub(".*\\:", "", i)
bg = sub("\\:.*", "", i)
if (mt != 'WT') {
for (tj in names(pca_object_list)) {
print(paste0(bg,'_',mt,'_',tj))
seurat_obj = full_seurat_obj[,full_seurat_obj$main_trajectory==tj]
pca_object = pca_object_list[[tj]]
pca_dims = dim(pca_object)[2]
seurat_obj[['pca']] = Seurat::CreateDimReducObject(embeddings=as.matrix(pca_object), key='PC_', assay='RNA')
seurat_object = seurat_obj[,seurat_obj$Background_Mutant==i | seurat_obj$Mutant == "WT"]
kadj = round(0.5 * sqrt(ncol(seurat_object)))
####PCA space
seurat_object = seurat_object %>%Seurat::L2Dim(reduction='pca')
df = data.frame()
for (j in unique(seurat_object$RT_group)) {
mt_cells = colnames(seurat_object)[seurat_object$RT_group==j]
wt_cells = colnames(seurat_object)[seurat_object$RT_group!=j]
mt.mtx = as.matrix(t(as.data.frame(seurat_object@reductions$pca.l2@cell.embeddings)[mt_cells, ]))
wt.mtx = as.matrix(t(as.data.frame(seurat_object@reductions$pca.l2@cell.embeddings)[wt_cells, ]))
# Run KNN and assign labels based on majority of NNs
knn = get.knnx(t(wt.mtx), t(mt.mtx), k=kadj, algorithm="kd_tree")
knn$label <- mapvalues(knn$nn.index,from=seq(1, ncol(wt.mtx)),to=as.character(seurat_object[,wt_cells]$Mutant),warn_missing=FALSE)
# MT scores
mt_counts = as.data.frame(table(seurat_object[,wt_cells]$Mutant))
mutant_mask = ifelse(knn$label=="WT", 0, 1)
mutant_neighbors = rowSums(mutant_mask)
pca_mutant_score = as.data.frame((mutant_neighbors / kadj) / (1 - mt_counts$Freq[mt_counts$Var1=='WT'] / sum(mt_counts$Freq)) - 1)
rownames(pca_mutant_score) = mt_cells
colnames(pca_mutant_score) = c("raw_mt_score")
df = rbind(df,pca_mutant_score)
}
#plotting
df_umap = as.data.frame(seurat_object[['umap']]@cell.embeddings)
df_umap$raw_mt_score = df[rownames(df_umap),]
df_umap$main_trajectory = seurat_object$main_trajectory
saveRDS(df_umap,paste0(sd,bg,'_',mt,'_',tj,'_lochness_excludeself_pool.rds'))
}
}
}
df_list = list()
for (i in unique(full_seurat_obj$Background_Mutant)) {
mt = sub(".*\\:", "", i)
bg = sub("\\:.*", "", i)
if (mt != 'WT') {
for (tj in names(pca_object_list)) {
df_umap=readRDS(paste0(sd,bg,'_',mt,'_',tj,'_lochness_excludeself_pool.rds'))
df_umap$mt = mt
df_umap$bg = bg
df_umap$cell_id = rownames(df_umap)
df_list[[paste0(tj,i)]] = df_umap
}}}
dff = do.call(rbind,df_list)
saveRDS(dff,paste0(sd,'lochness_excludeself_pool.rds'))