}
```
```{r}
silh_fun <- function(matr, label, nGenes = 1000, nPCs = 3){
t_matr <- t(matr)
HVG <- BrenneckeGetVariableGenes(t_matr, suppress.plot = T, fdr = 2)
matr <- matr[,HVG$Gene[1:nGenes]]
pca.data <- prcomp(matr, rank. = nPCs, center=TRUE)
dd <- dist(pca.data$x[, seq_len(nPCs)])
tibble(
label = unique(label$label)
) %>%
dplyr::mutate(num = 1:nrow(.)) -> tmp
label <- label %>% dplyr::left_join(tmp, by = c("label"))
summary(silhouette(label$num, dd))
}
```
```{r}
SEplot <- function(.x, point_size = 1, geom_line = T, p.adj = F, cutoff = 0.05){
if(isFALSE(p.adj)){
if(geom_line){
.x %>%
ggplot(aes(mean.expr, entropy)) +
geom_point(colour = '#1E90FF', size = point_size) +
geom_line(aes(mean.expr, fit), lwd = 0.7) +
theme_bw() +
theme(
axis.title = element_text(size = 15,color="black"),
axis.text = element_text(size = 15,color="black"),
legend.title = element_text(size = 0),
legend.text = element_text(size = 0),
axis.text.y = element_text(color="black"),
axis.text.x = element_text(color="black")
) +
labs(
x = "log(mean expression)",
y = "expression entropy"
) -> p
}
else{
.x %>%
ggplot(aes(mean.expr, entropy)) +
geom_point(colour = '#1E90FF', size = point_size) +
#geom_line(aes(mean.expr, fit), lwd = 0.7) +
theme_bw() +
theme(
axis.title = element_text(size = 15,color="black"),
axis.text = element_text(size = 15,color="black"),
legend.title = element_text(size = 0),
legend.text = element_text(size = 0),
axis.text.y = element_text(color="black"),
axis.text.x = element_text(color="black")
) +
labs(
x = "log(mean expression)",
y = "expression entropy"
) -> p
}
}
if(isTRUE(p.adj)){
.x <- .x %>% dplyr::mutate(sig = ifelse(p.adj <= cutoff, 1, 0))
if(geom_line){
.x %>%
ggplot(aes(mean.expr, entropy)) +
geom_point(aes(colour = factor(sig)), size = point_size) +
geom_line(aes(mean.expr, fit), lwd = 0.7) +
scale_color_manual(values = c("#1E90FF", "red")) +
theme_bw() +
theme(
legend.position = "none",
axis.title = element_text(size = 15,color="black"),
axis.text = element_text(size = 15,color="black"),
legend.title = element_text(size = 0),
legend.text = element_text(size = 0),
axis.text.y = element_text(color="black"),
axis.text.x = element_text(color="black")
) +
labs(
x = "log(mean expression)",
y = "expression entropy"
) -> p
}
else{
.x %>%
ggplot(aes(mean.expr, entropy)) +
geom_point(aes(colour = factor(sig)), size = point_size) +
#geom_line(aes(mean.expr, fit), lwd = 0.7) +
scale_color_manual(values = c("#1E90FF", "red")) +
theme_bw() +
theme(
legend.position = "none",
axis.title = element_text(size = 15,color="black"),
axis.text = element_text(size = 15,color="black"),
legend.title = element_text(size = 0),
legend.text = element_text(size = 0),
axis.text.y = element_text(color="black"),
axis.text.x = element_text(color="black")
) +
labs(
x = "log(mean expression)",
y = "expression entropy"
) -> p
}
}
return(p)
}
```
#rogue calculation for kmeans result
```{r}
kmeans_rouge <- function(matr, info){
tibble(Cluster = unique(info$Cluster)) %>%
dplyr::mutate(
matr = purrr::map(
.x = Cluster,
.f = function(.x){
info %>%
dplyr::filter(Cluster == .x) %>%
dplyr::pull(Barcode) -> barcode
tmp_matr <- matr[barcode,]
tmp_matr[is.na(tmp_matr)] <- 0
return(tmp_matr)
}
)
) -> tmp_mda
tmp_mda %>%
dplyr::mutate(
ent = purrr::map(
.x = matr,
.f = function(.x){
tmp <- SE_fun(expr = .x, window = 1, span = 0.2)
return(tmp)
}
)
) -> tmp_mda
tmp_mda <- tmp_mda %>%
dplyr::mutate(rogue = 1-purrr::map_dbl(ent, cal_rogue)) %>%
dplyr::select(-ent, -matr)
return(tmp_mda)
}
```
```{r}
get_mean <- function(matr, info){
info <- info %>% dplyr::mutate(ID = 1:nrow(.))
n_row <- length(unique(info$Cluster))
mean_expr <- matrix(data = NA, nrow = ncol(matr), ncol = n_row)
n <- 0
for (i in unique(info$Cluster)) {
n <- n + 1
ID <- info %>% dplyr::filter(Cluster == i) %>% dplyr::pull(ID)
mean_expr[,n] <- colMeans(matr[ID,])
}
rownames(mean_expr) <- colnames(matr)
colnames(mean_expr) <- paste("Cluster_",unique(info$Cluster), sep = "")
row_sum <- as.numeric(rowMeans(mean_expr))
row_sum <- log(row_sum+1)
tmp_expr <- log(mean_expr+1)
var1 <- rowMeans(tmp_expr)
mean_expr <- as.data.frame(mean_expr) %>% dplyr::mutate(ds = row_sum-var1)
rownames(mean_expr) <- colnames(matr)
return(mean_expr)
}
```
#mix cell types
```{r}
rogue_of_mixtures <- function(matr, info, cell_typs, n = 2000, rep = 20){
info1 <- info %>% dplyr::filter(label == cell_typs[1])
info2 <- info %>% dplyr::filter(label == cell_typs[2])
tibble(ratio = c(0,1/51,1/21,1/11,1/2)) %>%
dplyr::mutate(
rogue = purrr::map(
.x = ratio,
.f = function(.x){
num1 = ceiling(n*.x)
num2 = n-num1
print(num1)
tibble(Rep = 1:rep) %>%
dplyr::mutate(
rogue = purrr::map_dbl(
.x = Rep,
.f = function(.x){
index1 <- info1 %>% dplyr::sample_n(num1) %>% dplyr::pull(ID)
index2 <- info2 %>% dplyr::sample_n(num2) %>% dplyr::pull(ID)
index <- c(index1, index2)
expr <- matr[index,]
ent <- SE_fun(expr, span = 0.1)
rogue <- cal_rogue(ent)
return(rogue)
}
)
) -> rogue_rep
}
)
) -> tmp1
tibble(ratio = c(0,1/51,1/21,1/11)) %>%
dplyr::mutate(
rogue = purrr::map(
.x = ratio,
.f = function(.x){
num2 = ceiling(n*.x) #### cell number
num1 = n-num2  ####
print(num1)
tibble(Rep = 1:rep) %>%
dplyr::mutate(
rogue = purrr::map_dbl(
.x = Rep,
.f = function(.x){
index1 <- info1 %>% dplyr::sample_n(num1) %>% dplyr::pull(ID)
index2 <- info2 %>% dplyr::sample_n(num2) %>% dplyr::pull(ID)
index <- c(index1, index2)
expr <- matr[index,]
ent <- SE_fun(expr, span = 0.1)
rogue <- cal_rogue(ent)
return(rogue)
}
)
) -> rogue_rep
}
)
) -> tmp2
tmp2 <- tmp2[c(4:1),]
tmp <- tmp1 %>% dplyr::bind_rows(tmp2)
prop <- c("0:1","1:50","1:20","1:10","1:1","10:1","20:1","50:1","1:0")
tmp <- tmp %>% dplyr::mutate(prop = prop)
tmp <- tmp %>%
dplyr::mutate(
rogue = purrr::map2(
.x = rogue,
.y = prop,
.f = function(.x, .y){
.x %>% dplyr::mutate(prop = .y)
}
)
)
tmp <- tmp %>% dplyr::mutate(ID = 1:nrow(.))
tmp <- tmp %>%
dplyr::mutate(
rogue = purrr::map2(
.x = rogue,
.y = ID,
.f = function(.x, .y){
.x %>% dplyr::mutate(ID = .y)
}
)
)
return(tmp)
}
si_of_mixtures <- function(pca_da, info, cell_typs, n = 2000, rep = 20){
info1 <- info %>% dplyr::filter(label == cell_typs[1])
info2 <- info %>% dplyr::filter(label == cell_typs[2])
info_sub <- info %>% dplyr::filter(!(label %in% cell_typs))
tibble(ratio = c(0,1/21,1/11,1/4,1/2)) %>%
dplyr::mutate(
si = purrr::map(
.x = ratio,
.f = function(.x){
num1 = ceiling(n*.x)
num2 = n-num1
print(num1)
tibble(Rep = 1:rep) %>%
dplyr::mutate(
rogue = purrr::map_dbl(
.x = Rep,
.f = function(.x){
index1 <- info1 %>% dplyr::sample_n(num1) %>% dplyr::pull(ID)
index2 <- info2 %>% dplyr::sample_n(num2) %>% dplyr::pull(ID)
index <- c(index1, index2)
pda <- pca_da[c(index, info_sub$ID),]
info_tmp <- info[c(index, info_sub$ID),]
rownames(pda) <- 1:nrow(pda)
tibble(
label = unique(info_tmp$label)
) %>%
dplyr::mutate(num = 1:nrow(.)) -> tmp
info_tmp <- info_tmp %>% dplyr::left_join(tmp, by = c("label"))
dd <- dist(pda)
si <- summary(silhouette(info_tmp$num, dd))
return(as.numeric(si$clus.avg.widths[1]))
}
)
)
}
)
) -> tmp1
tibble(ratio = c(0,1/21,1/11,1/4)) %>%
dplyr::mutate(
si = purrr::map(
.x = ratio,
.f = function(.x){
num2 = ceiling(n*.x)
num1 = n-num2
print(num1)
tibble(Rep = 1:rep) %>%
dplyr::mutate(
rogue = purrr::map_dbl(
.x = Rep,
.f = function(.x){
index1 <- info1 %>% dplyr::sample_n(num1) %>% dplyr::pull(ID)
index2 <- info2 %>% dplyr::sample_n(num2) %>% dplyr::pull(ID)
index <- c(index1, index2)
pda <- pca_da[c(index, info_sub$ID),]
info_tmp <- info[c(index, info_sub$ID),]
rownames(pda) <- 1:nrow(pda)
tibble(
label = unique(info_tmp$label)
) %>%
dplyr::mutate(num = 1:nrow(.)) -> tmp
info_tmp <- info_tmp %>% dplyr::left_join(tmp, by = c("label"))
dd <- dist(pda)
si <- summary(silhouette(info_tmp$num, dd))
return(as.numeric(si$clus.avg.widths[1]))
}
)
)
}
)
) -> tmp2
tmp2 <- tmp2[c(4:1),]
tmp <- tmp1 %>% dplyr::bind_rows(tmp2)
prop <- c("0:1","1:20","1:10","1:3","1:1","3:1","10:1","20:1","1:0")
tmp <- tmp %>% dplyr::mutate(prop = prop)
tmp <- tmp %>%
dplyr::mutate(
si = purrr::map2(
.x = si,
.y = prop,
.f = function(.x, .y){
.x %>% dplyr::mutate(prop = .y)
}
)
)
tmp <- tmp %>% dplyr::mutate(ID = 1:nrow(.))
tmp <- tmp %>%
dplyr::mutate(
si = purrr::map2(
.x = si,
.y = ID,
.f = function(.x, .y){
.x %>% dplyr::mutate(ID = .y)
}
)
)
return(tmp)
}
```
```{r}
get_hd <- function(.x, .y, cell_type){
colnames(.y) <- colnames(.x)
.x <- .x %>%
dplyr::mutate(
mean.rogue = purrr::map_dbl(
.x = rogue,
.f = function(.x){
mean(.x$rogue)
}
)
)
.y <- .y %>%
dplyr::mutate(
mean.rogue = purrr::map_dbl(
.x = rogue,
.f = function(.x){
mean(.x$rogue)
}
)
)
.x <- .x %>%
dplyr::mutate(hd1 = 100*(.x$mean.rogue[1]-mean.rogue)/.x$mean.rogue[1]) %>%
dplyr::mutate(hd2 = 100*(.x$mean.rogue[9]-mean.rogue)/.x$mean.rogue[9]) %>%
dplyr::mutate(method = "rogue")
.y <- .y %>%
dplyr::mutate(hd1 = 100*(.y$mean.rogue[1]-mean.rogue)/2) %>%
dplyr::mutate(hd2 = 100*(.y$mean.rogue[9]-mean.rogue)/2) %>%
dplyr::mutate(method = "Silhouette")
.x$prop[6:9] <- .x$prop[4:1]
.y$prop[6:9] <- .y$prop[4:1]
tmp1 <- .x[2:5,] %>%
dplyr::bind_rows(.y[2:5,]) %>%
dplyr::mutate(label = paste(cell_type[1],"-",cell_type[2], sep = "")) %>%
dplyr::select(-hd2)
tmp2 <- .x[5:8,] %>%
dplyr::bind_rows(.y[5:8,]) %>%
dplyr::mutate(label = paste(cell_type[2],"-",cell_type[1], sep = "")) %>%
dplyr::rename(hd1 = hd2)
tmp1 %>% dplyr::bind_rows(tmp2)
}
```
```{r}
matr_filter <- function(.x, min.cells = 10, min.genes = 0){
cell_count <- colSums(.x > 0, na.rm = T)
gene_count <- rowSums(.x > 0, na.rm = T)
lq1 <- cell_count < min.cells
lq2 <- gene_count < min.genes
return(.x[!lq2, !lq1])
}
matr_toli <- function(ent, expr, n = 10, span = 0.1, r = 0.01){
sig.gene <- ent %>% dplyr::filter(p.adj < 0.05) %>% dplyr::pull(Gene)
ng <- length(sig.gene)
expr <- expr[,sig.gene]
mean.v <- c()
entr.v <- c()
for (i in 1:ng) {
.x <- as.numeric(expr[,i])
.x <- base::sort(.x, decreasing = T)
.x <- .x[-c(1:n)]
mean.v[i] <- log(mean(.x)+r)
entr.v[i] <- mean(log(.x+1))
}
mean.cut <- min(ent$mean.expr)
ent$mean.expr[1:ng] <- mean.v
ent$entropy[1:ng] <- entr.v
ent <- ent %>% dplyr::select(-p.adj) %>% dplyr::filter(mean.expr > mean.cut)
ent <- entropy_fit(ent, span = span)
return(ent)
}
```
```{r}
get_rogue_inc <- function(.x){
n_row <- nrow(.x)
.x$rogue[[1]] <- .x$rogue[[1]] %>% dplyr::mutate(anc = rogue)
tt <- .x
res1 <- .x$rogue[[1]]
a <- .x$cells[[1]]
for (i in 2:n_row) {
res2 <- .x$rogue[[i]]
b <- .x$cells[[i]]
res1 %>%
dplyr::mutate(
anc = purrr::map_dbl(
.x = cluster,
.f = function(.x){
tmp <- table(Var1 = b[a == .x]) %>% as.tibble()
tmp$n <- tmp$n/sum(tmp$n)
tmp %>%
dplyr::left_join(res2[,c(1,3)], by = c("Var1" = "cluster")) %>%
dplyr::mutate(rogue = n*rogue) %>%
dplyr::pull(rogue) %>%
sum()
}
)
) -> res1
tt$rogue[[i]] <- res1
}
average.rogue_ <- c()
for (i in 1:nrow(.x)) {
average.rogue_[i] <- mean(tt$rogue[[i]]$anc)
}
.x <- .x %>% dplyr::mutate(average.rogue = average.rogue_)
return(.x)
}
```
```{r}
bio.gene <- function(matr, batch, gene_num = 1500){
comb.ent <- SE_fun(matr)
comb.gene <- comb.ent$Gene[comb.ent$p.value < 0.01]
uni.batch <- unique(batch)
gene.list <- list()
for (i in 1:length(uni.batch)) {
tmp.expr <- matr[batch == uni.batch[i],]
tmp.expr <- matr_filter(tmp.expr, min.cells = 10)
tmp.ent <- SE_fun(tmp.expr)
gene.list[[i]] <- tmp.ent$Gene[tmp.ent$p.value < 0.01]
}
bio.gene <- gene.list[[1]]
for (i in 2:length(gene.list)) {
bio.gene <- unique(c(bio.gene, gene.list[[i]]))
}
bio.gene <- intersect(comb.gene, bio.gene)
return(bio.gene)
}
```
use_python("/home/heyao/data/tools/basic/anaconda3/bin/python")
anndata = import("anndata",convert=FALSE)
sc = import("scanpy.api",convert=FALSE)
np = import("numpy",convert=FALSE)
bbknn = import("bbknn", convert=FALSE)
package_version("reticulate")
package_version("reticulate")
package.version("reticulate")
package.version("ROGUE")
package.version("M3Drop")
package.version("GSVA")
ls()
anndata
isTRUE(T)
isTRUE(F)
ROGUE::matr.filter
?ROGUE::matr.filter
?ROGUE::matr.filter
?ROGUE::rogue
library(ROGUE)
rogue
?rogue
rogue
library(ROGUE)
?rogue
source('~/projects/04_SEmodel/06_R_package/ROGUE/R/ROGUE.R')
?rogue
rogue
library(roxygen2)
library(ROGUE)
rogue
?rogue
library(ROGUE)
library(roxygen2)
library(ROGUE)
library(ROGUE)
?rogue
?rogue.boxplot
library(ROGUE)
?rogue
library(ROGUE)
?SE_fun