.Rhistory

mu_genes <- sample(0:1000, ngene, replace = T)
bmu <- sample(1:4000, ngene, replace = T)
bmu <- bmu / 1000
var_genes <- mu_genes*bmu
plot(mu_genes, var_genes)
plot(sqrt(var_genes))
sd_genes <- sqrt(var_genes)
coreid <- nconnid[1]
# 1. start with matrix with no DE
mat0 <- matrix(0, nrow = ngene, ncol = nsample)
for (g in 1:ngene){
if (g %in% nconnid){
mat0[g,] <- rnorm(nsample, mean = mu_genes[coreid]*sample(1:50,1)/10, sd = sd_genes[coreid]*sample(1:50,1)/10)
} else {
mat0[g,] <- rnorm(nsample, mean = mu_genes[g], sd = sd_genes[g])
}
}
mat0[1:15,1:5]
mat0[1:nDE,ptime > cut1] <- 0
# 2. part of genes DE from cut point = cut1
DE_mu <- sample(500:1000, nDE, replace = T)
bmu2 <- sample(1:4000, nDE, replace = T)
bmu2 <- bmu2 / 1000
DE_sd <- sqrt(DE_mu*bmu2)
# assume they are all upregulated
mat1 <- matrix(0, nrow = ngene, ncol = nsample)
for (g in 1:nDE){
mat1[g,ptime > cut1] <- rnorm(sum(ptime > cut1), mean = DE_mu[g], sd = DE_sd[g])
}
mat2 <- mat0 + mat1
mat2[mat2 <0] <- 0
pheatmap(mat2[400:600, ], cluster_rows = F, cluster_cols = F)
rownames(mat2) <- paste("gene_", 1:ngene, sep = "")
colnames(mat2) <- paste("sample_", 1:nsample, sep = "")
# adjust for library size for the samples
boxplot(colSums(mat2))
library(SCDC)
mat.norm <- log2(getCPM0(mat2)*1e6+1)
boxplot(colSums(mat.norm))
pheatmap(mat.norm[400:600, 100:150], cluster_rows = F, cluster_cols = F)
simu.list <- list()
# samples by genes
for (l in 1:cutnum){
simu.list[[l]] <- t(mat.norm[,ptime >(l-1)/cutnum & ptime <= l/cutnum])
}
str(simu.list)
# begin modeling
test_k_div <- function(listk, alpha = 0.05){
t1 <- proc.time()
# add parameter selection procedure by BIC...
gnet <- JGL_k(Y=listk, penalty = "fused", lambda1 = 0.1, lambda2 = 0.1)
t2 <- proc.time()
message("JGL takes ", (t2-t1)[3], " seconds.")
K <- length(listk)
p <- ncol(gnet$theta[[1]])
theta1 <- gnet$theta[[K-1]]
theta2 <- gnet$theta[[K]]
wij <- matrix(nrow = p, ncol = p)
for (i in 1:p){
for (j in 1:p){
if (i <= j){
s1 <- (theta1[i,i]*theta1[j,j] + theta1[j,j]^2)/gnet$n[1]
s2 <- (theta2[i,i]*theta2[j,j] + theta2[j,j]^2)/gnet$n[2]
wij[i,j] <- (theta1[i,j] - theta2[i,j])/sqrt(s1+s2)
}
}
}
M <- max(wij^2, na.rm = T)
reject <- (M >= qgumbel(alpha) + 4*log(p) - log(log(p)))
pval <- pgumbel(M)
return(c(pval, reject))
}
# begin modeling
listk <- simu.list
t1 <- proc.time()
# add parameter selection procedure by BIC...
gnet <- JGL_k(Y=listk, penalty = "fused", lambda1 = 0.1, lambda2 = 0.1)
t2 <- proc.time()
install.packages('devtools')
devtools::install_github('xzhoulab/SPARK')
library(Rcpp)
sourceCpp("H:/994/papers/gene_networks/scripts/pairR2.cpp")
install.packages("pkgbuild") # pkgbuild is not available (for R version 3.5.0)
install.packages("devtools") # make sure you have the latest version from CRAN
library(devtools) # load package
# remove.packages("fansi")
devtools::install_github("r-lib/pkgbuild") # install updated version of pkgbuild from GitHub
install.packages("devtools")
library(devtools) # load package
# remove.packages("fansi")
devtools::install_github("r-lib/pkgbuild") # install updated version of pkgbuild from GitHub
# remove.packages("fansi")
devtools::install_github("r-lib/pkgbuild") # install updated version of pkgbuild from GitHub
install.packages("glue")
library(devtools) # load package
# remove.packages("fansi")
devtools::install_github("r-lib/pkgbuild") # install updated version of pkgbuild from GitHub
library(pkgbuild) # load package
find_rtools() # should be TRUE, assuming you have Rtools 3.5
library(devtools)
assignInNamespace("version_info", c(devtools:::version_info, list("3.5" = list(version_min = "3.3.0", version_max = "99.99.99", path = "bin"))), "devtools")
library(Rcpp)
sourceCpp("H:/994/papers/gene_networks/scripts/pairR2.cpp")
library(Seurat)
if (!require("devtools")) {
install.packages("devtools")
}
devtools::install_github("ne1s0n/MTGOsc")
install.packages("dplyr")
if (!require("devtools")) {
install.packages("devtools")
}
devtools::install_github("ne1s0n/MTGOsc")
remove.packages(dplyr)
remove.packages("dplyr")
if (!require("devtools")) {
install.packages("devtools")
}
devtools::install_github("ne1s0n/MTGOsc")
Sys.getenv("R_LIBS_USER")
Sys.getenv("R_LIBS_USER")
if (!require("devtools")) {
install.packages("devtools")
}
devtools::install_github("ne1s0n/MTGOsc")
remove.packages("dplyr")
library(mosum)
mosum
library(SCDC)
setwd("H:/994/SCDC/issues/data")
library(SCDC)
## read in ExpressionSet object from Perou lab, and Tabula Muris
bulk10x <- readRDS("peroubulk10x_fvb34.rds")
qc.perou<- readRDS("qc_perou.rds")
qc.tmouse<- readRDS("qc_tmouse.rds")
## here, tree-guided two-level deconvolution is performed (takes around 2 min):
perou <- qc.perou$sc.eset.qc
tmouse <- qc.tmouse$sc.eset.qc
perou$metacluster2[perou$md_cluster %in% c( "immune")] <- "immune"
perou$metacluster2[perou$md_cluster %in% c("basal","luminal","fibroblast","endothelial")] <- "BaLuFibEndo"
perou.subcl <- SCDC_prop_subcl_marker(bulk.eset = bulk10x, sc.eset = perou, ct.varname = "md_cluster",
fl.varname = "metacluster2", sample = "subj", ct.sub = c("endothelial","fibroblast","immune","luminal","basal"),
ct.fl.sub = unique(perou$metacluster2), select.marker = T, LFC.lim = 5)
tmouse$metacluster2[tmouse$md_cluster %in% c("immune")] <- "immune"
tmouse$metacluster2[tmouse$md_cluster %in% c("endothelial", "fibroblast","luminal","basal")] <- "BaLuFibEndo"
tmouse.subcl <- SCDC_prop_subcl_marker(bulk.eset = bulk10x, sc.eset = tmouse, ct.varname = "md_cluster",
fl.varname = "metacluster2", sample = "subj", ct.sub = c("endothelial","fibroblast","immune","luminal","basal"),
ct.fl.sub = unique(tmouse$metacluster2), select.marker = T, LFC.lim = 5)
ens_subcl_perou10x <- SCDC_ENSEMBLE(bulk.eset = bulk10x, prop.input = list(tmouse.subcl = tmouse.subcl, perou.subcl = perou.subcl),
ct.sub = c("endothelial","fibroblast","immune","luminal","basal"), search.length = 0.01, grid.search = T)
head(perou.subcl$prop.est)
SCDC_ENSEMBLE
bulk.eset = bulk10x
sc.eset.list = NULL
ct.sub = c("endothelial","fibroblast","immune","luminal","basal")
grid.search = F
search.length = 0.05
iter.max = 2000
nu = 1e-04
epsilon = 0.001
truep = NULL
weight.basis = T
prop.input = list(tmouse.subcl = tmouse.subcl, perou.subcl = perou.subcl)
Transform_bisque = F
if (!is.null(prop.input)) {
message("Using user-input estimated proportion list ...")
prop.list <- prop.input
}
row.list <- sapply(1:length(prop.list), function(x) {
rownames(prop.list[[x]]$yhat)
})
gene.prop <- Reduce("intersect", row.list)
head(row.list[[1]])
rownames(prop.list[[1]]$yhat)
str(prop.list[[1]])
SCDC_prop_subcl_marker
#' @param weight.basis logical, use basis matrix adjusted by MVW, default is T.
#' @param select.marker logical, select marker genes to perform deconvolution in tree-guided steps. Default is T.
#' @param markers A set of marker gene that input manully to be used in deconvolution. If NULL, then
#' @param marker.varname variable name of cluster groups when selecting marker genes. If NULL, then use ct.varname.
#' @param allgenes.fl logical, use all genes in the first-level deconvolution
#' @param pseudocount.use a constant number used when selecting marker genes, default is 1.
#' @param LFC.lim a threshold of log fold change when selecting genes as input to perform Wilcoxon's test.
#' @param truep true cell-type proportions for bulk samples if known
#' @return Estimated proportion, basis matrix, predicted gene expression levels for bulk samples
#' @export
SCDC_prop_subcl_marker <- function (bulk.eset, sc.eset, ct.varname, fl.varname, sample,
ct.sub = NULL, ct.fl.sub, iter.max = 3000, nu = 1e-04, epsilon = 0.001,
weight.basis = T, truep = NULL, select.marker = T, markers = NULL,
marker.varname = NULL, allgenes.fl = F, pseudocount.use = 1,
LFC.lim = 0.5, ...)
{
if (is.null(ct.sub)) {
ct.sub <- unique(sc.eset@phenoData@data[, ct.varname])[!is.na(unique(sc.eset@phenoData@data[,
ct.varname]))]
}
ct.sub <- ct.sub[!is.na(ct.sub)]
ct.fl.sub <- ct.fl.sub[!is.na(ct.fl.sub)]
bulk.eset <- bulk.eset[rowSums(exprs(bulk.eset)) > 0, , drop = FALSE]
sc.eset <- sc.eset[, sc.eset@phenoData@data[, ct.varname] %in%
ct.sub]
sc.basis <- SCDC_basis(x = sc.eset, ct.sub = ct.sub, ct.varname = ct.varname,
sample = sample)
sc.fl.basis <- SCDC_basis(x = sc.eset, ct.sub = ct.fl.sub[!is.na(ct.fl.sub)],
ct.varname = fl.varname, sample = sample)
if (select.marker) {
if (is.null(marker.varname)) {
marker.varname <- ct.varname
}
countmat <- exprs(sc.eset)
ct.group <- sc.eset@phenoData@data[, marker.varname]
markers.wilcox <- NULL
for (u in 1:length(unique(ct.group))) {
ct.group.temp <- ct.group == unique(ct.group)[u]
group.1 <- apply(X = countmat[, ct.group.temp], MARGIN = 1,
FUN = function(x) log(x = mean(x = expm1(x = x)) +
pseudocount.use))
group.2 <- apply(X = countmat[, !ct.group.temp],
MARGIN = 1, FUN = function(x) log(x = mean(x = expm1(x = x)) +
pseudocount.use))
genes.diff <- rownames(sc.eset)[(group.1 - group.2) >
LFC.lim]
count.use <- countmat[rownames(sc.eset) %in% genes.diff,
]
p_val <- sapply(1:nrow(count.use), function(x) {
wilcox.test(count.use[x, ] ~ ct.group.temp)$p.value
})
p_val_adj <- p.adjust(p = p_val, method = "bonferroni",
n = nrow(count.use))
markers.temp <- rownames(count.use)[p_val_adj < 0.05]
markers.wilcox <- c(markers.wilcox, markers.temp)
}
markers <- unique(markers.wilcox)
message("Selected ", length(markers), " marker genes by Wilcoxon test...")
}
if (weight.basis) {
basis <- sc.basis$basis.mvw
basis.fl <- sc.fl.basis$basis.mvw
}
else {
basis <- sc.basis$basis
basis.fl <- sc.fl.basis$basis
}
if (!is.null(markers)) {
commongenes <- Reduce(intersect, list(rownames(basis),
rownames(bulk.eset), markers))
commongenes.fl <- Reduce(intersect, list(rownames(basis.fl),
rownames(bulk.eset), markers))
}
else {
commongenes <- intersect(rownames(basis), rownames(bulk.eset))
commongenes.fl <- intersect(rownames(basis.fl), rownames(bulk.eset))
if (length(commongenes) < 0.2 * min(dim(sc.eset)[1],
dim(bulk.eset)[1])) {
stop("Too few common genes!")
}
}
message(paste("Used", length(commongenes), "common genes for all cell types, \n",
"Used", length(commongenes.fl), "common genes for first level cell types..."))
basis.mvw <- basis[commongenes, ct.sub]
basis.mvw.fl <- basis.fl[commongenes.fl, ct.fl.sub]
xbulk0 <- getCPM0(exprs(bulk.eset)[commongenes, ])
xbulk <- as.matrix(xbulk0)
colnames(xbulk) <- colnames(bulk.eset)
xbulk1 <- getCPM0(exprs(bulk.eset)[commongenes.fl, ])
xbulk.fl <- as.matrix(xbulk1)
ALS.S <- sc.basis$sum.mat[ct.sub]
N.bulk <- ncol(bulk.eset)
valid.ct <- (colSums(is.na(basis.mvw)) == 0) & (!is.na(ALS.S))
ALS.S.fl <- sc.fl.basis$sum.mat[ct.fl.sub]
valid.ct.fl <- (colSums(is.na(basis.mvw.fl)) == 0) & (!is.na(ALS.S.fl))
if (sum(valid.ct) <= 1) {
stop("Not enough valid cell type!")
}
message(paste("Used", sum(valid.ct), "cell types in deconvolution...\n",
"Used", sum(valid.ct.fl), "first level cell types ..."))
basis.mvw <- basis.mvw[, valid.ct]
ALS.S <- ALS.S[valid.ct]
basis.mvw.fl <- basis.mvw.fl[, valid.ct.fl]
ALS.S.fl <- ALS.S[valid.ct.fl]
prop.est <- NULL
rsquared <- NULL
for (i in 1:N.bulk) {
xbulk.temp <- xbulk[, i]
message(paste(colnames(xbulk)[i], "has common genes",
sum(xbulk[, i] != 0), "..."))
if (allgenes.fl) {
markers.fl <- names(xbulk.temp)
}
else {
markers.fl <- Reduce(intersect, list(markers, names(xbulk.temp)))
}
lm <- nnls::nnls(A = basis.mvw.fl[markers.fl, ], b = xbulk.temp[markers.fl])
delta <- lm$residuals
wt.gene <- 1/(nu + delta^2)
x.wt <- xbulk.temp[markers.fl] * sqrt(wt.gene)
b.wt <- sweep(basis.mvw.fl[markers.fl, ], 1, sqrt(wt.gene),
"*")
lm.wt <- nnls::nnls(A = b.wt, b = x.wt)
prop.wt.fl <- lm.wt$x/sum(lm.wt$x)
delta <- lm.wt$residuals
for (iter in 1:iter.max) {
wt.gene <- 1/(nu + delta^2)
x.wt <- xbulk.temp[markers.fl] * sqrt(wt.gene)
b.wt <- sweep(basis.mvw.fl[markers.fl, ], 1, sqrt(wt.gene),
"*")
lm.wt <- nnls::nnls(A = b.wt, b = x.wt)
delta.new <- lm.wt$residuals
prop.wt.fl.new <- lm.wt$x/sum(lm.wt$x)
if (sum(abs(prop.wt.fl.new - prop.wt.fl)) < epsilon) {
prop.wt.fl <- prop.wt.fl.new
delta <- delta.new
message("WNNLS for First level clusters Converged at iteration ",
iter)
break
}
prop.wt.fl <- prop.wt.fl.new
delta <- delta.new
}
names(prop.wt.fl) <- colnames(basis.mvw.fl)
rt <- table(sc.eset@phenoData@data[, ct.varname], sc.eset@phenoData@data[,
fl.varname])
rt <- rt[, ct.fl.sub]
rt.list <- list()
prop.wt <- NULL
for (j in 1:ncol(rt)) {
rt.list[[j]] <- rownames(rt)[rt[, j] > 0]
names(rt.list)[j] <- colnames(rt)[j]
sub.cl <- rownames(rt)[rt[, j] > 0]
if (length(sub.cl) > 1 & prop.wt.fl[colnames(rt)[j]] >
0) {
if (is.null(dim(prop.wt.fl))) {
xbulk.j <- basis.mvw.fl[, j] * prop.wt.fl[j] +
(xbulk.temp - basis.mvw.fl %*% lm.wt$x) *
prop.wt.fl[j]
}
else {
xbulk.j <- basis.mvw.fl[, j] * prop.wt.fl[,
j] + (xbulk.temp - basis.mvw.fl %*% lm.wt$x) *
prop.wt.fl[, j]
}
markers.sl <- Reduce(intersect, list(markers,
rownames(xbulk.j)))
basis.sl <- basis.mvw[markers.sl, rownames(rt)[rt[,
j] > 0]]
lm.sl <- nnls::nnls(A = basis.sl, b = xbulk.j[markers.sl,
])
delta.sl <- lm.sl$residuals
wt.gene.sl <- 1/(nu + delta.sl^2)
x.wt.sl <- xbulk.j[markers.sl, ] * sqrt(wt.gene.sl)
b.wt.sl <- sweep(basis.sl, 1, sqrt(wt.gene.sl),
"*")
lm.wt.sl <- nnls::nnls(A = b.wt.sl, b = x.wt.sl)
prop.wt.sl <- lm.wt.sl$x/sum(lm.wt.sl$x)
delta.sl <- lm.wt.sl$residuals
for (iter in 1:iter.max) {
wt.gene.sl <- 1/(nu + delta.sl^2)
x.wt.sl <- xbulk.j[markers.sl, ] * sqrt(wt.gene.sl)
b.wt.sl <- sweep(basis.sl, 1, sqrt(wt.gene.sl),
"*")
lm.wt.sl <- nnls::nnls(A = b.wt.sl, b = x.wt.sl)
delta.sl.new <- lm.wt.sl$residuals
prop.wt.sl.new <- lm.wt.sl$x/sum(lm.wt.sl$x)
if (sum(abs(prop.wt.sl.new - prop.wt.sl)) <
epsilon) {
prop.wt.sl <- prop.wt.sl.new
delta.sl <- delta.sl.new
cat("WNNLS for Second level clusters",
rownames(rt)[rt[, j] > 0], "Converged at iteration ",
iter)
break
}
prop.wt.sl <- prop.wt.sl.new
delta.sl <- delta.sl.new
}
names(prop.wt.sl) <- sub.cl
prop.wt <- c(prop.wt, prop.wt.sl * prop.wt.fl[colnames(rt)[j]])
}
else if (length(sub.cl) == 1) {
prop.wt <- c(prop.wt, prop.wt.fl[colnames(rt)[j]])
}
else if (length(sub.cl) > 1 & prop.wt.fl[colnames(rt)[j]] ==
0) {
prop.wt.sl <- rep(0, length(sub.cl))
names(prop.wt.sl) <- sub.cl
prop.wt <- c(prop.wt, prop.wt.sl)
}
}
prop.est <- rbind(prop.est, prop.wt)
}
rownames(prop.est) <- colnames(bulk.eset)
peval <- NULL
if (!is.null(truep)) {
peval <- SCDC_peval(ptrue = truep, pest = prop.est, pest.names = c("SCDC"),
select.ct = ct.sub)
}
# calculate yhat after deconv
yhat <- sc.basis$basis.mvw %*% t(prop.est)[colnames(sc.basis$basis.mvw),]
return(list(prop.est = prop.est, prop.wt.fl = prop.wt.fl,
basis.mvw = basis.mvw, peval = peval, sc.basis = sc.basis,
sc.fl.basis = sc.fl.basis, yhat = yhat))
}
perou.subcl <- SCDC_prop_subcl_marker(bulk.eset = bulk10x, sc.eset = perou, ct.varname = "md_cluster",
fl.varname = "metacluster2", sample = "subj", ct.sub = c("endothelial","fibroblast","immune","luminal","basal"),
ct.fl.sub = unique(perou$metacluster2), select.marker = T, LFC.lim = 5)
bulk10x <- readRDS("peroubulk10x_fvb34.rds")
qc.perou<- readRDS("qc_perou.rds")
qc.tmouse<- readRDS("qc_tmouse.rds")
## you could try the SCDC_prop() function for deconvolution without tree-guided structure first,
## but based on biological knowledge and the single cell information we have now, the proportions estimation result
## might not reflect the proportions that are close to truth.
## here, tree-guided two-level deconvolution is performed (takes around 2 min):
perou <- qc.perou$sc.eset.qc
tmouse <- qc.tmouse$sc.eset.qc
perou$metacluster2[perou$md_cluster %in% c( "immune")] <- "immune"
perou$metacluster2[perou$md_cluster %in% c("basal","luminal","fibroblast","endothelial")] <- "BaLuFibEndo"
perou.subcl <- SCDC_prop_subcl_marker(bulk.eset = bulk10x, sc.eset = perou, ct.varname = "md_cluster",
fl.varname = "metacluster2", sample = "subj", ct.sub = c("endothelial","fibroblast","immune","luminal","basal"),
ct.fl.sub = unique(perou$metacluster2), select.marker = T, LFC.lim = 5)
library(Biobase)
perou.subcl <- SCDC_prop_subcl_marker(bulk.eset = bulk10x, sc.eset = perou, ct.varname = "md_cluster",
fl.varname = "metacluster2", sample = "subj", ct.sub = c("endothelial","fibroblast","immune","luminal","basal"),
ct.fl.sub = unique(perou$metacluster2), select.marker = T, LFC.lim = 5)
head(perou.subcl$yhat)
tmouse$metacluster2[tmouse$md_cluster %in% c("immune")] <- "immune"
tmouse$metacluster2[tmouse$md_cluster %in% c("endothelial", "fibroblast","luminal","basal")] <- "BaLuFibEndo"
tmouse.subcl <- SCDC_prop_subcl_marker(bulk.eset = bulk10x, sc.eset = tmouse, ct.varname = "md_cluster",
fl.varname = "metacluster2", sample = "subj", ct.sub = c("endothelial","fibroblast","immune","luminal","basal"),
ct.fl.sub = unique(tmouse$metacluster2), select.marker = T, LFC.lim = 5)
ens_subcl_perou10x <- SCDC_ENSEMBLE(bulk.eset = bulk10x, prop.input = list(tmouse.subcl = tmouse.subcl, perou.subcl = perou.subcl),
ct.sub = c("endothelial","fibroblast","immune","luminal","basal"), search.length = 0.01, grid.search = T)
ens_subcl_perou10x$w_table
setwd("H:/994/SCDC/package/SCDC2/SCDC") # THE OLD PATH HAS ERROR TO UPDATE TO GITHUB
#longleaf: ------  setwd("/pine/scr/m/e/meichen/SCDC/SCDC")
devtools::document()
setwd("..")
devtools::install("SCDC")
#### update the vignettes website by:
setwd("H:/994/SCDC/package/SCDC2/SCDC")
pkgdown::build_site()
3
getPearson <- function(dt1, dt2, ct, subj){
cor(c(dt1[subj,ct]), c(dt2[subj,ct]))
}
combo.10x <- lapply(ens_subcl_perou10x$prop.list, function(x){
x$prop.est[,c("endothelial","fibroblast","luminal","basal","immune")]
})
ens_res <- function(wt, proplist, subject = c("FVB3","FVB4")){
ens <- wt_prop(wt, proplist = proplist)
p.ens <- getPearson(perou_pseudo_all$truep, ens ,
ct = c("endothelial","fibroblast","luminal","basal","immune"), subj = subject)
return(list(prop = ens, p = p.ens))
}
library(ggplot2)
library(reshape2)
ct1 <- c("mediumorchid1","mediumpurple1","lightskyblue","seagreen1","yellow","tan1","azure3")
getENSplot <- function(enstable, combo, subject, figname, subcat = NULL){
if (!is.null(subcat)){
enstable <- enstable[subcat,]
}
wtres.list <- list()
for (i in 1:nrow(enstable)){
wtres.list[[i]] <- ens_res(enstable[i,1:2], combo)$prop
}
names(wtres.list) <- rownames(enstable)
pcor.list <- NULL
for (i in 1:nrow(enstable)){
pcor.list[[i]] <- ens_res(enstable[i,1:2], combo, subject = subject)$p
}
names(pcor.list) <- rownames(enstable)
est.all <- rbind(perou_pseudo_all$truep[,c("endothelial","fibroblast","luminal","basal","immune")],
do.call(rbind, combo),
do.call(rbind, wtres.list))
est.all.3 <- est.all[rownames(est.all) == subject,]
rownames(est.all.3) <- as.factor(c("Seurat", "Tabula Muris", "Perou",
paste("ENSEMBLE:",names(pcor.list))))
p_t <- getPearson(perou_pseudo_all$truep, combo[[1]],
ct = c("endothelial","fibroblast","luminal","basal","immune"), subj = c(subject))
p_p <- getPearson(perou_pseudo_all$truep, combo[[2]],
ct = c("endothelial","fibroblast","luminal","basal","immune"), subj = c(subject))
dat_text <- data.frame(label = c("",paste("R=", round(c(p_t, p_p, pcor.list),2), sep = "")),
Var1 = as.factor(c("Seurat", "Tabula Muris", "Perou",
paste("ENSEMBLE:",names(pcor.list)))),
value =rep(0.9, length(pcor.list)+3),
Var2 = rep("endothelial",length(pcor.list)+3))
meltdt <- melt(est.all.3)
P <- ggplot(meltdt, aes(x=Var1, y=value, fill = factor(Var2, levels = c("endothelial","fibroblast","luminal","basal","immune")))) +
geom_bar(stat = "identity") + xlab("") + ylab("Percentage") +
theme(axis.text.x = element_text(angle = 30, hjust = 1, size=10),
axis.text.y = element_text(size = 10),
text = element_text(size = 10),
plot.title = element_text(size=10, face = "bold"),
plot.margin=unit(c(1,1,-5,0), "mm"),
legend.position="top",legend.title = element_blank(),
legend.text = element_text(size=8),
legend.box.spacing = unit(0, "mm"),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+
scale_fill_manual(values=ct1) +
guides(fill = guide_legend(nrow = 2))+
# facet_grid(cols = vars(ref)) +
geom_text(data = dat_text,
mapping = aes(x=Var1, y=value,label = label),
hjust   = 0.5, vjust   = 0)
png(file.path(vpath,  "mammarygland.png"), res = 80, height = 300, width = 300)
print(P)
dev.off()
}
setwd("H:/994/SCDC/package/SCDC2/SCDC")
pkgdown::build_site()
setwd("H:/994/SCDC/package/SCDC2/SCDC")
pkgdown::build_site()