#
import_network_data=function(file=NULL,type="protein",species="human"){
if(species=="human"){
if((type=="gene" || type=="protein") && is.null(file)){
x <- readRDS(system.file("extdata/corum.rds",package = "OmicsEV"))
}else if(type=="phosphopeptide"){
x <- readRDS(system.file("extdata/psp.rds",package = "OmicsEV"))
}
}else if(species=="drosophila"){
if((type=="gene" || type=="protein") && is.null(file)){
x <- readRDS(system.file("extdata/flybase_ppi.rds",package = "OmicsEV"))
}
}else{
cat("species:",species,"\n")
stop("The species is not supported!")
x <- NULL
}
return(x)
}
import_network_data_from_table=function(file){
a <- readr::read_tsv(file)
x <- a$complex
n_p <- sapply(x,function(y){
a=unique(unlist(strsplit(split=";",x=y)));
length(a)})
# filter complex only having one protein
x <- x[n_p>=2]
x <- sapply(x, function(y){
a=unique(unlist(strsplit(split=";",x=y)));
a=sort(a);paste(a,sep=";",collapse = ";")}) %>%
unique
aa <- sapply(x, function(y){
t(combn(unlist(strsplit(split = ";",x=y)),2))})
names(aa) <- NULL
aa <- do.call(rbind,aa)
res <- list()
res$complex <- apply(aa,1,function(xx){
sort(str_replace_all(xx,pattern=" ",replacement=""))}) %>%
t %>%
tibble::as_tibble() %>%
distinct() %>%
filter(V1!="",V2!="")
res$complex <- res$complex %>% filter(V1!=V2)
## random complex
rnd <- c(res$complex$V1,res$complex$V2) %>%
unique %>%
sort %>%
combn(m=2) %>%
t %>%
tibble::as_tibble()
yyy <- paste(res$complex$V1,res$complex$V2,sep="|")
rnd_p <- paste(rnd$V1,rnd$V2,sep="|")
rnd <- rnd[ !(rnd_p %in% yyy), ]
res$rnd_complex <- rnd
names(res) <- c("network","rnd_network")
return(res)
}
import_data_from_corum=function(file){
a <- read.csv(file,stringsAsFactors = FALSE)
x <- a$subunits.Gene.name.
n_p <- sapply(x,function(y){
a=unique(unlist(strsplit(split=";",x=y)));
length(a)})
# filter complex only having one protein
x <- x[n_p>=2]
x <- sapply(x, function(y){
a=unique(unlist(strsplit(split=";",x=y)));
a=sort(a);paste(a,sep=";",collapse = ";")}) %>%
unique
aa <- sapply(x, function(y){
t(combn(unlist(strsplit(split = ";",x=y)),2))})
names(aa) <- NULL
aa <- do.call(rbind,aa)
res <- list()
res$complex <- apply(aa,1,function(xx){
sort(str_replace_all(xx,pattern=" ",replacement=""))}) %>%
t %>%
tibble::as_tibble() %>%
distinct() %>%
filter(V1!="",V2!="")
res$complex <- res$complex %>% filter(V1!=V2)
## random complex
rnd <- c(res$complex$V1,res$complex$V2) %>%
unique %>%
sort %>%
combn(m=2) %>%
t %>%
tibble::as_tibble()
yyy <- paste(res$complex$V1,res$complex$V2,sep="|")
rnd_p <- paste(rnd$V1,rnd$V2,sep="|")
rnd <- rnd[ !(rnd_p %in% yyy), ]
res$rnd_complex <- rnd
names(res) <- c("network","rnd_network")
return(res)
}
# wget http://www.droidb.org/data/DroID_v2018_08/flybase_ppi.txt
# file = "inst/extdata/flybase_ppi.txt"
# import_data_from_flybase_ppi(file)
import_data_from_flybase_ppi=function(file){
a <- read.delim(file,stringsAsFactors = FALSE)
x <- a %>% select(SYMBOL1,SYMBOL2)
names(x) <- c("V1","V2")
x <- x %>% filter(V1!=V2) %>% distinct() %>%
apply(1,function(y) {sort(y)}) %>%
t %>%
as.data.frame(stringsAsFactors=FALSE) %>%
distinct()
all_genes <- sort(unique(c(x$V1,x$V2)))
res <- list()
res$network <- x
## random complex
rnd <- sort(all_genes) %>%
combn(m=2) %>%
t %>%
as.data.frame(stringsAsFactors=FALSE)
yyy <- paste(res$network$V1,res$network$V2,sep="|")
rnd_p <- paste(rnd$V1,rnd$V2,sep="|")
rnd <- rnd[ !(rnd_p %in% yyy), ]
res$rnd_network <- rnd
return(res)
}
calc_network_corr=function(x,y,sample_class=NULL,missing_value_ratio=0.00001,cpu=0){
res <- list()
res$use_sample_class <- sample_class
if(!is.null(sample_class)){
cat("Only use samples from class: ",
paste(sample_class,collapse = ","),"\n")
x <- lapply(x, function(xx){
rm_samples <- xx@sampleList %>% filter(!(class %in% sample_class))
yy <- removeSample(xx,rsamples = rm_samples$sample)
cat("Use samples:\n")
out <- yy@peaksData %>% dplyr::select(sample,class) %>%
dplyr::distinct() %>%
dplyr::group_by(class) %>%
dplyr::summarise(n=n())
print(out)
return(yy)
})
}else{
cat("Use all samples!\n")
}
res$missing_value_ratio <- missing_value_ratio
if(missing_value_ratio>0){
cat("Filter missing value with cutoff: ",missing_value_ratio,"\n")
x <- lapply(x, function(xx){
cat("dataset:",xx@ID,"\n")
yy <- xx %>% filterPeaks(ratio=missing_value_ratio)
return(yy)
})
}else{
cat("No missing value filtering!\n")
}
for(i in 1:length(x)){
if(i == 1){
common_ids <- x[[1]]@peaksData$ID
}else{
common_ids <- intersect(common_ids,x[[i]]@peaksData$ID)
}
}
common_ids <- unique(common_ids)
res$common_ids <- common_ids
res$n_used_pairs <- nrow(y$network)
cat("Common genes:",length(common_ids),"\n")
cat("Used complexes:",nrow(y$network),"\n")
if(nrow(y$rnd_network) > 10*nrow(y$network)){
set.seed(1900)
n_rnd <- nrow(y$rnd_network)
n <- nrow(y$network)
ind <- sample(1:n_rnd,2*n)
y$rnd_network <- y$rnd_network[ind,]
cat("Total random complexes:",n_rnd,"\n")
cat("Used random complexes:",nrow(y$rnd_network),"\n")
}
if(cpu==0){
cpu <- detectCores()
}
if(cpu > length(x)){
cpu <- length(x)
}
cat("Used cpus:",cpu,"\n")
cl <- makeCluster(getOption("cl.cores", cpu))
clusterExport(cl, c("calc_pair_cor"),envir=environment())
clusterExport(cl, c("calcCor"),envir=environment())
#clusterEvalQ(cl,library("dplyr"))
tmp_res <- parLapply(cl,x,fun = calc_pair_cor,pdat=y,idset=common_ids)
res$result <- tmp_res
stopCluster(cl)
return(res)
}
calc_pair_cor=function(para,pdat,idset=NULL){
## filter genes
valid_ids <- para@peaksData$ID %>% unique
pdat$network <- pdat$network %>%
filter(V1 %in% valid_ids, V2 %in% valid_ids)
pdat$rnd_network <- pdat$rnd_network %>%
filter(V1 %in% valid_ids, V2 %in% valid_ids)
## only consider the genes provided by users
if(!is.null(idset)){
valid_ids <- idset
pdat$network <- pdat$network %>%
filter(V1 %in% valid_ids, V2 %in% valid_ids)
pdat$rnd_network <- pdat$rnd_network %>%
filter(V1 %in% valid_ids, V2 %in% valid_ids)
}
pdat$network <- pdat$network %>% mutate(nid=1:nrow(pdat$network))
dat <- pdat$network %>% gather(key="network",value="ID",-nid)
m <- inner_join(dat,para@peaksData,by="ID")
res <- list()
res$network <- m %>% group_by(nid) %>%
do(cor=calcCor(.)) %>%
ungroup() %>%
as.data.frame() %>%
mutate(cor=unlist(cor))
## random complex
pdat$rnd_network <- pdat$rnd_network %>% mutate(nid=1:nrow(pdat$rnd_network))
dat <- pdat$rnd_network %>% gather(key="network",value="ID",-nid)
m <- inner_join(dat,para@peaksData,by="ID")
res$rnd_network <- m %>% group_by(nid) %>%
do(cor=calcCor(.)) %>%
ungroup() %>%
as.data.frame() %>%
mutate(cor=unlist(cor))
return(res)
}
calcCor=function(x){
x <- x %>% dplyr::select(ID,sample,value,nid) %>%
spread(key=ID,value=value) %>%
as.data.frame()
res <- cor(x[,3],x[,4],use = "com",method = "spearman")
return(res)
}
plot_network_cor=function(x, out_dir="./",prefix="test"){
res <- list()
dataset_name <- names(x)
for(i in 1:length(dataset_name)){
x[[i]]$network$nid <- dataset_name[i]
x[[i]]$rnd_network$nid <- dataset_name[i]
}
ks_res <- x %>% lapply(function(x){
x1 <- unlist(x$network$cor)
x2 <- unlist(x$rnd_network$cor)
# rr <- ks.test(x1,x2)
# rr$statistic
roc_input <- dplyr::bind_rows(data.frame(class="network",cor=x1,stringsAsFactors = FALSE),
data.frame(class="rnd_network",cor=x2,stringsAsFactors = FALSE))
roc.obj <- pROC::roc(roc_input$class,roc_input$cor,percent = FALSE)
auroc_value <- roc.obj$auc
return(auroc_value)
})
res$ks <- ks_res
#ks_file <- paste(out_dir,"/",prefix,"-ks.txt",sep="")
#write(ks_res,file = ks_file)
#ppi_cor_res$res11$complex$nid <- "RNA"
network_data_list <- lapply(x, function(y){
return(y$network)
})
network_data_df <- bind_rows(network_data_list) %>% mutate(type="IntraComplex")
rnd_network_data_list <- lapply(x, function(y){
return(y$rnd_network)
})
rnd_network_data_df <- bind_rows(rnd_network_data_list) %>% mutate(type="InterComplex")
dat <- rbind(network_data_df,rnd_network_data_df)
save(dat,file = "dat.rda")
fig <- paste(out_dir,"/",prefix,"-complex_all_boxplot.png",sep="")
res$network_boxplot <- fig
png(fig,width = 1000,height = 500,res=160)
aaaa <- dat %>% group_by(nid) %>%
filter(type=="InterComplex") %>%
summarize(c=median(cor))
#arrange(c)
dat -> xxxx
dat$nid <- factor(dat$nid,levels = aaaa$nid )
gg <- ggplot(dat,aes(x=nid,y=cor,colour=type))+
theme_bw()+
#geom_violin(position="dodge")+
geom_boxplot(position="dodge")+
xlab("data table")+
ylab("Correlation")+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
print(gg)
dev.off()
save(dat,file="complex_res.rda")
cor_res <- dat %>% group_by(nid,type) %>%
summarise(cor=median(cor)) %>%
ungroup() %>%
spread(key="type",value="cor") %>%
mutate(diff=IntraComplex-InterComplex) %>%
rename(dataSet=nid)
ks_df <- data.frame(dataSet=names(ks_res),ks=ks_res %>% unlist())
res$cor <- merge(cor_res,ks_df,sort=FALSE)
## save plot data
fig_data_file <- paste(out_dir,"/",prefix,"-complex_all_boxplot.rds",sep="")
fig_data <- list()
fig_data$plot_object <- gg
fig_data$plot_data <- dat
saveRDS(fig_data,file = fig_data_file)
return(res)
}
## protein rna correlation
# "pearson", "kendall", "spearman"
# This function is used to calculate the correlation between protein and mRNA.
calcCorBetweenProteinAndRNA=function(para1,para2,log=TRUE,select_by=1,top_n=1000,
geneset = NULL,
outdir = "./",use_class=NULL,
reg=FALSE,
cor_method="spearman",valueID="value"){
prefix <- para1@ID
if(!is.null(use_class)){
cat("Only use class ",use_class,"\n")
if(reg == FALSE){
para1@peaksData <- para1@peaksData %>% dplyr::filter(class==use_class)
para2@peaksData <- para2@peaksData %>% dplyr::filter(class==use_class)
}else{
para1@peaksData <- para1@peaksData %>% dplyr::filter(grepl(pattern = use_class,x = class))
para2@peaksData <- para2@peaksData %>% dplyr::filter(grepl(pattern = use_class,x = class))
}
cat("Used samples:",length(unique(para1@peaksData$sample)),length(unique(para2@peaksData$sample)),"\n")
}else{
cat("Use all samples!\n")
}
if(!is.null(geneset)){
cat("Only use specified gene set\n")
para1@peaksData <- para1@peaksData %>% dplyr::filter(ID %in% geneset)
para2@peaksData <- para2@peaksData %>% dplyr::filter(ID %in% geneset)
}
x1 <- para1@peaksData %>% dplyr::select(ID,sample,!!valueID)
x2 <- para2@peaksData %>% dplyr::select(ID,sample,!!valueID)
m <- merge(x1,x2,by=c("ID","sample"))
## firstly, plot a scatterplot: x axis is sd, y axis is cor
names(m)[3:4] <- c("x","y")
m$x[m$x <= 0] <- NA
m$y[m$y <= 0] <- NA
if(log==TRUE){
m$x <- log2(m$x)
m$y <- log2(m$y)
}
## filter genes with all values are missing
na_row <- apply(m[,3:4],1,function(x){any(is.na(x))})
m <- m[!na_row,]
## filter genes occurring less than 3 samples
n3 <- m %>% group_by(ID) %>% summarise(n=sum(!is.na(x) & !is.na(y))) %>% filter(n>=3)
m <- m %>% filter(ID %in% n3$ID)
## Sample wise correlation
save(m,file = "m.rda")
sample_wise_cor <- m %>% group_by(sample) %>%
dplyr::summarise(cor=cor(x,y,use = "com",method = cor_method)) %>%
dplyr::mutate(label=prefix)
## Gene wise correlation
res <- m %>% group_by(ID) %>% dplyr::summarise(sd_x=sd(x,na.rm = TRUE),
sd_y=sd(y,na.rm=TRUE),
cor=cor(x,y,use = "com",method = cor_method),
pvalue=get_cor_pvalue(x,y,cor_method = cor_method))
#save(para1,para2,valueID,na_row,m,res,file="test111.rda")
png(paste(outdir,"/",prefix,"-cor-sd.png",sep=""),width = 600,height = 1200,res=120)
cor_name <- paste("Cor (",cor_method,")",sep="")
par(mfrow=c(4,2),mar=c(3,3,2,1),mgp=c(1.6,0.6,0))
plot(res$sd_x,res$cor, xlab="SD (1)",ylab=cor_name,col=rgb(0.5,0.4,0.7,0.4),pch=15,cex=0.6)
plot(res$sd_y,res$cor, xlab="SD (2)",ylab=cor_name,col=rgb(0.5,0.4,0.7,0.4),pch=15,cex=0.6)
plot(res$sd_x,res$sd_y, xlab="SD (1)",ylab="SD (2)",col=rgb(0.5,0.4,0.7,0.4),pch=15,cex=0.6)
hist(res$sd_x,nclass=50,xlab="SD (1)",main="")
hist(res$sd_y,nclass=50,xlab="SD (2)",main="")
cor_mean <- sprintf("%.4f",mean(res$cor,na.rm=TRUE))
cor_md <- sprintf("%.4f",median(res$cor,na.rm=TRUE))
hist(res$cor,nclass=50,xlab=cor_name,main=paste("Data point: ",nrow(res),"\nmean = ",cor_mean,", median = ",cor_md,sep = ""))
dev.off()
rr <- data.frame(n=nrow(res),
n5=sum(res$cor>=0.5,na.rm = TRUE),
n6=sum(res$cor>=0.6,na.rm = TRUE),
n7=sum(res$cor>=0.7,na.rm = TRUE),
n8=sum(res$cor>=0.8,na.rm = TRUE),
median_cor=median(res$cor,na.rm=TRUE),
stringsAsFactors = FALSE)
fres <- list()
fres$feature_wise_cor_data <- res
fres$feature_wise <- rr
fres$sample_wise <- sample_wise_cor
fres$hist_fig <- plot_mRNA_protein_hist(fres$feature_wise_cor_data,out_dir = outdir,prefix = prefix)
return(fres)
}
get_cor_pvalue=function(x,y,cor_method="spearman"){
cor_res <- cor.test(x,y,method = cor_method,use = "com")
return(cor_res$p.value)
}
plot_mRNA_protein_hist=function(x,out_dir="./",prefix="test"){
br <- c(seq(min(x$cor,na.rm = TRUE)-0.1,-0.000000001,by=0.01),
seq(0,max(x$cor,na.rm = TRUE),by=0.01))
max_count <- max(table(cut(x$cor,breaks = br)))
positive_cor_ratio <- sum(x$cor>0,na.rm = TRUE)/nrow(x)
x$pvalue <- p.adjust(x$pvalue,method = "BH")
sig_positive_cor_ratio <- sum(x$pvalue <= 0.01,na.rm = TRUE)/nrow(x)
mean_cor <- mean(x$cor,na.rm = TRUE)
n_pairs <- nrow(x)
positive_cor_ratio <- sprintf("%.2f%%",100*positive_cor_ratio)
sig_positive_cor_ratio <- sprintf("%.2f%%",100*sig_positive_cor_ratio) #%>% paste("*\\'%\\'~",sep = "")
gg <- x %>% mutate(col=ifelse(cor>0,"red","blue")) %>%
ggplot(aes(x=cor,fill=col)) +
geom_histogram(breaks=br,color="black",size=0.1)+
xlab("Spearman’s correlation")+
ylab("Frequency")+
ggpubr::theme_pubr()+
theme(legend.position = "none")+
scale_fill_manual(breaks = c("red", "blue"),
values=c("red", "blue"))+
geom_vline(xintercept = mean_cor,linetype=2)+
annotate("text",x=mean_cor,y=1.05*max_count,label=paste("Mean = ",sprintf("%.2f",mean_cor),sep = ""),
hjust=-0.05,size=3.5)+
annotate("text",x=min(br)+0.01,y=1.05*max_count,
label=paste0(n_pairs," pairs\n",
positive_cor_ratio," positive correlation\n",
sig_positive_cor_ratio," significant\npositive correlation\n(adjusted P <= 0.01)"),
vjust=1,hjust=0,size=3.5)
fig <- paste(out_dir,"/",prefix,"-mRNA-protein-hist.png",sep = "")
png(fig,width = 650,height = 500,res=150)
print(gg)
dev.off()
gg_obj <- paste(out_dir,"/",prefix,"-mRNA-protein-hist_ggplot.rds",sep = "")
saveRDS(gg,file = gg_obj)
return(fig)
}
calc_protein_rna_corr=function(x,rna,sample_class=NULL,out_dir="./",cpu=0,
missing_value_ratio=0.5){
x <- lapply(x, function(xx){
yy <- metaX::filterPeaks(xx,ratio = missing_value_ratio)
return(yy)
})
if(!is.null(sample_class)){
cat("Only use samples from class: ",
paste(sample_class,collapse = ","),"\n")
x <- lapply(x, function(xx){
rm_samples <- xx@sampleList %>% filter(!(class %in% sample_class))
yy <- removeSample(xx,rsamples = rm_samples$sample)
cat("The number of samples used in correlation analysis:",yy@peaksData$sample %>% unique %>% length,"\n")
return(yy)
})
}else{
cat("Use all samples!\n")
}
for(i in 1:length(x)){
if(i == 1){
common_ids <- x[[1]]@peaksData$ID
}else{
common_ids <- intersect(common_ids,x[[i]]@peaksData$ID)
}
}
common_ids <- intersect(common_ids,rna@peaksData$ID)
common_ids <- unique(common_ids)
cat("Common IDs:", length(common_ids),"\n")
if(cpu==0){
cpu <- detectCores()
}
if(cpu > length(x)){
cpu <- length(x)
}
cat("Use cpus:",cpu,"\n")
cl <- makeCluster(getOption("cl.cores", cpu))
res <- parLapply(cl,x,fun = calcCorBetweenProteinAndRNA,
para2=rna,
log=TRUE,
select_by=1,
top_n=1000,
cor_method="spearman",
valueID="value")
stopCluster(cl)
dataset_name <- names(res)
for(i in 1:length(dataset_name)){
res[[i]]$feature_wise$dataSet <- dataset_name[i]
res[[i]]$sample_wise$dataSet <- dataset_name[i]
res[[i]]$feature_wise_cor_data$dataSet <- dataset_name[i]
}
## feature wise
cor_res <- lapply(res, function(y){y$feature_wise}) %>% bind_rows()
cor_res <- bind_cols(cor_res %>% select(dataSet),cor_res %>% select(-dataSet))
## sample wise
sample_wise_cor_res <- lapply(res, function(y){y$sample_wise}) %>% bind_rows()
sample_wise_cor_res_table <- sample_wise_cor_res %>% group_by(dataSet) %>%
summarise(median_cor=median(cor,na.rm = TRUE))
fres <- list(feature_wise_cor_table = cor_res,
sample_wise_cor_table = sample_wise_cor_res_table,
data=res,
use_sample_class = sample_class)
## plot
gg <- ggplot(sample_wise_cor_res,aes(x=dataSet,y=cor,color=dataSet)) +
geom_violin()+
geom_boxplot(width=0.1)+
geom_text(data = sample_wise_cor_res_table, aes(x = dataSet, y = median_cor, label = sprintf("%.4f",median_cor)),
size = 5, vjust = -0.2, color = "black",)+
xlab("data table")+
ylab("Sample wise correlation")+
theme(legend.position="none")
sample_wise_cor_fig <- paste(out_dir,"/sample_wise_cor.png",sep="")
png(sample_wise_cor_fig,width = 800,height = 400,res=120)
print(gg)
dev.off()
fres$sample_wise_cor_fig <- sample_wise_cor_fig
# generate CDF and boxplot plots for gene-wise correlation analysis
## boxplot
gene_wise_cor_res_all <- lapply(res, function(y){y$feature_wise_cor_data}) %>% bind_rows()
fres$gene_wise_cor_boxplot_fig <- plot_feature_wise_cor_boxplot(gene_wise_cor_res_all,out_dir = out_dir,prefix = "OmicsEV")
## CDF
fres$gene_wise_cor_cdf_fig <- plot_feature_wise_cor_cdf(gene_wise_cor_res_all,out_dir = out_dir,prefix = "OmicsEV")
return(fres)
}
plot_feature_wise_cor_cdf=function(x, out_dir="./",prefix="test"){
fig <- paste(out_dir,"/",prefix,"-feature_wise_cor_cdf.png",sep="")
png(fig,width = 550,height = 550,res=150)
gg <- ggplot(x,aes(cor,group=dataSet,color=dataSet))+
stat_ecdf(geom = "step",size=0.3)+
theme(legend.position = c(0, 1),
legend.justification = c(0, 1),
legend.key = element_blank(),
legend.background=element_blank())+
ylab("CDF")+
xlab("Correlation")
print(gg)
dev.off()
## save plot data
fig_data_file <- paste(out_dir,"/",prefix,"-feature_wise_cor_cdf.rds",sep="")
fig_data <- list()
fig_data$plot_object <- gg
fig_data$plot_data <- x
saveRDS(fig_data,file = fig_data_file)
return(fig)
}
plot_feature_wise_cor_boxplot=function(x, out_dir="./",prefix="test"){
median_cor <- x %>% group_by(dataSet) %>% summarise(median_cor=median(cor,na.rm = TRUE)) %>% arrange(desc(median_cor))
x$dataSet <- factor(x$dataSet,levels = median_cor$dataSet)
fig <- paste(out_dir,"/",prefix,"-feature_wise_cor_boxplot.png",sep="")
png(fig,width = 800,height = 400,res=120)
gg <- ggplot(x,aes(x=dataSet,y=cor))+
geom_boxplot(width=0.5,outlier.size=0.2)+
ylab("Correlation")+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
print(gg)
dev.off()
return(fig)
}
calc_basic_metrics=function(x,class_color=NULL,out_dir="./",cpu=0){
if(!is.null(class_color)){
class_col <- read_tsv(class_color) %>% filter(!is.na(col))
}else{
class_col <- NULL
}
dataset_name <- names(x)
xx <- list()
for(i in 1:length(dataset_name)){
xx[[i]] <- list(name=dataset_name[i],data=x[[i]])
}
names(xx) <- dataset_name
## output
fres <- list()
plist <- list()
plist$scale <- "pareto"
plist$t <- 3
plist$missValueRatioQC <- 0.5
plist$missValueRatioSample <- 0.5
plist$classCol <- class_col
plist$pcaLabel <- NA
plist$center <- TRUE
plist$pngWidth <- 6.5
plist$pngHeight <- 6.5
plist$pdfWidth <- 6
plist$pdfHeight <- 6
plist$legendRowBatch <- 4
plist$legendRowClass <- 4
if(length(dataset_name) >=2){
fig <- plot_upset(x,out_dir = out_dir, prefix = "total_identification_upset")
fres$datasets_id_overlap <- fig
}
density_fig <- plot_density(x,out_dir = out_dir, prefix = "sample_wise_density")
if(cpu==0){
cpu <- detectCores()
}
ncpu = cpu
if(cpu > length(dataset_name)){
cpu <- length(dataset_name)
}
cat("Use cpu:",cpu,"\n")
if(cpu > 1){
cl <- makeCluster(getOption("cl.cores", cpu))
clusterExport(cl, c("metaXpipe"),envir=environment())
res <- parLapply(cl,xx,fun = run_basic_metrics,plist=plist,out_dir=out_dir)
stopCluster(cl)
}else{
res <- lapply(xx,run_basic_metrics,plist=plist,out_dir=out_dir)
}
fres$parameters <- plist
fres$density_plot <- density_fig
fres$datasets <- res
## distribution
data_dv <- calc_metrics_for_data_distribution_roc(x,cpu=ncpu)
fres$quant_median_metric <- data_dv
return(fres)
}
min_max_scale = function(x){
return((x- min(x)) /(max(x)-min(x)))
}
plot_upset=function(x,out_dir="./",prefix="test"){
a <- lapply(x, function(y){
return(y@peaksData$ID %>% unique())
})
fig <- paste(out_dir,"./",prefix,".png",sep = "")
png(fig, width = 9.5, height = 5.5,res=120, units = "in")
m <- ComplexHeatmap::make_comb_mat(a)
ss = ComplexHeatmap::set_size(m)
cs = ComplexHeatmap::comb_size(m)
ht <- ComplexHeatmap::UpSet(m,
set_order = order(ss),
comb_order = order(ComplexHeatmap::comb_degree(m), -cs),
top_annotation = HeatmapAnnotation(
"Intersection size" = anno_barplot(cs,
ylim = c(0, max(cs)*1.1),
border = FALSE,
gp = grid::gpar(fill = "black"),
height = grid::unit(4, "cm")
),
annotation_name_side = "left",
annotation_name_rot = 90))
ht <- ComplexHeatmap::draw(ht)
od <- ComplexHeatmap::column_order(ht)
ComplexHeatmap::decorate_annotation("Intersection size", {
grid::grid.text(cs[od], x = seq_along(cs), y = grid::unit(cs[od], "native") + grid::unit(3, "pt"),
default.units = "native", just = "bottom", gp = grid::gpar(fontsize = 9))
})
dev.off()
return(fig)
}
run_basic_metrics=function(x,plist,out_dir="./"){
res <- list()
para <- x$data
dat_name <- x$name
c_out_dir <- paste(out_dir,"/",dat_name,sep="")
para@outdir <- c_out_dir
para@prefix <- dat_name
if(is.null(para@sampleList) || is.na(para@sampleList) ||
nrow(para@sampleList) ==0){
sampleList <- read.delim(para@sampleListFile,stringsAsFactors = FALSE)
}else{
sampleList <- para@sampleList
}
checkSampleList(sampleList)
if(is.null(para@ratioPairs)){
stop("Please set the value of ratioPairs!")
}
##
plsdaPara <- new("plsDAPara")
plsdaPara@scale <- plist$scale
plsdaPara@t <- plist$t
plsdaPara@do <- FALSE
missValueRatioQC <- plist$missValueRatioQC
missValueRatioSample <- plist$missValueRatioSample
classCol <- plist$classCol
pcaLabel <- plist$pcaLabel
center <- plist$center
## directory structure, data (all figures and other files), report.html
raw_outdir <- para@outdir
## all figures and other files
para@outdir <- paste(para@outdir,"/data",sep="")
makeDirectory(para)
if(is.null(para@sampleList) || is.na(para@sampleList) ||
nrow(para@sampleList) ==0){
expDesign <- read.delim(para@sampleListFile,stringsAsFactors = FALSE)
}else{
expDesign <- para@sampleList
}
expDesign$class <- as.character(expDesign$class)
expDesign$class[is.na(expDesign$class)] <- "QC"
expClass <- expDesign %>% group_by(class) %>% summarise(n=length(sample))
expBatch <- expDesign %>% group_by(batch) %>% summarise(n=length(sample))
message(date(),"\treset the peaksData...")
para <- reSetPeaksData(para = para)
##
res$total_features <- para@peaksData$ID %>% unique() %>% length()
## 2. pre-processing, filter noise peaks
if(hasQC(para)){
message(date(),"\tfilter peaks in QC sample:>=",missValueRatioQC)
pre_para <- para
para <- filterQCPeaks(para = para,ratio = missValueRatioQC)
rpeak <- length(unique(pre_para@peaksData$ID)) - length(unique(para@peaksData$ID))
paste("Remove peaks in QC samples with missing
value greater than ",100*missValueRatioQC, " percent:",rpeak,".",sep="")
rm(pre_para)
}
message(date(),"\tfilter peaks in non-QC sample:>=",missValueRatioSample)
pre_para <- para
para <- filterPeaks(para = para,ratio = missValueRatioSample)
rpeak <- length(unique(pre_para@peaksData$ID)) - length(unique(para@peaksData$ID))
##
res$valid_features <- para@peaksData$ID %>% unique() %>% length()
paste("Remove peaks in non-QC samples with missing
value greater than ",100*missValueRatioSample," percent:",rpeak,".")
rm(pre_para)
message("The CV distribution before normalization: ")
printCV(para,valueID = "value")
## 3. Quality control
message("plot peak number distribution...")
fig <- plotPeakNumber(para,legendRowBatch=10)
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
##
res$features_number_distribution <- fig$fig
message("plot CV distribution...")
class_sample <- para@peaksData %>% dplyr::select(class,sample) %>% distinct()
class_2 <- which(table(class_sample$class)>=2) %>% names()
cv_para <- para
cv_para@peaksData <- cv_para@peaksData %>% dplyr::filter(class %in% class_2)
cv_para@sampleList <- cv_para@sampleList %>% dplyr::filter(class %in% class_2)
fig <- metaX::plotCV(cv_para)
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
##
res$cv_distribution <- fig$fig
res$cv_stat <- get_cv(cv_para@peaksData) %>% mutate(dataSet=cv_para@ID)
message("plot missing value distribution...")
fig <- plotMissValue(para,height = 3.7,width = 9.15)
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
##
res$missing_value_distribution <- fig$fig
res$non_missing_value_ratio <- data.frame(dataSet=para@ID,
ratio=get_missing_value_ratio(para@peaksData),
stringsAsFactors = FALSE)
message("plot peak intensity distribution...")
fig <- plotIntDistr(para,width = 9.15)
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
##
res$features_quant_distribution <- fig$fig
message("plot TIC distribution...")
fig <- plotPeakSumDist(para,valueID = "value",height = 3.7,width = 5.5)
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
ticTable <- para@peaksData %>% group_by(sample,batch,class) %>%
dplyr::summarise(tic=sum(value,na.rm = TRUE)) %>%
group_by(batch) %>%
dplyr::summarise(n00=quantile(tic)[1],
n25=quantile(tic)[2],
n50=quantile(tic)[3],
n75=quantile(tic)[4],
n100=quantile(tic)[5])
print(ticTable)
message("plot average intensity distribution...")
fig <- plotPeakMeanDist(para,valueID = "value")
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
meanTable <- para@peaksData %>% group_by(sample,batch,class) %>%
dplyr::summarise(meanIntensity=mean(value,na.rm = TRUE)) %>%
group_by(batch) %>%
dplyr::summarise(n00=quantile(meanIntensity)[1],
n25=quantile(meanIntensity)[2],
n50=quantile(meanIntensity)[3],
n75=quantile(meanIntensity)[4],
n100=quantile(meanIntensity)[5])
print(meanTable)
message(date(),"\tmissing value inputation...")
para <- missingValueImpute(para,cpu=1)
paste("The missing value were imputed by ",para@missValueImputeMethod,".",sep="")
save(para,file = "para.rds")
fig <- plotCorHeatmap(para = para,valueID = "value",anno = TRUE,
samples = NULL,
height = 6,width = 6,sortBy="batch",cluster = FALSE,
classCol=classCol)
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
##
res$sample_wise_cor_heatmap <- fig$fig
if(hasQC(para)){
message("plot correlation heatmap...")
#saveRDS(para,file = "para.rda")
prefix_w <- para@prefix
para@prefix <- paste(prefix_w,"-qc",sep = "")
fig <- plotCorHeatmap(para = para,valueID = "value",anno = TRUE,
samples = NA, ## QC sample
height = 6,width = 6,cluster = FALSE,
classCol=classCol)
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
para@prefix <- prefix_w
res$qc_sample_wise_cor_heatmap <- fig$fig
}
#save(para,t,file="test.rda")
save(para,plist,scale,center,pcaLabel,classCol,file=paste(para@outdir,"/","pca.rda",sep=""))
ppca <- transformation(para,method = plist$t,valueID = "value")
ppca <- metaX::preProcess(ppca,scale = plist$scale,center = plist$center,
valueID = "value")
fig <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = TRUE,
rmQC = FALSE,label = pcaLabel,classColor = classCol,
legendRowClass = plist$legendRowClass,
legendRowBatch = plist$legendRowBatch,
pngWidth = plist$pngWidth,
pngHeight = plist$pngHeight,
pdfWidth = plist$pdfWidth,
pdfHeight = plist$pdfHeight)
prefix_bak <- ppca@prefix
ppca@prefix <- paste(ppca@prefix,"-nobatch",sep="")
fig_nobatch <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = FALSE,
rmQC = FALSE,label = pcaLabel,classColor = classCol,
legendRowClass = plist$legendRowClass,
legendRowBatch = plist$legendRowBatch,
pngWidth = plist$pngWidth,
pngHeight = plist$pngHeight,
pdfWidth = plist$pdfWidth,
pdfHeight = plist$pdfHeight)
##
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
res$pca_with_batch <- fig$fig
res$pca_plot_data <- fig$plotdata
## plot heatmap
save(ppca,classCol,file = "heatmap.rda")
fig <- plotHeatMap(ppca,valueID="value",log=FALSE,rmQC=FALSE,
#scale="row",
#clustering_distance_rows="euclidean",
#clustering_distance_cols="euclidean",
#clustering_method="ward.D2",
seriation = FALSE,
classCol=classCol,
show_colnames=FALSE)
ppca@prefix <- prefix_bak
#fig1 <- paste("data/",basename(fig$fig),sep="")
#fig2 <- paste("data/",basename(fig$highfig),sep="")
res$cluster_heatmap <- fig$fig
## remove QC
ppca@prefix <- paste(ppca@prefix,"-noqc",sep="")
fig <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = TRUE,
rmQC = TRUE,label = pcaLabel,classColor = classCol,
legendRowClass = plist$legendRowClass,
legendRowBatch = plist$legendRowBatch,
pngWidth = plist$pngWidth,
pngHeight = plist$pngHeight,
pdfWidth = plist$pdfWidth,
pdfHeight = plist$pdfHeight)
## no QC, no batch
ppca@prefix <- paste(ppca@prefix,"-nobatch",sep="")
fig <- metaX::plotPCA(ppca,valueID = "value",scale = "none",batch = FALSE,
rmQC = TRUE,label = pcaLabel,classColor = classCol,
legendRowClass = plist$legendRowClass,
legendRowBatch = plist$legendRowBatch,
pngWidth = plist$pngWidth,
pngHeight = plist$pngHeight,
pdfWidth = plist$pdfWidth,
pdfHeight = plist$pdfHeight)
#plotPLSDA(para)
if(!is.null(para@ratioPairs) && para@ratioPairs != ""){
## only perform this analysis when a user provides comparison group
## information.
para@peaksData$valueNorm <- para@peaksData$value
#para <- peakStat(para = para,plsdaPara = plsdaPara,doROC=FALSE,
# pcaLabel=pcaLabel,classColor = pcaColor)
qf <- paste("data/",para@prefix,"-quant.txt",sep="")
}
saveRDS(para,file = paste(para@outdir,"/",para@prefix,"-result.rds",
sep=""))
message("Print information about the current R session:")
writeLines(capture.output(sessionInfo()),
paste(para@outdir,"/",para@prefix,"-sessionInfo.txt",sep=""))
return(res)
}
get_sample_data=function(x){
a <- lapply(x, function(y){
yy <- y@peaksData %>% dplyr::select(sample,class,batch) %>%
dplyr::distinct() %>%
dplyr::mutate(dataSet=y@ID)
yy$sample <- as.character(yy$sample)
yy$class <- as.character(yy$class)
return(yy)
})
xx <- dplyr::bind_rows(a)
class_samples <- xx %>% group_by(class,dataSet) %>%
dplyr::summarise(samples=dplyr::n()) %>%
ungroup() %>%
tidyr::spread(key=dataSet,value=samples)
return(class_samples)
}
get_identification_summary_table=function(x,format=TRUE){
x1 <- get_metrics(x$basic_metrics$datasets,metric="total_features")
x2 <- get_metrics(x$basic_metrics$datasets,metric="valid_features")
a <- data.frame(dataSet=names(x1),x1=x1,x2=x2,stringsAsFactors = FALSE)
if(format==TRUE){
a <- a %>% mutate(x1=cell_spec(x1, bold = ifelse(x1 >= max(x1), TRUE, FALSE), color = ifelse(x1 >= max(x1), "red", "black"))) %>%
mutate(x2=cell_spec(x2, bold = ifelse(x2 >= max(x2), TRUE, FALSE), color = ifelse(x2 >= max(x2), "red", "black")))
# dplyr::mutate(x1 = formattable::color_bar("lightgreen")(x1),
# x2 = formattable::color_bar("lightgreen")(x2)) #%>%
#dplyr::mutate(x11 = ifelse(x1 >= max(x1),
# cell_spec(x1, color = "red", bold = T),
# cell_spec(x1, color = "green", italic = T)),
# x22 = ifelse(x2 >= max(x2),
# cell_spec(x2, color = "red", bold = T),
# cell_spec(x2, color = "green", italic = T))) %>%
#rename(x1=NULL,x2=NULL)
}
names(a) <- c("dataSet","#identified features","#quantifiable features")
rownames(a) <- NULL
return(a)
}
get_full_path=function(x){
for(i in 1:length(x)){
x[i] <- normalizePath(x[i])
}
return(x)
}
run_kbet=function(x,out_dir="./",prefix="test"){
res <- lapply(x, function(y){
xx <- metaX:::getPeaksTable(y,value="value")
png(paste(y@ID,".png",sep=""))
if(min(table(xx$batch)) <= 10){
dat <- kBET(df=xx[,-c(1:4)],batch = xx$batch,plot = TRUE,
heuristic = FALSE,k0 = 3)
}else{
dat <- kBET(df=xx[,-c(1:4)],batch = xx$batch,plot = TRUE)
}
ff <- list(kbet=dat,name=y@ID)
dev.off()
return(ff)
})
res_table <- lapply(res, function(a){
y <- a$kbet$summary[1,] %>% mutate(dataSet=a$name) %>% select(dataSet,everything())
})
res_table <- bind_rows(res_table)
dd <- lapply(res, function(a){
y <- a$kbet
dat <- data.frame(value=c(y$stats$kBET.expected,y$stats$kBET.observed),
Test=c(rep("Expected(random)",length(y$stats$kBET.expected)),
rep("Observed(kBET)",length(y$stats$kBET.observed))),
dataSet=a$name,
stringsAsFactors = FALSE)
})
gdat <- bind_rows(dd)
fig <- paste(out_dir,"/",prefix,"-kBET-boxplot.png",sep="")
png(fig,width = 800,height = 500,res=130)
gg <- ggplot(gdat,aes(x=Test,y=value,color=dataSet))+
geom_boxplot()+
ylab("Rejection rate")
#facet_grid(.~dataSet)
print(gg)
dev.off()
fres <- list(kBET_boxplot=fig,table=res_table,kbet=res)
return(fres)
}
calc_silhouette_width=function(x,transform_method="3",scale_method="pareto",
center=TRUE,missing_value_ratio=0.5){
res <- lapply(x, function(y){
y <- metaX::filterPeaks(y,ratio = missing_value_ratio)
y@missValueImputeMethod <- "knn"
y <- metaX::missingValueImpute(y,valueID = "value",method="knn")
ppca <- metaX::transformation(y,method = transform_method,valueID = "value")
ppca <- metaX::preProcess(ppca,scale = scale_method,center = center,
valueID = "value")
xx <- metaX:::getPeaksTable(y,value="value")
#data: a matrix (rows: samples, columns: features (genes))
#batch: vector or factor with batch label of each cell
pca.data <- prcomp(xx[,-c(1:4)], center=FALSE) #compute PCA representation of the data
batch.silhouette <- batch_sil(pca.data, xx$batch)
return(batch.silhouette)
})
res <- res %>% unlist()
return(res)
}
calc_pca_batch_regression=function(x,transform_method=3,scale_method="pareto",
center=TRUE,missing_value_ratio=0.5,format=FALSE,
top_pc=10){
res <- lapply(x, function(y){
y <- metaX::filterPeaks(y,ratio = missing_value_ratio)
y@missValueImputeMethod <- "knn"
y <- metaX::missingValueImpute(y,valueID = "value",method="knn")
ppca <- metaX::transformation(y,method = transform_method,valueID = "value")
ppca <- metaX::preProcess(ppca,scale = scale_method,center = center,
valueID = "value")
xx <- metaX:::getPeaksTable(ppca,value="value")
#data: a matrix (rows: samples, columns: features (genes))
#batch: vector or factor with batch label of each cell
pca.data <- prcomp(xx[,-c(1:4)], center=FALSE) #compute PCA representation of the data
yy <- pcRegression(pca.data,batch = xx$batch,n_top = top_pc)
save(pca.data,yy,xx,top_pc,file = "pca_batch_reg.rda")
return(list(pcr=yy,name=y@ID))
return(yy)
})
save(res,file="res.rda")
fres <- lapply(res, function(a){
r2 <- as.data.frame(a$pcr$r2)
r2 <- r2[1:min(top_pc,nrow(r2)),]
#r2$R.squared <- cell_spec(r2$R.squared, bold = T,
# color = ifelse(r2$p.value.lm <= 0.05, "red", "black"),
# font_size = spec_font_size(r2$R.squared))
#r2$R.squared <- cell_spec(r2$R.squared, )
r2$PC <- row.names(r2)
r2$PC <- as.integer(str_replace_all(r2$PC,pattern = "PC",replacement = ""))
r2$dataSet <- a$name
r2$ExplainedVar <- a$pcr$ExplainedVar[1:nrow(r2)]
return(r2)
})
fres <- bind_rows(fres)
fres$R.squared <- sprintf("%.3f",fres$R.square) %>% as.numeric()
fres$R.squared <- cell_spec(fres$R.squared, bold = T,
color = ifelse(fres$p.value.lm <= 0.05, "red", "black"),
font_size = spec_font_size(fres$R.squared))
fres2 <- fres %>% select(PC,R.squared,dataSet) %>% spread(key=dataSet,value=R.squared)
## explained_var
explained_var <- fres %>% select(PC,ExplainedVar,dataSet) %>% spread(key=dataSet,value=ExplainedVar)
## pcRegscale
pcRegscale <- sapply(res, function(a){
return(a$pcr$pcRegscale)
})
pcRegscale <- data.frame(dataSet=names(pcRegscale),
pcRegscale=pcRegscale,
stringsAsFactors = FALSE)
row.names(pcRegscale) <- NULL
return(list(pcr=res,
table=fres2,
explained_var=explained_var,
pcRegscale=pcRegscale,
top_pc=top_pc))
}
get_pcr_table=function(x,top_pc=10){
fres <- lapply(x, function(a){
r2 <- as.data.frame(a$pcr$r2)
r2 <- r2[1:min(top_pc,nrow(r2)),]
#r2$R.squared <- cell_spec(r2$R.squared, bold = T,
# color = ifelse(r2$p.value.lm <= 0.05, "red", "black"),
# font_size = spec_font_size(r2$R.squared))
#r2$R.squared <- cell_spec(r2$R.squared, )
r2$PC <- row.names(r2)
r2$dataSet <- a$name
return(r2)
})
fres <- bind_rows(fres)
return(fres)
}
calc_batch_effect_metrics=function(x,out_dir="./",prefix="test",missing_value_cutoff=0.5){
cat("Missing value imputation ...\n")
x <- lapply(x, function(y){
y <- metaX::filterPeaks(y,ratio = missing_value_cutoff)
y <- missingValueImpute(y)
return(y)
})
#kbet_res <- run_kbet(x)
sil_res <- calc_silhouette_width(x)
pcr_res <- calc_pca_batch_regression(x)
#res <- list(kbet=kbet_res,sil=sil_res,pcr=pcr_res)
res <- list(kbet=NULL,sil=sil_res,pcr=pcr_res)
return(res)
}
plot_pca=function(x,out_dir="./",prefix="test",pointSize=0.8,labelSize=4,
legendRowBatch = 10,
legendRowClass = NULL,show_class=FALSE,
pc="12",ncol_sub_fig=3,fig_res=120, pngWidth=6.5, legend_pos=NULL){
# pc = "12" or "13"
plotData <- lapply(x,function(y){y$pca_plot_data}) %>% bind_rows()
if(pc == "12"){
ggobj <- ggplot(data = plotData,aes(x=x,y=y,color=batch))
ggobj <- ggobj + geom_hline(yintercept=0,colour="gray")+
geom_vline(xintercept=0,colour="gray")+
#geom_point()+
#xlab(paste("PC1"," (",sprintf("%.2f%%",100*pca.res@R2[1]),") ",sep=""))+
#ylab(paste("PC2"," (",sprintf("%.2f%%",100*pca.res@R2[2]),") ",sep=""))+
xlab("PC1")+
ylab("PC2")+
theme_bw()+
theme(#legend.position = legendPosition,
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
}else{
ggobj <- ggplot(data = plotData,aes(x=x,y=z,color=batch))
ggobj <- ggobj + geom_hline(yintercept=0,colour="gray")+
geom_vline(xintercept=0,colour="gray")+
#geom_point()+
#xlab(paste("PC1"," (",sprintf("%.2f%%",100*pca.res@R2[1]),") ",sep=""))+
#ylab(paste("PC2"," (",sprintf("%.2f%%",100*pca.res@R2[2]),") ",sep=""))+
xlab("PC1")+
ylab("PC3")+
theme_bw()+
theme(#legend.position = legendPosition,
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
}
ggobj <- ggobj + stat_ellipse(aes(),geom = "path")
if(show_class==FALSE){
ggobj <- ggobj + geom_point(size=pointSize)
}else{
ggobj <- ggobj + geom_point(aes(shape=class),size=pointSize)+
scale_shape_manual(values=1:n_distinct(plotData$class))
}
ggobj <- ggobj + facet_wrap(.~dataSet,ncol = ncol_sub_fig,scales="free")
if(!is.null(legendRowBatch)){
ggobj <- ggobj + guides(color = guide_legend(nrow = legendRowBatch))
}
if(!is.null(legend_pos)){
ggobj <- ggobj + theme(legend.position=legend_pos)
}
#pdf(file = highfig,width = pdfWidth,height = pdfHeight)
#pngWidth <- 6.5
ht <- ceiling(length(unique(plotData$dataSet))/ncol_sub_fig) * pngWidth/ncol_sub_fig
fig <- paste(out_dir,"/",prefix,"-pca_batch.png",sep="")
png(fig,width = pngWidth,height = ht,res = fig_res,units = "in")
print(ggobj)
dev.off()
return(list(fig=fig,data=plotData))
}
plot_density=function(x,out_dir="./",prefix="test",filter_by_quantile=0){
a <- lapply(x, function(y){
dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
as_tibble()
dat$value[dat$value<=0] <- NA
return(dat)
})
if(filter_by_quantile>0){
a <- lapply(a, function(y){
b <- quantile(y$value,probs = c(filter_by_quantile,1-filter_by_quantile),na.rm = TRUE)
y <- filter(y,value >= b[1], value <= b[2])
return(y)
})
}
dat <- bind_rows(a)
#dat$value[dat$value<=0] <- NA
gg <- dat %>% ggplot(aes(x=log2(value),color=sample)) +
#stat_density(geom="line")+
geom_line(stat="density",alpha=0.4)+
#geom_density()+
facet_wrap(.~dataSet,ncol = 3,scales="free")+
guides(color=FALSE)+
theme_bw()+
theme(#legend.position = legendPosition,
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
pngWidth <- 6.5
ht <- ceiling(length(unique(dat$dataSet))/3) * pngWidth/3
fig <- paste(out_dir,"/",prefix,"-density.png",sep="")
png(fig,width = pngWidth,height = ht,res = 120,units = "in")
print(gg)
dev.off()
return(fig)
}
calc_ml_metrics=function(x,sample_list=NULL,sample_class=NULL,use_all=TRUE,
missing_value_cutoff=0.5, valueID="value",cpu=0,
n_repeats=20){
if(missing_value_cutoff > 0){
x <- lapply(x, filterPeaks, ratio=missing_value_cutoff)
}
if(use_all==FALSE){
## use common genes/proteins
for(i in 1:length(x)){
if(i == 1){
common_ids <- x[[1]]@peaksData$ID
}else{
common_ids <- intersect(common_ids,x[[i]]@peaksData$ID)
}
}
common_ids <- unique(common_ids)
cat("Common IDs:", length(common_ids),"\n")
x <- lapply(x, function(y){
y@peaksData <- y@peaksData %>% filter(ID %in% common_ids)
return(y)
})
}
if(!is.null(sample_list)){
sam <- read.delim(sample_list,stringsAsFactors = FALSE) %>%
dplyr::select(sample,class)
x <- lapply(x, function(y){
y@sampleList <- read.delim(sample_list,stringsAsFactors = FALSE)
y@peaksData$class<- NULL
y@peaksData <- merge(y@peaksData,sam,by="sample")
return(y)
})
class_group <- sam$class %>% unique()
}else{
if(!is.null(sample_class)){
class_group <- sample_class
}else{
class_group <- x[[1]]@peaksData$class %>% unique()
}
}
cat("Use class:",paste(class_group,collapse = ","),"\n")
res <- lapply(x, function(y){
if(requireNamespace("doMC",quietly = TRUE)){
if(cpu==0){
library(doMC)
cpu = detectCores()
registerDoMC(cores = cpu)
}else if(cpu>=2){
library(doMC)
registerDoMC(cores = cpu)
}
}
if(any(is.na(y@peaksData$value))){
y <- metaX::missingValueImpute(y)
}else{
cat("No missing value!\n")
}
dd <- lapply(1:n_repeats,function(i){
cat("Repeat:", i,"\n")
dat <- featureSelection(y,group=class_group,method = "rf",
valueID = valueID, fold = 5,
resampling_method = "LOOCV",
repeats = 10,
plot_roc = FALSE,
ratio = 2/3, k = 100,
metric = "ROC",
sizes = length(unique(y@peaksData$ID)),
plotCICurve = FALSE,verbose=TRUE)
dat$results$Repeat_ID = i
return(dat)
})
eval_res <- lapply(dd,function(x){ x$results %>% select(everything())}) %>%
bind_rows()
return(eval_res)
#return(dat)
})
name_datasets <- names(res)
for(i in 1:length(res)){
res[[i]]$dataSet <- name_datasets[i]
}
res_table <- bind_rows(res)
ftable <- res_table %>% group_by(dataSet) %>% dplyr::summarise(mean_ROC=mean(ROC),median_ROC=median(ROC),sd_ROC=sd(ROC))
fres <- list(data=res_table,table=ftable,class_group=class_group)
return(fres)
}
plot_ml_boxplot=function(x, out_dir="./",prefix="test"){
mean_roc <- x %>% group_by(dataSet) %>% summarise(mean_ROC=mean(ROC)) %>% arrange(desc(mean_ROC))
x$dataSet <- factor(x$dataSet,levels = mean_roc$dataSet)
fig <- paste(out_dir,"/",prefix,"-ml_boxplot.png",sep="")
png(fig,width = 800,height = 400,res=120)
gg <- ggplot(x,aes(x=dataSet,y=ROC))+
geom_boxplot(width=0.5,outlier.size=0.2)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
stat_summary(fun.y=mean, colour="darkred", geom="point",
shape=18, size=3, show.legend = FALSE) +
geom_text(data = mean_roc %>% mutate(mean_roc_label=sprintf("%.2f",mean_ROC)),
aes(label = mean_roc_label, y = mean_ROC + 0.1*(max(x$ROC)-min(x$ROC))),colour="darkred")
print(gg)
dev.off()
## save plot data
fig_data_file <- paste(out_dir,"/",prefix,"-ml_boxplot.rds",sep="")
fig_data <- list()
fig_data$plot_object <- gg
fig_data$plot_data <- x
saveRDS(fig_data,file = fig_data_file)
return(fig)
}
calc_metrics_for_data_distribution=function(x,cpu=0){
# x is list of datasets
if(cpu==0){
cpu = detectCores()
}
a <- lapply(x, function(y){
# cat("Dataset:",y@ID,"\n")
dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
as_tibble()
dat$value[dat$value<=0] <- NA
# log2 scale
dat$value <- log2(dat$value)
all_samples <- dat$sample %>% unique()
ss <- combn(all_samples,m=2)
cl <- makeCluster(getOption("cl.cores", cpu))
clusterExport(cl, c("get_ks_statistic"),envir=environment())
ksv <- parLapply(cl,1:ncol(ss),fun = get_ks_statistic,ss=ss,dat=dat)
stopCluster(cl)
ksv <- unlist(ksv)
ks <- median(ksv,na.rm = TRUE)
res <- data.frame(dataSet=y@ID,quant_median_ks=ks,n=length(ksv),stringsAsFactors = FALSE) %>%
as_tibble()
return(res)
}) %>% bind_rows()
a$raw_quant_median_ks <- a$quant_median_ks
#a$quant_median_ks <- 1 - min_max_scale(a$quant_median_ks)
a$quant_median_ks <- 1 - a$quant_median_ks/(a$quant_median_ks+1)
return(a)
}
calc_metrics_for_data_distribution_roc=function(x,cpu=0){
# x is list of datasets
if(cpu==0){
cpu = detectCores()
}
a <- lapply(x, function(y){
# cat("Dataset:",y@ID,"\n")
dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
as_tibble()
dat$value[dat$value<=0] <- NA
# log2 scale
dat$value <- log2(dat$value)
all_samples <- dat$sample %>% unique()
ss <- combn(all_samples,m=2)
cl <- makeCluster(getOption("cl.cores", cpu))
clusterExport(cl, c("get_two_sample_auroc"),envir=environment())
ksv <- parLapply(cl,1:ncol(ss),fun = get_two_sample_auroc,ss=ss,dat=dat)
stopCluster(cl)
ksv <- unlist(ksv)
ks <- median(ksv,na.rm = TRUE)
res <- data.frame(dataSet=y@ID,quant_median_metric=ks,n=length(ksv),stringsAsFactors = FALSE) %>%
as_tibble()
return(res)
}) %>% bind_rows()
#a$raw_quant_median_ks <- a$quant_median_ks
#a$quant_median_ks <- 1 - min_max_scale(a$quant_median_ks)
#a$quant_median_ks <- 1 - a$quant_median_ks/(a$quant_median_ks+1)
return(a)
}
get_ks_statistic = function(s,ss,dat){
ii <- as.character(ss[,s])
x1 <- dat$value[dat$sample==ii[1]]
x2 <- dat$value[dat$sample==ii[2]]
b <- ks.test(x1,x2)
return(abs(b$statistic))
}
get_two_sample_auroc = function(s,ss,dat){
ii <- as.character(ss[,s])
x <- dat %>% dplyr::filter(sample %in% ii[1:2]) %>%
dplyr::filter(!is.na(value))
x$sample <- factor(x$sample,levels = ii[1:2])
roc.obj <- pROC::roc(x$sample,x$value,percent = FALSE)
auroc_value <- roc.obj$auc
metric <- 1-2*abs(auroc_value-0.5)
return(metric)
}
# generate overview table and overview radar figure
generate_overview_table=function(x,highlight_top_n=3,min_auc=0.8){
## generate an overview table which contains different metrics
## input: x =>
## two data tables, one for radar plot real value, one for overview table
total_features <- x$input_parameters$total_features
dat <- get_identification_summary_table(x,format = FALSE)
show_dat <- dat
if("#identified features" %in% names(dat) && "#quantifiable features" %in% names(dat)){
dat <- dat %>% dplyr::mutate(`#identified features`=`#identified features`/total_features,
`#quantifiable features`=`#quantifiable features`/total_features)
show_dat <- show_dat %>% dplyr::mutate(`#identified features` = paste(`#identified features`,"\n(",format_number(dat$`#identified features`),")",sep=""),
`#quantifiable features`=paste(`#quantifiable features`,"\n(",format_number(dat$`#quantifiable features`),")",sep=""))
}
## add missing value metric
missing_value_table <- lapply(x$basic_metrics$datasets,function(y)y[["non_missing_value_ratio"]]) %>%
dplyr::bind_rows() %>%
dplyr::rename(non_missing_value_ratio=ratio)
dat <- merge(dat,missing_value_table %>% dplyr::select(dataSet,non_missing_value_ratio),by="dataSet")
show_dat <- merge(show_dat,missing_value_table %>%
dplyr::select(dataSet,non_missing_value_ratio),by="dataSet") %>%
dplyr::mutate(non_missing_value_ratio=format_number(non_missing_value_ratio))
## data distribution metric
if(!("quant_median_metric" %in% names(x$basic_metrics)) && length(x$input_parameters$datasets) >= 2){
ncpu <- ifelse("cpu" %in% names(x$input_parameters),x$input_parameters$cpu,0)
data_dv <- calc_metrics_for_data_distribution_roc(x$input_parameters$datasets,cpu=ncpu)
dat <- merge(dat,data_dv %>%
dplyr::select(dataSet,quant_median_metric) %>%
dplyr::rename(data_dist_similarity=quant_median_metric),by="dataSet")
show_dat <- merge(show_dat,data_dv %>%
dplyr::select(dataSet,quant_median_metric) %>%
dplyr::mutate(data_dist_similarity=format_number(quant_median_metric)) %>%
dplyr::select(-quant_median_metric),by="dataSet")
}else if("quant_median_metric" %in% names(x$basic_metrics)){
data_dv <- x$basic_metrics$quant_median_metric
dat <- merge(dat,data_dv %>%
dplyr::select(dataSet,quant_median_metric) %>%
dplyr::rename(data_dist_similarity=quant_median_metric),by="dataSet")
show_dat <- merge(show_dat,data_dv %>%
dplyr::select(dataSet,quant_median_metric) %>%
dplyr::mutate(data_dist_similarity=format_number(quant_median_metric)) %>%
dplyr::select(-quant_median_metric),by="dataSet")
}
if(!is.null(x$batch_effect_metrics)){
# kBET
if(!is.null(x$batch_effect_metrics$kbet)){
dat <- merge(dat,x$batch_effect_metrics$kbet$table %>%
select(dataSet,kBET.observed))
dat$kBET <- 1 - dat$kBET.observed
dat$kBET.observed <- NULL
show_dat <- merge(show_dat,x$batch_effect_metrics$kbet$table %>%
select(dataSet,kBET.observed))
show_dat$kBET <- paste(format_number(show_dat$kBET.observed),"\n(",format_number(1-show_dat$kBET.observed),")",sep = "")
show_dat$kBET.observed <- NULL
}
#dat$kBET.observed <- cell_spec(dat$kBET.observed,
# color = ifelse(y >= max(y), "red", "black"))
# Silhouette width
sil <- data.frame(dataSet=names(x$batch_effect_metrics$sil),
silhouette_width=x$batch_effect_metrics$sil)
dat <- merge(dat,sil)
dat$silhouette_width <- 1 - abs(dat$silhouette_width)
show_dat <- merge(show_dat,sil)
show_dat$silhouette_width <- paste(format_number(show_dat$silhouette_width),"\n(",format_number(1 - abs(show_dat$silhouette_width)),")",sep = "")
# PCR
# pcr <- get_pcr_table(x$batch_effect_metrics$pcr$pcr)
# sig_pc <- pcr %>% dplyr::filter(p.value.lm<=0.05) %>%
# dplyr::select(PC) %>%
# dplyr::distinct()
# pcr <- pcr %>% dplyr::filter(PC %in% sig_pc$PC) %>%
# dplyr::select(dataSet,PC,R.squared) %>%
# dplyr::mutate(PC=paste("batch_effect_",PC,sep="")) %>%
# tidyr::spread(key = PC,value = R.squared)
pcRegscale <- x$batch_effect_metrics$pcr$pcRegscale
#pcRegscale$pcRegscale <- 1 - pcRegscale$pcRegscale
dat <- merge(dat,pcRegscale)
dat$pcRegscale <- 1 - dat$pcRegscale
show_dat <- merge(show_dat,pcRegscale)
show_dat$pcRegscale <- paste(format_number(show_dat$pcRegscale),"\n(",format_number(1 - show_dat$pcRegscale),")",sep = "")
}
## complex
if(!is.null(x$network_table)){
dat <- merge(dat,x$network_table$cor %>% dplyr::select(dataSet,ks) %>%
dplyr::rename(complex_auc=ks))
show_dat <- merge(show_dat,x$network_table$cor %>% dplyr::select(dataSet,ks) %>%
dplyr::rename(complex_auc=ks)) %>%
dplyr::mutate(complex_auc=format_number(complex_auc))
}
## function prediction
if(!is.null(x$fun_pred)){
f_res <- get_func_pred_meanAUC(x$fun_pred$data,min_auc = min_auc)
dat <- merge(dat,f_res)
show_dat <- merge(show_dat,f_res) %>%
dplyr::mutate(func_auc=format_number(func_auc))
}
## class prediction
if(!is.null(x$ml)){
dat <- merge(dat,x$ml$table %>%
dplyr::select(dataSet,mean_ROC) %>%
dplyr::rename(class_auc=mean_ROC))
show_dat <- merge(show_dat,x$ml$table %>%
dplyr::select(dataSet,mean_ROC) %>%
dplyr::rename(class_auc=mean_ROC)) %>%
dplyr::mutate(class_auc=format_number(class_auc))
}
## median CV: QC samples
cv_table <- get_cv_table(x$basic_metrics$datasets,"cv_stat")
if("QC" %in% cv_table$class){
cv_table <- cv_table %>%
dplyr::filter(class=="QC") %>%
dplyr::select(dataSet,cv30) %>%
dplyr::rename(CV30=cv30) # %>%
#dplyr::mutate(scaled_median_CV=1-median_CV/max(median_CV))
dat <- merge(dat,cv_table,by="dataSet") # %>%
#dplyr::select(-median_CV) %>%
#dplyr::rename(median_CV=scaled_median_CV)
show_dat <- merge(show_dat,cv_table,by="dataSet") %>%
dplyr::mutate(CV30=format_number(CV30)) #%>%
#dplyr::select(-scaled_median_CV)
}
## mRNA-protein correlation
if(!is.null(x$input_parameters$x2)){
dat <- merge(dat,x$protein_rna$feature_wise_cor_table %>%
dplyr::select(dataSet,median_cor) %>%
dplyr::rename(gene_wise_cor=median_cor))
show_dat <- merge(show_dat,x$protein_rna$feature_wise_cor_table %>%
dplyr::select(dataSet,median_cor) %>%
dplyr::rename(gene_wise_cor=median_cor)) %>%
dplyr::mutate(gene_wise_cor=format_number(gene_wise_cor))
dat <- merge(dat,x$protein_rna$sample_wise_cor_table %>%
dplyr::select(dataSet,median_cor) %>%
dplyr::rename(sample_wise_cor=median_cor))
show_dat <- merge(show_dat,x$protein_rna$sample_wise_cor_table %>%
dplyr::select(dataSet,median_cor) %>%
dplyr::rename(sample_wise_cor=median_cor)) %>%
dplyr::mutate(sample_wise_cor=format_number(sample_wise_cor))
}
return(list(dat=dat,show_dat=show_dat))
}
plot_radar=function(x){
x <- x[,c(1:2,ncol(x):3)]
radar_dat <- x
p <- plot_ly(
type = 'scatterpolargl',
#fill = 'toself',
#alpha = 0.1,
#line = list(smoothing = 1, shape = "spline"),
mode = 'lines+markers'#,
#marker = list(color = "white",
# size = 4,
# line = list(width = 2)),
)
#
for(i in 1:nrow(radar_dat)){
p <- p %>% add_closed_trace(
r = radar_dat[i,-1] %>% as.numeric(),
theta = names(radar_dat)[-1],
#line_close = TRUE,
name = x$dataSet[i]
#line = list(color = "#d3d3d3", dash = "3px")
#mode = "lines"
)
}
p <- p %>%
layout(
polar = list(
radialaxis = list(
visible = T,
range = c(0,1)
)
)
)
return(p)
}
add_closed_trace <- function(p, r, theta, ...){
plotly::add_trace(p, r = c(r, r[1]), theta = c(theta, theta[1]), ...)
}
## Without missing value imputation
noise_signal_analysis=function(x, qc_sample=NULL,bio_sample=NULL, out_dir="./"){
a <- lapply(x, function(y){
dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
as_tibble()
dat$value[dat$value<=0] <- NA
dat$value <- log2(dat$value)
dat <- dat %>% filter(class %in% c(qc_sample, bio_sample))
dat$class <- ifelse(dat$class %in% qc_sample,"QC","Bio")
b <- dat %>% dplyr::group_by(ID,class) %>%
dplyr::summarise(sd=sd(value,na.rm = TRUE),.groups = 'drop') %>%
tidyr::spread(key = class, value = sd)
b <- b %>% dplyr::mutate(nsr=QC/Bio, dataSet=y@ID)
return(b)
})
snr_res <- dplyr::bind_rows(a)
res <- list()
res$table <- snr_res %>% group_by(dataSet) %>%
dplyr::summarise(nsr=median(nsr,na.rm=TRUE),n=n()) %>%
dplyr::arrange(nsr)
res$raw_data <- snr_res
snr_res$dataSet <- factor(snr_res$dataSet,levels = res$table$dataSet)
fig <- paste(out_dir,"/snr.png",sep="")
png(fig,height = 500,width = 700,res=130)
gg <- ggplot(snr_res,aes(x=dataSet,y=nsr,color=dataSet)) +
geom_boxplot(outlier.size = 0.5) +
ylab("Noise to signal ratio")+
theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust=0.5),legend.position = "none")
print(gg)
dev.off()
res$fig <- fig
## MAD
mad_figs <- lapply(x, function(y){
dat <- y@peaksData %>% mutate(dataSet=y@ID,sample=as.character(sample)) %>%
as_tibble()
dat$value[dat$value<=0] <- NA
dat$value <- log2(dat$value)
b <- dat %>% group_by(sample,class) %>%
dplyr::summarise(mad=mad(value,na.rm = TRUE)) %>%
ungroup() %>%
dplyr::arrange(mad)
b$i <- 1:nrow(b)
fig <- paste(out_dir,"/mad",y@ID,".png",sep = "")
png(fig,width = 1000,height = 400,res=150)
gg <- ggplot(b,aes(x=i,y=mad,color=class)) +
geom_point() +
geom_hline(yintercept = median(b$mad)) +
xlab("Rank") +
ylab("MAD (sample)")
print(gg)
dev.off()
return(fig)
})
names(mad_figs) <- names(x)
res$mad_figs <- mad_figs
return(res)
}
get_cv=function(peaksData){
cvstat <- peaksData %>%
dplyr::group_by(class,ID) %>%
dplyr::summarize(cv=sd(value,na.rm=TRUE)/mean(value,na.rm=TRUE)) %>%
dplyr::ungroup() %>%
dplyr::group_by(class) %>%
dplyr::summarize(median_cv=median(cv,na.rm=TRUE),
mean_cv=mean(cv,na.rm=TRUE),
n=n(),
cv30=sum(cv<=0.3,na.rm = TRUE)/n()) %>%
as.data.frame()
return(cvstat)
}
get_cv_table=function(x, metric = "total_features"){
res <- lapply(x,function(y)y[[metric]]) %>% dplyr::bind_rows()
return(res)
}
get_missing_value_ratio=function(x){
n_samples <- length(unique(x$sample))
n_features <- length(unique(x$ID))
n_v <- n_samples * n_features
n_valid_v <- sum(!is.na(x$value))
ratio <- n_valid_v/n_v
return(ratio)
}
get_cv_table=function(x, metric = "total_features"){
res <- lapply(x,function(y)y[[metric]]) %>% dplyr::bind_rows()
return(res)
}
format_overview_table=function(ov_table){
a <- ov_table$dat
show_a <- ov_table$show_dat
nm <- names(ov_table)
for(i in 2:ncol(a)){
y <- a[,i]
y_show <- show_a[,i]
y_show <- cell_spec(y_show, bold = ifelse(y >= max(y), TRUE, FALSE),color = ifelse(y >= max(y), "red", "black"))
show_a[,i] <- y_show
}
show_ov_table <- as.data.frame(t(show_a[,-1]))
names(show_ov_table) <- show_a$dataSet
show_ov_table$metric <- rownames(show_ov_table)
rownames(show_ov_table) <- NULL
show_ov_table <- show_ov_table %>% dplyr::select(metric,everything())
return(show_ov_table)
}
format_number=function(x,digit=4){
format_str <- paste("%.",digit,"f",sep = "")
y <- sprintf(format_str,x)
return(y)
}
get_formated_id_table=function(final_res){
id_table <- get_identification_summary_table(final_res,format = FALSE) %>% dplyr::rename(`data table`=dataSet)
total_features <- final_res$input_parameters$total_features
show_id_table <- id_table
if("#identified features" %in% names(id_table) && "#quantifiable features" %in% names(id_table)){
id_table <- id_table %>% dplyr::mutate(`#identified features`=`#identified features`/total_features,
`#quantifiable features`=`#quantifiable features`/total_features)
show_id_table <- show_id_table %>% dplyr::mutate(`#identified features` = paste(`#identified features`,"\n(",sprintf("%.2f%%",100*id_table$`#identified features`),")",sep=""),
`#quantifiable features` = paste(`#quantifiable features`,"\n(",sprintf("%.2f%%",100*id_table$`#quantifiable features`),")",sep=""))
}
y <- id_table$`#identified features`
show_id_table$`#identified features` <- cell_spec(show_id_table$`#identified features`, bold = ifelse(y >= max(y), TRUE, FALSE),color = ifelse(y >= max(y), "red", "black"))
y <- id_table$`#quantifiable features`
show_id_table$`#quantifiable features` <- cell_spec(show_id_table$`#quantifiable features`, bold = ifelse(y >= max(y), TRUE, FALSE),color = ifelse(y >= max(y), "red", "black"))
return(list(id_table=id_table,show_id_table=show_id_table))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.