Protein phosphatase 2A regulates senescence and immunogenicity in medulloblastoma models.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match

[1] § Methods › In vitro experiments › RNA-Seq. ↔ R/pipeline_functions.R, lines 2929–3031 · score 0.51 · MSigDB, NetBID2, Pathway, regulatory, RNA, genes

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 3,700 lines · 195 KB · Apache-2.0 · 1 match

#' @import Biobase limma tximport igraph biomaRt openxlsx msigdbr ConsensusClusterPlus kableExtra
#' @importFrom GEOquery getGEO
#' @importFrom RColorBrewer brewer.pal
#' @importFrom plot3D scatter3D
#' @importFrom plotrix draw.ellipse draw.circle
#' @importFrom impute impute.knn
#' @importFrom umap umap umap.defaults
#' @importFrom rhdf5 H5Fopen H5Fclose
#' @importFrom DESeq2 DESeqDataSetFromTximport DESeq
#' @importFrom ComplexHeatmap Heatmap
#' @importFrom graphics plot
#' @importFrom aricode clustComp
#' @importFrom GSVA gsva
#' @importFrom MCMCglmm MCMCglmm
#' @importFrom arm bayesglm
#' @importFrom reshape melt
#' @importFrom ordinal clm clmm
#' @importFrom rmarkdown render pandoc_available html_document
#' @importFrom Matrix rowSums
#' @importFrom SummarizedExperiment assay
#' @importFrom lme4 lmer
#' @importFrom grDevices col2rgb colorRampPalette dev.off pdf rgb
#' @importFrom graphics abline arrows axis barplot boxplot hist image layout legend lines mtext par points polygon rect segments strheight stripchart strwidth text
#' @importFrom stats IQR aggregate as.dendrogram as.dist cutree density dist fisher.test gaussian glm hclust kmeans ks.test lm median model.matrix na.omit order.dendrogram p.adjust pchisq pnorm prcomp pt qnorm quantile sd splinefun
#' @importFrom utils read.delim write.table
################################
library(Biobase) ## basic functions for bioconductor
library(GEOquery) ## for samples from GEO
library(limma) ## for data normalization for micro-array
library(DESeq2) ## for data normalization for RNASeq
library(tximport) ## for data import from Salmon/sailfish/kallisto/rsem/stringtie output
library(RColorBrewer) ## for color scale
library(colorspace) ## for color scale
library(plot3D) ## for 3D plot
library(igraph) ## for network related functions
library(plotrix) ## for draw.ellipse
library(biomaRt) ## for gene id conversion
library(openxlsx) ## for output into excel
library(impute) ## for impute
library(msigdbr) ## for msigDB gene sets
library(ComplexHeatmap) ## for complex heatmap
library(umap) ## for umap visualization
library(rhdf5) ## for read in MICA results
library(GSVA)
library(MCMCglmm)
library(arm)
library(reshape)
library(ordinal)
library(rmarkdown)
library(aricode)
##
check_para <- function(para_name,envir){
if(base::exists(para_name,envir=envir)==FALSE){message(sprintf('%s missing !',para_name));return(0)}
if(is.null(base::get(para_name,envir=envir))==TRUE){message(sprintf('%s is NULL !',para_name));return(0)}
return(1)
}
check_option <- function(para_name,option_list,envir){
if(!base::get(para_name,envir=envir) %in% option_list){
message(sprintf('Only accept %s set at: %s !',para_name,base::paste(option_list,collapse=';')));return(0)
}
return(1)
}
clean_charVector <- function(x){
x1 <- names(x)
x <- as.character(x);
x[which(x=='')] <- 'NULL';
x[which(is.null(x)==TRUE)] <- 'NULL'
x[which(is.na(x)==TRUE)] <- 'NA'
names(x) <- x1
x
}
##
#
#' Preload database files into R workspace for NetBID2
#'
#' \code{db.preload} is a pre-processing function for NetBID2. It preloads needed data into R workspace,
#' and saves it locally under db/ directory with specified species name and analysis level.
#'
#' Users need to set the species name (e.g. human, mouse) and
#' analysis level (transcript or gene level). TF list and SIG list are optional, if not specified, list from package data will be used as default.
#'
#' @param use_level character, users can choose "transcript" or "gene". Default is "gene".
#' @param use_spe character, the name of an interested species (e.g. "human", "mouse", "rat"). Default is "human".
#' @param update logical, if TRUE, previous loaded RData will be updated. Default is FALSE.
#' @param TF_list a character vector, the list of TF (Transcription Factor) names. If NULL, the pre-defined list in the package will be used.
#' Default is NULL.
#' @param SIG_list a character vector, the list of SIG (Signaling Factor) names. If NULL, the pre-defined list in the package will be use.
#' Default is NULL.
#' @param input_attr_type character, the type of the TF_list and SIG_list.
#' Details please check biomaRt, \url{https://bioconductor.org/packages/release/bioc/vignettes/biomaRt/inst/doc/biomaRt.html}.
#' If TF_list and SIG_list are not specified, the list in the NetBID2 package will be used.
#' This only support "external_gene_name" and "ensembl_gene_id".
#' Default is "external_gene_name".
#' @param main.dir character, the main directory for NetBID2.
#' If NULL, will be \code{system.file(package = "NetBID2")}. Default is NULL.
#' @param db.dir character, a path for saving the RData.
#' Default is \code{db} directory under the \code{main.dir}, if \code{main.dir} is provided.
#' @param useCache Boolean, parameter pass to \code{getBM()} indicating whether the results cache should be used. Setting to FALSE will disable reading and writing of the cache.
#'
#' @return Return TRUE if loading is successful, otherwise return FALSE. Two variables will be loaded into R workspace, \code{tf_sigs} and \code{db_info}.
#' @examples
#' db.preload(use_level='gene',use_spe='human',update=FALSE)
#'
#' \dontrun{
#' db.preload(use_level='transcript',use_spe='human',update=FALSE)
#' db.preload(use_level='gene',use_spe='mouse',update=FALSE)
#' }
#' @export
db.preload <- function(use_level='transcript',use_spe='human',update = FALSE,
TF_list=NULL,SIG_list=NULL,input_attr_type='external_gene_name',
main.dir=NULL,
db.dir=sprintf("%s/db/",main.dir),useCache = TRUE){
all_input_para <- c('use_level','use_spe','update')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('use_level',c('transcript','gene'),envir=environment()),
check_option('update',c(TRUE,FALSE),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
## load annotation info, including: TF/Sig list, gene info
if(is.null(main.dir)==TRUE){
main.dir <- system.file(package = "NetBID2")
message(sprintf('main.dir not set, will use package directory: %s',main.dir))
}
if(is.null(db.dir)==TRUE){
db.dir <- sprintf("%s/db/",main.dir)
}
message(sprintf('Will use directory %s as the db.dir',db.dir))
message(sprintf('Your setting for species is %s, with level at %s',use_spe,use_level))
use_spe <- toupper(use_spe)
output.db.dir <- sprintf('%s/%s',db.dir,use_spe)
RData.file <- sprintf('%s/%s_%s.RData', output.db.dir,use_spe,use_level)
if (!file.exists(RData.file) | update==TRUE) { ## not exist or need to update
## get info from use_spe
ensembl <- biomaRt::useMart("ensembl")
all_ds <- biomaRt::listDatasets(ensembl)
w1 <- grep(sprintf("^%s GENES",use_spe),toupper(all_ds$description))
if(base::length(w1)==0){
tmp_use_spe <- unlist(strsplit(use_spe,' ')); tmp_use_spe <- tmp_use_spe[base::length(tmp_use_spe)]
w1 <- grep(sprintf(".*%s_GENE_ENSEMBL",toupper(tmp_use_spe)),toupper(all_ds$dataset))
if(base::length(w1)==1){use_spe <- toupper(strsplit(all_ds[w1,2],' ')[[1]][1]); output.db.dir <- sprintf('%s/%s',db.dir,use_spe)}
}
if(base::length(w1)==1){
w2 <- all_ds[w1,1]
mart <- biomaRt::useMart(biomart="ensembl", dataset=w2) ## get id for input spe
message(sprintf('Read in ensembl annotation file for %s and output all db files in %s/%s !',use_spe,db.dir,use_spe))
}
if(base::length(w1)==0){
message(sprintf('Check input use_spe parameter: %s, not included in the ensembl database',use_spe))
return(FALSE)
}
if(base::length(w1)>1){
w2 <- base::paste(all_ds[w1,2],collapse=';')
message(sprintf('Check input use_spe parameter: %s, more than one species match in ensembl database : %s,
please check and re-try',use_spe,w2))
return(FALSE)
}
}
RData.file <- sprintf('%s/%s_%s.RData', output.db.dir,use_spe,use_level)
if (update == TRUE | !file.exists(RData.file)) {
if(!file.exists(output.db.dir)){
dir.create(output.db.dir)
}
## get attributes for mart
filters <- biomaRt::listFilters(mart)
attributes <- biomaRt::listAttributes(mart)
ensembl.attr.transcript <- c('ensembl_transcript_id','ensembl_gene_id',
'external_transcript_name','external_gene_name',
'gene_biotype','gene_biotype',
'chromosome_name','strand','start_position','end_position','band','transcript_start','transcript_end',
'description','phenotype_description','refseq_mrna')
ensembl.attr.gene <- c('ensembl_gene_id','external_gene_name',
'gene_biotype',
'chromosome_name','strand','start_position','end_position','band',
'description','phenotype_description','refseq_mrna')
if(use_spe=='HUMAN'){
ensembl.attr.transcript <- c(ensembl.attr.transcript,'hgnc_symbol','entrezgene_id')
ensembl.attr.gene <- c(ensembl.attr.gene,'hgnc_symbol','entrezgene_id')
}
## do not output hgnc in non-human species
if(use_spe != 'HUMAN')
ensembl.attr.transcript <- base::setdiff(ensembl.attr.transcript,'hgnc_symbol')
if(use_spe != 'HUMAN')
ensembl.attr.gene <- base::setdiff(ensembl.attr.gene,'hgnc_symbol')
## if too much: Query ERROR: caught BioMart::Exception::Usage: Too many attributes selected for External References
# judge input type if TF_list, Sig_list not equal to NULL
if(is.null(TF_list)==FALSE | is.null(SIG_list)==FALSE){
if(!input_attr_type %in% filters$name){
message(sprintf('%s not in the filter name, please retry !',input_attr_type));return(FALSE)
}
if(!input_attr_type %in% ensembl.attr.transcript)
ensembl.attr.transcript <- c(ensembl.attr.transcript,input_attr_type)
if(!input_attr_type %in% ensembl.attr.gene)
ensembl.attr.gene <- c(ensembl.attr.gene,input_attr_type)
}
## get TF/SIG list and output to output.db.dir, if not defined by user, will use in db/ (human)
# for TF list
filter_attr <- input_attr_type
if(is.null(TF_list)){ ## use TF.txt in db/
if(use_spe != 'HUMAN' & use_spe != 'MOUSE'){ ## if spe not human/mouse
TF_f <- sprintf('%s/%s_TF_%s.txt',db.dir,'HUMAN',filter_attr)
message(sprintf('Will use %s file as the input TF_list!',TF_f))
TF_list <- read.delim(TF_f,stringsAsFactors=FALSE,header=F)$V1
filter_attr <- 'external_gene_name'
tmp1 <- biomaRt::getBM(attributes=c('hsapiens_homolog_associated_gene_name','external_gene_name'),values=TRUE,mart=mart,filters='with_hsapiens_homolog',useCache = useCache)
TF_list <- base::unique(tmp1[which(tmp1[,1] %in% TF_list),2])
}else{
TF_f <- sprintf('%s/%s_TF_%s.txt',db.dir,use_spe,filter_attr)
message(sprintf('Will use %s file as the input TF_list!',TF_f))
TF_list <- read.delim(TF_f,stringsAsFactors=FALSE,header=F)$V1
}
}
if(use_level=='transcript'){
message(sprintf('Begin read TF list information from ensembl for %s !',use_spe))
TF_info <- biomaRt::getBM(attributes = ensembl.attr.transcript,values=TF_list, mart=mart, filters=filter_attr,useCache = useCache)
}
if(use_level=='gene'){
message(sprintf('Begin read TF list information from ensembl for %s !',use_spe))
TF_info <- biomaRt::getBM(attributes = ensembl.attr.gene,values=TF_list, mart=mart, filters=filter_attr,useCache = useCache)
}
# for SIG list
if(is.null(SIG_list)){
if(use_spe != 'HUMAN' & use_spe != 'MOUSE'){ ## if spe not human/mouse
SIG_f <- sprintf('%s/%s_SIG_%s.txt',db.dir,'HUMAN',filter_attr)
message(sprintf('Will use %s file as the input SIG_list!',SIG_f))
SIG_list <- read.delim(SIG_f,stringsAsFactors=FALSE,header=F)$V1
filter_attr <- 'external_gene_name'
tmp1 <- biomaRt::getBM(attributes=c('hsapiens_homolog_associated_gene_name','external_gene_name'),values=TRUE,mart=mart,filters='with_hsapiens_homolog',useCache = useCache)
SIG_list <- base::unique(tmp1[which(tmp1[,1] %in% SIG_list),2])
}else{
SIG_f <- sprintf('%s/%s_SIG_%s.txt',db.dir,use_spe,filter_attr)
message(sprintf('Will use %s file as the input SIG_list!',SIG_f))
SIG_list <- read.delim(SIG_f,stringsAsFactors=FALSE,header=F)$V1
}
}
if(use_level=='transcript'){
message(sprintf('Begin read SIG list information from ensembl for %s !',use_spe))
SIG_info <- biomaRt::getBM(attributes = ensembl.attr.transcript,values=SIG_list, mart=mart, filters=filter_attr,useCache = useCache)
}
if(use_level=='gene'){
message(sprintf('Begin read SIG list information from ensembl for %s !',use_spe))
SIG_info <- biomaRt::getBM(attributes = ensembl.attr.gene,values=SIG_list, mart=mart, filters=filter_attr,useCache = useCache)
}
# check input not in the list
miss_TF <- base::unique(base::setdiff(TF_list,TF_info[[filter_attr]]))
miss_SIG <- base::unique(base::setdiff(SIG_list,SIG_info[[filter_attr]]))
if(base::length(miss_TF)>0){message(sprintf("%d TFs could not match,please check and choose to re-try : %s",
base::length(miss_TF),base::paste(sort(miss_TF),collapse=';')))}
if(base::length(miss_SIG)>0){message(sprintf("%d SIGs could not match,please check and choose to re-try : %s",
base::length(miss_SIG),base::paste(sort(miss_SIG),collapse=';')))}
####### output full info
tf_sigs <- list();tf_sigs$tf <- list();tf_sigs$sig <- list();
tf_sigs$tf$info <- TF_info; tf_sigs$sig$info <- SIG_info;
for(each_id_type in base::intersect(c('ensembl_transcript_id','ensembl_gene_id',
'external_transcript_name','external_gene_name','hgnc_symbol',
'entrezgene_id','refseq_mrna'),colnames(TF_info))){
tf_sigs$tf[[each_id_type]] <- base::setdiff(base::unique(TF_info[[each_id_type]]),"")
tf_sigs$sig[[each_id_type]] <- base::setdiff(base::unique(SIG_info[[each_id_type]]),"")
}
db_info <- all_ds[w1,]
save(tf_sigs,db_info=db_info,file = RData.file)
}
load(RData.file,.GlobalEnv)
return(TRUE)
}
#' Get Transcription Factor (TF) and Signaling Factor (SIG) List
#'
#' \code{get.TF_SIG.list} is a function converts gene ID into the corresponding TF/SIG list,
#' with selected gene/transcript type.
#'
#' @param use_genes a vector of characters, genes will be used in the network construction.
#' If NULL, no filter will be performed to the TF/SIG list. Default is NULL.
#' @param use_gene_type character, the attribute name inherited from the biomaRt package.
#' Some options are, "ensembl_gene_id", "ensembl_gene_id_version", "ensembl_transcript_id", "ensembl_transcript_id_version" and "refseq_mrna".
#' All options can be accessed by calling \code{biomaRt::useMart} (e.g. mart <- biomaRt::useMart('ensembl',db_info[1]); biomaRt::listAttributes(mart)$name).
#'
#' The type must match the gene type from the input \code{use_genes}. Default is "external_gene_name".
#' @param ignore_version logical, if TRUE, the version "ensembl_gene_id_version" or "ensembl_transcript_id_version" will be ignored.
#' Default is FALSE.
#' @param dataset character, the dataset used for ID conversion (e.g. "hsapiens_gene_ensembl").
#' If NULL, use \code{db_info[1]} from \code{db.preload}. Default is NULL.
#'
#'
#' @return Return a list containing two elements. \code{tf} is the TF list, \code{sig} is the SIG list.
#'
#' @examples
#' db.preload(use_level='transcript',use_spe='human',update=FALSE)
#' use_genes <- c("ENST00000210187","ENST00000216083","ENST00000216127",
#' "ENST00000216416","ENST00000217233","ENST00000221418",
#' "ENST00000504956","ENST00000507468")
#' res_list <- get.TF_SIG.list(use_gene_type = 'ensembl_transcript_id',
#' use_genes=use_genes,
#' dataset='hsapiens_gene_ensembl')
#' print(res_list)
#'
#' \dontrun{
#' }
#'
#' @export
get.TF_SIG.list <- function(use_genes=NULL,
use_gene_type='external_gene_name',ignore_version=FALSE,
dataset=NULL){
#
all_input_para <- c('use_genes')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('ignore_version',c(TRUE,FALSE),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(dataset)==TRUE){
check_res <- check_para('db_info',envir=environment())
if(base::min(check_res)==0){message('Please use db.preload() to get db_info or set dataset, check and re-try!');return(FALSE)}
dataset <- db_info[1]
}
if(is.null(tf_sigs)==TRUE){
message('tf_sigs not loaded yet, please run db.preload() before processing !');return(FALSE);
}
n1 <- names(tf_sigs$tf)[-1]
if(use_gene_type %in% n1){
if(is.null(use_genes)==TRUE){
TF_list <- base::unique(tf_sigs$tf[[use_gene_type]])
SIG_list <- base::unique(tf_sigs$sig[[use_gene_type]])
}else{
TF_list <- base::unique(base::intersect(use_genes,tf_sigs$tf[[use_gene_type]]))
SIG_list <- base::unique(base::intersect(use_genes,tf_sigs$sig[[use_gene_type]]))
}
}else{
if(grepl('version$',use_gene_type)==TRUE & ignore_version==TRUE){
use_genes_no_v <- gsub('(.*)\\..*','\\1',use_genes)
transfer_tab <- data.frame(to_type=use_genes,from_type=use_genes_no_v,stringsAsFactors = FALSE)
print(str(transfer_tab))
TF_list <- transfer_tab[which(transfer_tab$from_type %in% tf_sigs$tf[[gsub('(.*)_version','\\1',use_gene_type)]]),'to_type']
SIG_list <- transfer_tab[which(transfer_tab$from_type %in% tf_sigs$sig[[gsub('(.*)_version','\\1',use_gene_type)]]),'to_type']
}else{
mart <- biomaRt::useMart(biomart="ensembl", dataset=dataset) ## get mart for id conversion !!!! db_info is saved in db RData
#filters <- biomaRt::listFilters(mart)
attributes <- biomaRt::listAttributes(mart)
if(!use_gene_type %in% attributes$name){
message(sprintf('%s not in the attributes for %s, please check and re-try !',use_gene_type,dataset));return(FALSE)
}
transfer_tab <- get_IDtransfer(from_type=use_gene_type,to_type=n1[1],ignore_version = ignore_version)
print(str(transfer_tab))
TF_list <- get_name_transfertab(use_genes=tf_sigs$tf[n1[1]][[1]],transfer_tab = transfer_tab,ignore_version = ignore_version,ignore_order=TRUE,from_type=n1[1],to_type=use_gene_type)
SIG_list <- get_name_transfertab(use_genes=tf_sigs$sig[n1[1]][[1]],transfer_tab = transfer_tab,ignore_version = ignore_version,ignore_order=TRUE,from_type=n1[1],to_type=use_gene_type)
}
TF_list <- base::unique(base::intersect(use_genes,TF_list))
SIG_list <- base::unique(base::intersect(use_genes,SIG_list))
}
message(sprintf('%d TFs and %s SIGs are included in the expression matrix !',base::length(TF_list),base::length(SIG_list)))
return(list(tf=TF_list,sig=SIG_list))
}
#' Creates Data Frame for ID Conversion
#'
#' \code{get_IDtransfer} creates a data frame for ID conversion using biomaRt. For example, to convert Ensembl ID into gene symbol.
#'
#' @param from_type character, the attribute name match the current ID type (the type of \code{use_genes}).
#' Such as "ensembl_gene_id", "ensembl_gene_id_version", "ensembl_transcript_id", "ensembl_transcript_id_version" or "refseq_mrna".
#' The "attribute" is inherited from the biomaRt package. For details, user can call \code{biomaRt::listAttributes()} to display all available attributes in the selected dataset.
#' @param to_type character, the attribute name to convert into.
#' @param add_type character, the additional attribute name to add into the conversion data frame.
#' @param use_genes a vector of characters, the genes for ID conversion.
#' If NULL, all genes will be selected.
#' @param dataset character, name of the dataset used for ID conversion. For example, "hsapiens_gene_ensembl".
#' If NULL, \code{db_info[1]} will be used. \code{db_info} requires the calling of \code{db.preload} in the previous steps.
#' Default is NULL.
#' @param ignore_version logical, if it is set to TRUE and \code{from_type} is "ensembl_gene_id_version" or "ensembl_transcript_id_version",
#' the version of the original ID will be ignored in ID mapping.
#' Default is FALSE.
#' @param useCache Boolean, parameter pass to \code{getBM()} indicating whether the results cache should be used. Setting to FALSE will disable reading and writing of the cache.
#'
#' @return
#' Return a data frame for ID conversion.
#'
#' @examples
#' use_genes <- c("ENST00000210187","ENST00000216083","ENST00000216127",
#' "ENST00000216416","ENST00000217233","ENST00000221418")
#' transfer_tab <- get_IDtransfer(from_type = 'ensembl_transcript_id',
#' to_type='external_gene_name',
#' use_genes=use_genes,
#' dataset='hsapiens_gene_ensembl')
#' ## get transfer table !!!
#' res1 <- get_name_transfertab(use_genes,transfer_tab=transfer_tab)
#' transfer_tab_withtype <- get_IDtransfer2symbol2type(from_type = 'ensembl_transcript_id',
#' use_genes=use_genes,
#' dataset='hsapiens_gene_ensembl')
#' ## get transfer table !!!
#' \dontrun{
#' }
#' @export
get_IDtransfer <- function(from_type=NULL,to_type=NULL,add_type=NULL,use_genes=NULL,dataset=NULL,ignore_version=FALSE,useCache = TRUE){
#
all_input_para <- c('from_type','to_type')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('ignore_version',c(TRUE,FALSE),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(dataset)==TRUE){
check_res <- check_para('db_info',envir=environment())
if(base::min(check_res)==0){message('Please use db.preload() to get db_info or set dataset, check and re-try!');return(FALSE)}
dataset <- db_info[1]
}
mart <- biomaRt::useMart(biomart="ensembl", dataset=dataset) ## get mart for id conversion !!!! db_info is saved in db RData
attributes <- biomaRt::listAttributes(mart)
if(!from_type %in% attributes$name){
message(sprintf('%s not in the attributes for %s, please check and re-try !',from_type,dataset));return(FALSE)
}
if(!to_type %in% attributes$name){
message(sprintf('%s not in the attributes for %s, please check and re-try !',to_type,dataset));return(FALSE)
}
ori_from_type <- from_type
ori_use_genes <- use_genes
if(from_type %in% c('ensembl_gene_id_version','ensembl_transcript_id_version')){
if(ignore_version==FALSE){
message(sprintf('Attention: %s in %s will be updated with new version number, please check the output.
If lots of missing, try to set ignore_version=TRUE and try again !',from_type,dataset));
}else{
from_type <- gsub('(.*)_version','\\1',from_type)
if(is.null(use_genes)==FALSE) use_genes <- gsub('(.*)\\..*','\\1',use_genes)
}
}
if(is.null(use_genes)==TRUE | base::length(use_genes)>100){
tmp1 <- biomaRt::getBM(attributes=c(from_type,to_type,add_type),values=1,mart=mart,filters='strand',useCache = useCache)
tmp2 <- biomaRt::getBM(attributes=c(from_type,to_type,add_type),values=-1,mart=mart,filters='strand',useCache = useCache)
tmp1 <- base::rbind(tmp1,tmp2)
if(is.null(use_genes)==FALSE){
tmp1 <- tmp1[which(tmp1[,1] %in% use_genes),]
}
}else{
tmp1 <- biomaRt::getBM(attributes=c(from_type,to_type,add_type),values=use_genes,mart=mart,filters=from_type,useCache = useCache)
}
if(ori_from_type %in% c('ensembl_gene_id_version','ensembl_transcript_id_version') & is.null(use_genes)==FALSE & ignore_version==TRUE){
tmp2 <- data.frame(ori_from_type=ori_use_genes,from_type=use_genes,stringsAsFactors=FALSE)
names(tmp2) <- c(ori_from_type,from_type)
tmp1 <- base::merge(tmp2,tmp1,by.y=from_type,by.x=from_type)[c(from_type,to_type,add_type,ori_from_type)]
}
w1 <- apply(tmp1,1,function(x)base::length(which(is.na(x)==TRUE | x=="")))
transfer_tab <- tmp1[which(w1==0),]
for(i in 1:ncol(transfer_tab)){
transfer_tab[,i] <- as.character(transfer_tab[,i])
}
return(transfer_tab)
}
#' Create Data Frame for ID Conversion Between Species
#'
#' \code{get_IDtransfer_betweenSpecies} creates a data frame to convert ID between species.
#'
#' @param from_spe character, name of the original species (e.g. "human", "mouse", "rat") that \code{use_genes} belongs to. Default is "human".
#' @param to_spe character, name of the target species (e.g. "human", "mouse", "rat"). Default is "mouse".
#' @param from_type character, the attribute name match the current ID type (the type of \code{use_genes}).
#' Such as "ensembl_gene_id", "ensembl_gene_id_version", "ensembl_transcript_id", "ensembl_transcript_id_version" or "refseq_mrna".
#' The "attribute" is inherited from the \code{biomaRt} package. For details, user can call \code{biomaRt::listAttributes()} function to display all available attributes in the selected dataset.
#' @param to_type character, the attribute name match the target ID type.
#' @param use_genes a vector of characters, the genes for ID conversion. Must be the genes with ID type of \code{from_type}, and from species \code{from_spe}.
#' If NULL, all the possible genes will be shown in the conversion table. Default is NULL.
#' @param useCache Boolean, parameter pass to \code{getBM()} indicating whether the results cache should be used. Setting to FALSE will disable reading and writing of the cache.
#'
#' @return Return a data frame for ID conversion, from one species to another.
#'
#' @examples
#' use_genes <- c("ENST00000210187","ENST00000216083","ENST00000216127",
#' "ENST00000216416","ENST00000217233","ENST00000221418")
#' transfer_tab <- get_IDtransfer_betweenSpecies(from_spe='human',
#' to_spe='mouse',
#' from_type = 'ensembl_transcript_id',
#' to_type='external_gene_name',
#' use_genes=use_genes)
#' ## get transfer table !!!
#' transfer_tab <- get_IDtransfer_betweenSpecies(from_spe='human',
#' to_spe='mouse',
#' from_type = 'ensembl_transcript_id',
#' to_type='ensembl_transcript_id_version',
#' use_genes=use_genes)
#' ## get transfer table !!!
#' \dontrun{
#' transfer_tab <- get_IDtransfer_betweenSpecies(from_spe='human',
#' to_spe='mouse',
#' from_type='refseq_mrna',
#' to_type='refseq_mrna')
#' }
#' @export
get_IDtransfer_betweenSpecies <- function(from_spe='human',to_spe='mouse',
from_type=NULL,to_type=NULL,
use_genes=NULL,useCache = TRUE){
#
all_input_para <- c('from_spe','to_spe','from_type','to_type')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
from_spe <- toupper(from_spe)
to_spe <- toupper(to_spe)
ensembl <- biomaRt::useMart("ensembl")
all_ds <- biomaRt::listDatasets(ensembl)
w1 <- grep(sprintf("^%s GENES",from_spe),toupper(all_ds$description))
if(base::length(w1)==1){
from_spe_ds <- all_ds[w1,1]
mart1 <- biomaRt::useMart(biomart="ensembl", dataset=from_spe_ds) ## get id for input spe
}
if(base::length(w1)==0){
message(sprintf('Check input from_spe parameter: %s, not included in the ensembl database',from_spe))
return(FALSE)
}
if(base::length(w1)>1){
w2 <- base::paste(all_ds[w1,2],collapse=';')
message(sprintf('Check input from_spe parameter: %s, more than one species match in ensembl database : %s,
please check and re-try',from_spe,w2))
return(FALSE)
}
w1 <- grep(sprintf("^%s GENES",to_spe),toupper(all_ds$description))
if(base::length(w1)==1){
to_spe_ds <- all_ds[w1,1]
mart2 <- biomaRt::useMart(biomart="ensembl", dataset=to_spe_ds) ## get id for input spe
}
if(base::length(w1)==0){
message(sprintf('Check input to_spe parameter: %s, not included in the ensembl database',to_spe))
return(FALSE)
}
if(base::length(w1)>1){
w2 <- base::paste(all_ds[w1,2],collapse=';')
message(sprintf('Check input to_spe parameter: %s, more than one species match in ensembl database : %s,
please check and re-try',to_spe,w2))
return(FALSE)
}
#### mart1 mart2
attributes <- biomaRt::listAttributes(mart1)
if(!from_type %in% attributes$name){
message(sprintf('%s not in the attributes for %s, please check and re-try !',from_type,from_spe));return(FALSE)
}
attributes <- biomaRt::listAttributes(mart2)
if(!to_type %in% attributes$name){
message(sprintf('%s not in the attributes for %s, please check and re-try !',to_type,to_spe));return(FALSE)
}
## get homolog between from_spe to to_spe
cn1 <- gsub('(.*)_gene_ensembl','\\1',from_spe_ds)
cn2 <- attributes$name ## attribute names in mart2
cn3 <- cn2[grep(sprintf('%s_homolog_associated_gene_name',cn1),cn2)]
if(base::length(cn3)!=1){
message('No homolog info found in Biomart, sorry !');return(FALSE)
}
tmp1 <- get_IDtransfer(from_type=from_type,to_type='external_gene_name',use_genes=use_genes,dataset=from_spe_ds)
tmp2 <- biomaRt::getBM(attributes=c(cn3,'external_gene_name'),values=TRUE,
mart=mart2,filters=sprintf('with_%s_homolog',cn1),useCache = useCache)
colnames(tmp1) <- sprintf('%s_%s',colnames(tmp1),from_spe)
colnames(tmp2) <- sprintf('%s_%s',colnames(tmp2),to_spe)
tmp3 <- base::merge(tmp1,tmp2,by.x=sprintf('external_gene_name_%s',from_spe),by.y=sprintf('%s_%s',cn3,to_spe))
transfer_tab <- tmp3[,c(2,3,1)]
if(to_type != 'external_gene_name'){
tmp4 <- get_IDtransfer(from_type='external_gene_name',to_type=to_type,use_genes=tmp3[,3],dataset=to_spe_ds)
colnames(tmp4) <- sprintf('%s_%s',colnames(tmp4),to_spe)
tmp5 <- base::merge(tmp3,tmp4)
transfer_tab <- tmp5[,c(3,4,2,1)]
}
return(transfer_tab)
}
#' Create Data Frame for ID Conversion With Biotype Information
#'
#' \code{get_IDtransfer2symbol2type} creates a data frame to convert original ID into gene symbol and gene biotype (gene level),
#' or into transcript symbol and transcript biotype (transcript level).
#'
#' @param from_type character, the attribute name matches the current ID type (the type of use_genes).
#' Such as "ensembl_gene_id", "ensembl_gene_id_version", "ensembl_transcript_id", "ensembl_transcript_id_version" or "refseq_mrna".
#' The "attribute" is inherited from the biomaRt package.
#' For details, user can call \code{biomaRt::listAttributes()} function to display all available attributes in the selected dataset.
#' @param use_genes a vector of characters, the genes for ID conversion.
#' If NULL, all genes will be selected.
#' @param dataset character, name of the dataset used for ID conversion.
#' For example, "hsapiens_gene_ensembl".
#' If NULL, \code{db_info[1]} will be used. \code{db_info} requires the calling of \code{db.preload} in the previous steps.
#' Default is NULL.
#' @param use_level character, users can chose between "transcript" and "gene". Default is "gene".
#' @param ignore_version logical, if it is set to TRUE and \code{from_type} is "ensembl_gene_id_version" or "ensembl_transcript_id_version",
#' the version of the original ID will be ignored in ID mapping.
#'
#' @return
#' Return a data frame for ID conversion, from ID to gene symbol and gene biotype.
#'
#' @examples
#' use_genes <- c("ENST00000210187","ENST00000216083","ENST00000216127",
#' "ENST00000216416","ENST00000217233","ENST00000221418")
#' transfer_tab <- get_IDtransfer(from_type = 'ensembl_transcript_id',
#' to_type='external_gene_name',use_genes=use_genes,
#' dataset='hsapiens_gene_ensembl')
#' ## get transfer table !!!
#' res1 <- get_name_transfertab(use_genes,transfer_tab=transfer_tab)
#' transfer_tab_withtype <- get_IDtransfer2symbol2type(from_type = 'ensembl_transcript_id',
#' use_genes=use_genes,
#' dataset='hsapiens_gene_ensembl',
#' use_level='transcript')
#' ## get transfer table !!!
#' \dontrun{
#' }
#' @export
get_IDtransfer2symbol2type <- function(from_type=NULL,use_genes=NULL,dataset=NULL,use_level='gene',ignore_version=FALSE){
#
all_input_para <- c('from_type')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('ignore_version',c(TRUE,FALSE),envir=environment()),
check_option('use_level',c('gene','transcript'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(dataset)==TRUE){
check_res <- check_para('db_info',envir=environment())
if(base::min(check_res)==0){message('Please use db.preload() to get db_info or set dataset, check and re-try!');return(FALSE)}
dataset <- db_info[1]
}
message(sprintf('Your setting is at %s level',use_level))
mart <- biomaRt::useMart(biomart="ensembl", dataset=dataset) ## get mart for id conversion !!!! db_info is saved in db RData
attributes <- biomaRt::listAttributes(mart)
if(!from_type %in% attributes$name){
message(sprintf('%s not in the attributes for %s, please check and re-try !',from_type,dataset));return(FALSE)
}
if(use_level=='gene') tmp1 <- get_IDtransfer(from_type=from_type,to_type='external_gene_name',add_type='gene_biotype',use_genes=use_genes,dataset=dataset,ignore_version=ignore_version)
if(use_level=='transcript') tmp1 <- get_IDtransfer(from_type=from_type,to_type='external_transcript_name',add_type='transcript_biotype',use_genes=use_genes,dataset=dataset,ignore_version=ignore_version)
transfer_tab <- tmp1
return(transfer_tab)
}
#' Convert Original Gene ID into Target Gene ID
#'
#' \code{get_name_transfertab} converts the original gene IDs into target gene IDs, with conversion table provided.
#'
#' @param use_genes a vector of characters, the genes for ID conversion.
#' @param transfer_tab data.frame, the conversion table. Users can create it by calling \code{get_IDtransfer}.
#' @param from_type character, the attribute name match the current ID type (the type of use_genes).
#' Such as "ensembl_gene_id", "ensembl_gene_id_version", "ensembl_transcript_id", "ensembl_transcript_id_version" or "refseq_mrna".
#' The "attribute" is inherited from the biomaRt package. For details, user can call \code{biomaRt::listAttributes()} to see all available attributes in the selected dataset.
#' If NULL, will use the first column of \code{transfer_tab}.
#' @param to_type character, the attribute name to convert into. If NULL, will use the second column of \code{transfer_tab}.
#' @param ignore_version logical, if TRUE and \code{from_type} is "ensembl_gene_id_version" or "ensembl_transcript_id_version", the version will be ignored. Default is FALSE.
#' @param ignore_order logical, whether need to ignore the output order to match the input list of \code{use_genes}. Default is FALSE.
#' @return Return a vector of converted gene IDs.
#'
#' @examples
#' use_genes <- c("ENST00000210187","ENST00000216083","ENST00000216127",
#' "ENST00000216416","ENST00000217233","ENST00000221418")
#' transfer_tab <- get_IDtransfer(from_type = 'ensembl_transcript_id',
#' to_type='external_gene_name',use_genes=use_genes,
#' dataset='hsapiens_gene_ensembl')
#' ## get transfer table !!!
#' res1 <- get_name_transfertab(use_genes=use_genes,transfer_tab=transfer_tab)
#' transfer_tab_withtype <- get_IDtransfer2symbol2type(from_type = 'ensembl_transcript_id',
#' use_genes=use_genes,
#' dataset='hsapiens_gene_ensembl')
#' ## get transfer table !!!
#' \dontrun{
#' }
#' @export
get_name_transfertab <- function(use_genes=NULL,transfer_tab=NULL,from_type=NULL,to_type=NULL,ignore_version=FALSE,ignore_order=FALSE){
#
all_input_para <- c('use_genes','transfer_tab')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('ignore_version',c(TRUE,FALSE),envir=environment()),
check_option('ignore_order',c(TRUE,FALSE),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(from_type)==TRUE){from_type=colnames(transfer_tab)[1];}
if(is.null(to_type)==TRUE){to_type=colnames(transfer_tab)[2];}
if(ignore_version==TRUE){
w1 <- which(colnames(transfer_tab)==from_type)
transfer_tab[,w1] <- gsub('(.*)\\..*','\\1',transfer_tab[,w1])
from_type <- gsub('(.*)_version','\\1',from_type)
colnames(transfer_tab)[w1] <- from_type
use_genes <- gsub('(.*)\\..*','\\1',use_genes)
}
transfer_tab <- base::unique(transfer_tab[,c(from_type,to_type)])
x <- use_genes;
t1 <- base::unique(transfer_tab[which(transfer_tab[,from_type] %in% x),])
c1 <- base::unique(t1[,from_type])
if(base::length(c1)<nrow(t1) & ignore_order==FALSE){
message('Gene ID in from type contain multiple items!');return(FALSE)
}
if(ignore_order==TRUE){
x1 <- t1[,to_type]
}else{
rownames(t1) <- t1[,from_type]
x1 <- t1[x,to_type]
w1 <- which(is.na(x1)==TRUE)
x1[w1] <- x[w1]
}
return(x1)
}
#' Manipulation of Working Directories for NetBID2 Network Construction Step
#'
#' \code{NetBID.network.dir.create} is used to help users create an organized working directory for the network construction step in NetBID2 analysis.
#' However, it is not essential for the analysis.
#' It creates a hierarchcial working directory and returns a list contains this directory information.
#'
#' This function needs users to define the main working directory and the project's name.
#' It creates a main working directory with a subdirectory of the project.
#' It also automatically creates three subfolders (QC, DATA and SJAR) within the project folder. QC/,
#' storing Quality Control related plots; DATA/, saving data in RData format;
#' SJAR/, storing files needed for running SJAracne command.
#' This function also returns a list object (example, \code{network.par} in the demo) with directory information wrapped inside.
#' This list is an essential for
#' network construction step, all the important intermediate data generated later will be wrapped inside.
#' @param project_main_dir character, name or absolute path of the main working directory.
#' @param project_name character, name of the project folder.
#'
#' @return \code{NetBID.network.dir.create} returns a list object, containing main.dir (path of the main working directory),
#' project.name (project name), out.dir (path of the project folder, which contains three subfolders), out.dir.QC,
#' out.dir.DATA and out.dir.SJAR.
#' @examples
#'
#' \dontrun{
#' # Creating a main working directory under the current working directory by folder name
#' network.par <- NetBID.network.dir.create("MyMainDir","MyProject")
#' # Or creating a main working directory under the current working directory by relative path
#' network.par <- NetBID.network.dir.create("./MyMainDir","MyProject")
#' # Or creating a main working directory to a specific path by absolute path
#' network.par <- NetBID.network.dir.create("~/Desktop/MyMainDir","MyProject")
#' }
#' @export
NetBID.network.dir.create <- function(project_main_dir=NULL,project_name=NULL){
#
if(base::exists('network.par')==TRUE){
stop('network.par is occupied in the current session,please manually run: rm(network.par) and re-try, otherwise will not change !');
}
#
all_input_para <- c('project_main_dir','project_name')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
network.par <- list()
network.par$main.dir <- project_main_dir
network.par$project.name <- project_name
network.par$out.dir <- sprintf('%s/%s',network.par$main.dir,network.par$project.name)
# create output directory
if (!dir.exists(project_main_dir)) {
dir.create(project_main_dir, recursive = TRUE)
}
if (!dir.exists(network.par$out.dir)) {
dir.create(network.par$out.dir, recursive = TRUE)
}
network.par$out.dir.QC <- paste0(network.par$out.dir, '/QC/')
if (!dir.exists(network.par$out.dir.QC)) {
dir.create(network.par$out.dir.QC, recursive = TRUE) ## directory for QC
}
network.par$out.dir.DATA <- paste0(network.par$out.dir, '/DATA/')
if (!dir.exists(network.par$out.dir.DATA)) {
dir.create(network.par$out.dir.DATA, recursive = TRUE) ## directory for DATA
}
network.par$out.dir.SJAR <- paste0(network.par$out.dir, '/SJAR/')
if (!dir.exists(network.par$out.dir.SJAR)) {
dir.create(network.par$out.dir.SJAR, recursive = TRUE) ## directory for SJARAcne
}
message(sprintf('Project space created, please check %s',network.par$out.dir))
return(network.par)
}
#' Manipulation of Working Directories for NetBID2 Driver Estimation Step
#'
#' \code{NetBID.analysis.dir.create} is used to help users create an organized working directory
#' for the driver estimation step in NetBID2 analysis.
#' However, it is not essential for the analysis.
#' It creates a hierarchcial working directory and returns a list contains this directory information.
#'
#' This function requires user to define the main working directory and the project’s name.
#' It creates a main working directory with a subdirectory of the project.
#' It also automatically creates three subfolders (QC, DATA and PLOT) within the project folder.
#' QC/, storing Quality Control related plots; DATA/, saving data in RData format; PLOT/, storing output plots.
#' This function also returns a list object (e.g. \code{analysis.par} in the demo) with directory information wrapped inside.
#' This list is an essential for driver construction step, all the important intermediate data generated later will be wrapped inside.
#'
#' @param project_main_dir character, name or absolute path of the main working directory for driver analysis.
#' @param project_name character, name of the project folder.
#' @param network_dir character, name or absolute path of the main working directory for network construction.
#' @param network_project_name character, the project name of network construction. Or use the project name of SJARACNe.
#' This parameter is optional. If one didn't run NetBID2 network construction part in the pipeline, he could set it to NULL.
#' If one like to follow the NetBID2 pipeline, he should set it to the path of the TF network file and the SIG network file.
#' @param tf.network.file character, the path of the TF network file (e.g. "XXX/consensus_network_ncol_.txt").
#' Default is the path of network_project_name.
#' @param sig.network.file character, the path of the SIG network file (e.g. "XXX/consensus_network_ncol_.txt").
#' Default is the path of network_project_name.
#'
#' @return Returns a list object, containing main.dir (path of the main working directory), project.name (project name),
#' out.dir (path of the project folder, which contains three subfolders), out.dir.QC, out.dir.DATA and out.dir.PLOT.
#'
#' @examples
#'
#' \dontrun{
#' network.dir <- sprintf('%s/demo1/network/',system.file(package = "NetBID2")) # use demo
#' network.project.name <- 'project_2019-02-14' #
#' project_main_dir <- 'demo1/'
#' project_name <- 'driver_test'
#' analysis.par <- NetBID.analysis.dir.create(project_main_dir=project_main_dir,
#' project_name=project_name,
#' network_dir=network.dir,
#' network_project_name=network.project.name)
#' }
#' @export
NetBID.analysis.dir.create <- function(project_main_dir=NULL,project_name=NULL,
network_dir=NULL,
network_project_name=NULL,
tf.network.file=NULL,
sig.network.file=NULL){
#
if(base::exists('analysis.par')==TRUE){
stop('analysis.par is occupied in the current session,please manually run: rm(analysis.par) and re-try, otherwise will not change !');
}
#
all_input_para <- c('project_main_dir','project_name')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
if((is.null(network_dir)==TRUE | is.null(network_project_name)==TRUE) & (is.null(tf.network.file)==TRUE | is.null(sig.network.file)==TRUE)){
message('Either network_dir,network_project_name or tf.network.file,sig.network.file is required, please check and re-try !')
return(FALSE);
}
#
analysis.par <- list()
analysis.par$main.dir <- project_main_dir
analysis.par$project.name <- project_name
analysis.par$out.dir <- sprintf('%s/%s/',analysis.par$main.dir,analysis.par$project.name)
analysis.par$tf.network.file <- ''
analysis.par$sig.network.file <- ''
if(is.null(tf.network.file)==FALSE){
analysis.par$tf.network.file <- tf.network.file
}else{
tf_net1 <- sprintf('%s/SJAR/%s/output_tf_sjaracne_%s_out_.final/consensus_network_ncol_.txt',
network_dir,network_project_name,network_project_name) ## old version of sjaracne
tf_net2 <- sprintf('%s/SJAR/SJARACNE_%s_TF/consensus_network_ncol_.txt',
network_dir,network_project_name) ## new version of sjaracne
if(file.exists(tf_net2)) analysis.par$tf.network.file <- tf_net2 else analysis.par$tf.network.file <- tf_net1
}
if(is.null(tf.network.file)==FALSE){
analysis.par$sig.network.file <- sig.network.file
}else{
sig_net1 <- sprintf('%s/SJAR/%s/output_sig_sjaracne_%s_out_.final/consensus_network_ncol_.txt',
network_dir,network_project_name,network_project_name) ## old version of sjaracne
sig_net2 <- sprintf('%s/SJAR/SJARACNE_%s_SIG/consensus_network_ncol_.txt',
network_dir,network_project_name) ## new version of sjaracne
if(file.exists(sig_net2)) analysis.par$sig.network.file <- sig_net2 else analysis.par$sig.network.file <- sig_net1
}
if(file.exists(analysis.par$tf.network.file)){
message(sprintf('TF network file found in %s',analysis.par$tf.network.file))
}else{
message(sprintf('TF network file not found in %s, please check and re-try !',analysis.par$tf.network.file))
return(FALSE)
}
if(file.exists(analysis.par$sig.network.file)){
message(sprintf('SIG network file found in %s',analysis.par$sig.network.file))
}else{
message(sprintf('SIG network file not found in %s, please check and re-try ',analysis.par$sig.network.file))
return(FALSE)
}
# create output directory
if (!dir.exists(analysis.par$out.dir)) {
dir.create(analysis.par$out.dir, recursive = TRUE)
}
analysis.par$out.dir.QC <- paste0(analysis.par$out.dir, '/QC/')
if (!dir.exists(analysis.par$out.dir.QC)) {
dir.create(analysis.par$out.dir.QC, recursive = TRUE) ## directory for QC
}
analysis.par$out.dir.DATA <- paste0(analysis.par$out.dir, '/DATA/')
if (!dir.exists(analysis.par$out.dir.DATA)) {
dir.create(analysis.par$out.dir.DATA, recursive = TRUE) ## directory for DATA
}
analysis.par$out.dir.PLOT <- paste0(analysis.par$out.dir, '/PLOT/')
if (!dir.exists(analysis.par$out.dir.PLOT)) {
dir.create(analysis.par$out.dir.PLOT, recursive = TRUE) ## directory for Result Plots
}
#
message(sprintf('Analysis space created, please check %s',analysis.par$out.dir))
return(analysis.par)
}
#' Save Data Produced by Corresponding NetBID2 Pipeline Step.
#'
#' \code{NetBID.saveRData} is a function to save complicated list object generated by certain steps of NetBID2's pipeline
#' (e.g. load gene expression file from GEO, 'exp-load').
#' This function is not essential, but it is highly suggested for easier pipeline step checkout and reference.
#'
#' There are two important steps in the NetBID2 pipeline, network construction and driver analysis.
#' User can save these two complicated list objects, network.par and analysis.par.
#' Assigning the \code{step} name to save the RData for easier reference.
#' Calling \code{NetBID.loadRData} to load the corresponding step RData, users can avoid repeating the former steps.
#'
#' @param network.par list, stores all related datasets from network construction pipeline step.
#' @param analysis.par list, stores all related datasets from driver analysis pipeline step.
#' @param step character, name of the pipeline step decided by user for easier reference.
#'
#' @examples
#' \dontrun{
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' NetBID.saveRData(analysis.par=analysis.par,step='ms-tab_test')
#' }
#' @export
NetBID.saveRData <- function(network.par=NULL,analysis.par=NULL,step='exp-load'){
#
all_input_para <- c('step')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(network.par)==FALSE & is.null(analysis.par)==FALSE){
message('Can not save network.par and analysis.par at once, please only use one !');return(FALSE)
}
if(is.null(network.par)==FALSE){
save(network.par,file=sprintf('%s/network.par.Step.%s.RData',network.par$out.dir.DATA,step))
message(sprintf('Successful save to %s',sprintf('%s/network.par.Step.%s.RData',network.par$out.dir.DATA,step)))
}
if(is.null(analysis.par)==FALSE){
save(analysis.par,file=sprintf('%s/analysis.par.Step.%s.RData',analysis.par$out.dir.DATA,step))
message(sprintf('Successful save to %s',sprintf('%s/analysis.par.Step.%s.RData',analysis.par$out.dir.DATA,step)))
}
}
#' Reload Saved RData Created by \code{NetBID.saveRData}.
#'
#' \code{NetBID.loadRData} is a function loads RData saved by \code{NetBID.saveRData} function.
#' It prevents user from repeating former pipeline steps.
#'
#' @param network.par list, stores all related datasets from network construction step.
#' @param analysis.par list, stores all related datasets from driver analysis step.
#' @param step character, name of the pipeline step. It should be previously assigned by user when calling \code{NetBID.saveRData} function.
#'
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#'
#' @export
NetBID.loadRData <- function(network.par=NULL,analysis.par=NULL,step='exp-load'){
#
all_input_para <- c('step')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(network.par)==FALSE & is.null(analysis.par)==FALSE){
message('Can not load network.par and analysis.par at once, please only use one !');return(FALSE)
}
if(is.null(network.par)==FALSE){
load(file=sprintf('%s/network.par.Step.%s.RData',network.par$out.dir.DATA,step),.GlobalEnv)
message(sprintf('Successful load from %s',sprintf('%s/network.par.Step.%s.RData',network.par$out.dir.DATA,step)))
}
if(is.null(analysis.par)==FALSE){
load(file=sprintf('%s/analysis.par.Step.%s.RData',analysis.par$out.dir.DATA,step),.GlobalEnv)
message(sprintf('Successful load from %s',sprintf('%s/analysis.par.Step.%s.RData',analysis.par$out.dir.DATA,step)))
}
}
#' Download Gene Expression Series From GEO Database with Platform Specified
#'
#' \code{load.exp.GEO} downloads user assigned Gene Expression Series (GSE file) along with its Platform from GEO dataset.
#' It returns an ExpressionSet class object and saves it as RData. If the GSE RData already exists, it will be loaded directly.
#' It also allows users to update the Gene Expression Series RData saved before.
#'
#' @param out.dir character, the file path used to save the GSE RData. If the data already exsits, it will be loaded from this path.
#' @param GSE character, the GEO Series Accession ID.
#' @param GPL character, the GEO Platform Accession ID.
#' @param getGPL logical, if TRUE, the corresponding GPL file will be downloaded. Default is TRUE.
#' @param update logical, if TRUE, the previous stored Gene ExpressionSet RData will be updated. Default is FALSE
#'
#' @return Return an ExpressionSet class object.
#' @examples
#'
#' \dontrun{
#' # Download the GSE116028 which performed on GPL6480 platform
#' # from GEO and save it to the current directory
#' # Assign this ExpressionSet object to net_eset
#' net_eset <- load.exp.GEO(out.dir='./',
#' GSE='GSE116028',
#' GPL='GPL6480',
#' getGPL=TRUE,
#' update=FALSE)
#' }
#' @export
load.exp.GEO <- function(out.dir = NULL,GSE = NULL,GPL = NULL,getGPL=TRUE,update = FALSE){
#
all_input_para <- c('out.dir','GSE','GPL')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('getGPL',c(TRUE,FALSE),envir=environment()),
check_option('update',c(TRUE,FALSE),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(!grepl('^GSE',GSE)){
message('Only support GSE ID')
return(FALSE)
}
expRData_dir <- sprintf('%s/%s_%s.RData', out.dir, GSE,GPL)
if (file.exists(expRData_dir) & update == FALSE) {
message(sprintf('RData exist in %s and update==TRUE, will directly load from RData .',expRData_dir))
load(expRData_dir)
} else{
eset <- GEOquery::getGEO(GSE, GSEMatrix = TRUE, getGPL = getGPL)
if (base::length(eset) > 1)
idx <- grep(GPL, attr(eset, "names"))
else
idx <- 1
eset <- eset[[idx]]
if(GPL!=annotation(eset)) {GPL <- annotation(eset); message(sprintf('Real GPL:%s',GPL))}
expRData_dir <- sprintf('%s/%s_%s.RData', out.dir, GSE,GPL)
save(eset, file = expRData_dir)
message(sprintf('RData for the eset is saved in %s .',expRData_dir))
}
return(eset)
}
#' Load Gene Expression Set from Salmon Output (demo version)
#'
#' \code{load.exp.RNASeq.demoSalmon} is a function to read in Salmon results and convert it to eSet/DESeqDataSet class object.
#'
#' This function helps users to read in results created by Salmon.
#' Due to the complicated manipulations (e.g. reference sequence) in processing Salmon, this demo function may not be suitable for all scenarios.
#'
#' @param salmon_dir character, the directory to save the results created by Salmon.
#' @param tx2gene data.frame or NULL, this parameter will be passed to \code{tximport}. For details, please check \code{tximport}.
#' If NULL, will read in one of the transcript names from Salmon's results. Note, it works when using e.g. "gencode.v29.transcripts.fa" from GENCODE as reference.
#' @param use_phenotype_info data.frame, the data frame contains phenotype information. It must have the columns \code{use_sample_col} and \code{use_design_col}.
#' @param use_sample_col character, the column name, indicating which column in \code{use_phenotype_info} should be used as the sample name.
#' @param use_design_col character, the column name, indicating which column in \code{use_phenotype_info} should be used as the design feature for samples.
#' @param return_type character, the class of the return object.
#' "txi" is the output of tximport. It is a list containing three matrices, abundance, counts and length.
#' "counts" is the matrix of raw count.
#' "tpm" is the raw tpm.
#' "fpm", "cpm" is the fragments/counts per million mapped fragments (fpm/cpm).
#' "raw-dds" is the DESeqDataSet class object, which is the original one without processing.
#' "dds" is the DESeqDataSet class object, which is processed by DESeq.
#' "eset" is the ExpressionSet class object, which is processed by DESeq and vst.
#' Default is "tpm".
#' @param merge_level character, users can choose between "gene" and "transcript".
#' "gene", the original salmon results will be mapped to the transcriptome and the expression matrix will be merged to the gene level.
#' This only works when using e.g. "gencode.vXX.transcripts.fa" from GENCODE as the reference.
#' @export
load.exp.RNASeq.demoSalmon <- function(salmon_dir = NULL,tx2gene=NULL,
use_phenotype_info = NULL,
use_sample_col=NULL,
use_design_col=NULL,
return_type='tpm',
merge_level='gene') {
#
all_input_para <- c('salmon_dir','use_phenotype_info','use_sample_col','use_design_col','return_type','merge_level')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('return_type',c('txi','counts','tpm','fpm','cpm','raw-dds','dds','eset'),envir=environment()),
check_option('merge_level',c('gene','transcript'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
files <- file.path(salmon_dir, list.files(salmon_dir), "quant.sf")
if(!grepl('_salmon',use_phenotype_info[1,use_sample_col])) sample_name <- gsub('(.*)_salmon', '\\1', list.files(salmon_dir)) ## if no _salmon is also ok
names(files) <- sample_name
w1 <- base::length(files)
message(sprintf('%d %s/quant.sf found !',w1,salmon_dir))
if(is.null(tx2gene)){
gene_info <- read.delim(file = files[1], stringsAsFactors = FALSE)[, 1]
gen1 <- sapply(gene_info, function(x)unlist(strsplit(x, '\\|')))
gen1 <- t(gen1)
if(merge_level=='gene'){
tx2gene <- data.frame('transcript' = gene_info,'gene' = gen1[,2],stringsAsFactors = FALSE)
}else{
tx2gene <- data.frame('transcript' = gene_info,'gene' = gen1[,1],stringsAsFactors = FALSE)
}
}
eset <- load.exp.RNASeq.demo(files,type='salmon',
tx2gene=tx2gene,
use_phenotype_info=use_phenotype_info,
use_sample_col=use_sample_col,
use_design_col=use_design_col,
return_type=return_type,
merge_level=merge_level)
return(eset)
}
#' Load Gene Expression Set from RNA-Seq Results (demo version)
#'
#' \code{load.exp.RNASeq.demo} is a function to read in RNA-Seq results and convert it to \code{eSet/DESeqDataSet} class object.
#'
#' This function helps users to read in RNA-Seq results from various sources.
#' Due to the complicated manipulations (e.g. reference sequence) in processing RNA-Seq, this demo function may not be suitable for all scenarios.
#'
#' @param files a vector of characters, the filenames for the transcript-level abundances. It will be passed to \code{tximport}.
#' For details, please check \code{tximport}.
#' @param type character, the type of software used to generate the abundances. It will be passed to \code{tximport}.
#' For details, please check \code{tximport}.
#' @param tx2gene data.frame or NULL, this parameter will be passed to \code{tximport}. For details, please check \code{tximport}.
#' @param use_phenotype_info data.frame, the data frame contains phenotype information. It must have the columns \code{use_sample_col} and \code{use_design_col}.
#' @param use_sample_col character, the column name, indicating which column in \code{use_phenotype_info} should be used as the sample name.
#' @param use_design_col character, the column name, indicating which column in \code{use_phenotype_info} should be used as the design feature for samples.
#' @param return_type character, the class of the return object.
#' "txi" is the output of \code{tximport}. It is a list containing three matrices, abundance, counts and length.
#' "counts" is the matrix of raw count.
#' "tpm" is the raw tpm.
#' "fpm", "cpm" is the fragments/counts per million mapped fragments.
#' "raw-dds" is the DESeqDataSet class object, which is the original one without processing.
#' "dds" is the DESeqDataSet class object, which is processed by \code{DESeq}.
#' "eset" is the ExpressionSet class object, which is processed by \code{DESeq} and \code{vst}.
#' Default is "tpm".
#' @param merge_level character, users can choose between "gene" and "transcript".
#' "gene", the original salmon results will be mapped to the transcriptome and the expression matrix will be merged to the gene level.
#' This only works when using e.g. "gencode.vXX.transcripts.fa" from GENCODE as the reference.
#' @export
load.exp.RNASeq.demo <- function(files,type='salmon',
tx2gene=NULL,
use_phenotype_info = NULL,
use_sample_col=NULL,
use_design_col=NULL,
return_type='tpm',
merge_level='gene') {
#
all_input_para <- c('files','type','tx2gene','use_phenotype_info','use_sample_col','use_design_col','return_type','merge_level')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('return_type',c('txi','counts','tpm','fpm','cpm','raw-dds','dds','eset'),envir=environment()),
check_option('merge_level',c('gene','transcript'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
n1 <- colnames(use_phenotype_info)
if(!use_sample_col %in% n1){
message(sprintf('%s not in the colnames of use_phenotype_info,
please check and re-try !',use_sample_col));return(FALSE)
}
if(!use_design_col %in% n1){
message(sprintf('%s not in the colnames of use_phenotype_info,
please check and re-try !',use_design_col));return(FALSE)
}
# get intersected samples
rownames(use_phenotype_info) <- use_phenotype_info[,use_sample_col]
w1 <- base::intersect(names(files),rownames(use_phenotype_info))
files <- files[w1]; use_phenotype_info <- use_phenotype_info[w1,]
if(base::length(w1)==0){
message(sprintf('No sample could match the %s in the use_phenotype_info, please check and re-try !',use_sample_col))
return(FALSE)
}
message(sprintf('%d samples could further processed !',base::length(w1)))
# import into txi
txi <- tximport::tximport(files, type = type, tx2gene = tx2gene) ## key step one, tximport
if(return_type=='counts'){
return(txi$counts)
}
if(return_type=='tpm'){
return(txi$abundance)
}
if(return_type=='txi'){
return(txi)
}
use_phenotype_info <- use_phenotype_info[colnames(txi$abundance), ]
tmp_phe <- base::cbind(group=use_phenotype_info[,use_design_col],use_phenotype_info,stringsAsFactors=FALSE)
# import into deseq2
dds <- DESeq2::DESeqDataSetFromTximport(txi, colData = tmp_phe, design = ~ group) ## key step two, DESeqDataSetFromTximport
if(return_type=='raw-dds'){
return(dds)
}
if(return_type=='fpm' | return_type=='cpm'){
return(DESeq2::fpm(dds))
}
dds <- DESeq2::DESeq(dds)
if(return_type=='dds'){
return(dds)
}else{
vsd <- DESeq2::vst(dds)
mat <- SummarizedExperiment::assay(vsd)
eset <- generate.eset(exp_mat=mat, phenotype_info = use_phenotype_info, feature_info = NULL, annotation_info='Salmon')
if(return_type=='eset') return(eset)
if(return_type=='both') return(list(eset=eset,dds=dds))
}
}
#' Generate ExpressionSet Object
#'
#' \code{generate.eset} generates ExpressionSet class object to contain and describe the high-throughput assays.
#' Users need to define its slots, which are expression matrix (required),
#' phenotype information and feature information (optional).
#' It is very useful when only expression matrix is available.
#'
#' @param exp_mat matrix, the expression data matrix. Each row represents a gene/transcript/probe, each column represents a sample.
#' @param phenotype_info data.frame, the phenotype information for all the samples in \code{exp_mat}.
#' In the phenotype data frame, each row represents a sample, each column represents a phenotype feature.
#' The row names must match the column names of \code{exp_mat}. If NULL, it will generate a single-column data frame.
#' Default is NULL.
#' @param feature_info data.frame, the feature information for all the genes/transcripts/probes in \code{exp_mat}.
#' In the feature data frame, each row represents a gene/transcript/probe and each column represents an annotation of the feature.
#' The row names must match the row names of \code{exp_mat}. If NULL, it will generate a single-column data frame.
#' Default is NULL.
#' @param annotation_info character, the annotation set by users for easier reference. Default is "".
#'
#' @return Return an ExressionSet object.
#'
#' @examples
#' mat1 <- matrix(rnorm(10000),nrow=1000,ncol=10)
#' colnames(mat1) <- paste0('Sample',1:ncol(mat1))
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' eset <- generate.eset(exp_mat=mat1)
#' @export
generate.eset <- function(exp_mat=NULL, phenotype_info=NULL, feature_info=NULL, annotation_info="") {
#
all_input_para <- c('exp_mat')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(dim(exp_mat))==TRUE){
exp_mat <- t(as.matrix(exp_mat));rownames(exp_mat) <- 'g1'
}
if (is.null(phenotype_info)) {
phenotype_info <- data.frame(group = colnames(exp_mat), stringsAsFactors = FALSE)
rownames(phenotype_info) <- colnames(exp_mat)
}
if (is.null(feature_info)) {
feature_info <- data.frame(gene = rownames(exp_mat), stringsAsFactors = FALSE)
rownames(feature_info) <- rownames(exp_mat)
}
if((class(phenotype_info)=='character' | is.null(dim(phenotype_info))==TRUE) & is.null(names(phenotype_info))==TRUE){
phenotype_info <- data.frame(group = phenotype_info, stringsAsFactors = FALSE)
rownames(phenotype_info) <- colnames(exp_mat)
}
if((class(feature_info)=='character' | is.null(dim(feature_info))==TRUE) & is.null(names(feature_info))==TRUE){
feature_info <- data.frame(gene = feature_info, stringsAsFactors = FALSE)
rownames(feature_info) <- rownames(exp_mat)
}
#
eset <-
new(
"ExpressionSet",
phenoData = new("AnnotatedDataFrame", phenotype_info),
featureData = new("AnnotatedDataFrame", feature_info),
annotation = annotation_info,
exprs = as.matrix(exp_mat)
)
return(eset)
}
#' Merge Two ExpressionSet Class Objects into One
#'
#' \code{merge_eset} merges two ExpressionSet class objects and returns one ExpresssionSet object.
#' If genes in the two ExpressionSet objects are identical, the expression matrix will be combined directly.
#' Otherwise, Z-transformation is strongly suggested to be performed before combination (set std=TRUE).
#'
#' @param eset1 ExpressionSet class, the first ExpressionSet.
#' @param eset2 ExpressionSet class, the second ExpressionSet.
#' @param group1 character, name of the first ExpressionSet.
#' @param group2 character, name of the second ExpressionSet.
#' @param use_col a vector of characters, the column names in the phenotype information to be kept.
#' If NULL, shared column names of \code{eset1} and \code{eset2} will be used. Default is NULL.
#' @param group_col_name character, name of the column which contains the names defined in \code{group1} and \code{group2}.
#' This column is designed to show which original ExpressionSet each sample comes from before combination.
#' Default name of this column is "original_group".
#' @param remove_batch logical, if TRUE, remove the batch effects from these two expression datasets. Default is FALSE.
#' @param std logical, whether to perform std to the original expression matrix. Default is FALSE.
#'
#' @return Return an ExressionSet class object.
#' @examples
#' mat1 <- matrix(rnorm(10000),nrow=1000,ncol=10)
#' colnames(mat1) <- paste0('Sample1_',1:ncol(mat1))
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' eset1 <- generate.eset(exp_mat=mat1)
#' mat2 <- matrix(rnorm(10000),nrow=1000,ncol=10)
#' colnames(mat2) <- paste0('Sample2_',1:ncol(mat1))
#' rownames(mat2) <- paste0('Gene',1:nrow(mat1))
#' eset2 <- generate.eset(exp_mat=mat2)
#' new_eset <- merge_eset(eset1,eset2)
#' @export
merge_eset <- function(eset1,eset2,
group1=NULL,group2=NULL,
group_col_name='original_group',
use_col = NULL,
remove_batch = FALSE,std=FALSE) {
#
all_input_para <- c('eset1','eset2','group_col_name')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('remove_batch',c(TRUE,FALSE),envir=environment()),
check_option('std',c(TRUE,FALSE),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
mat1 <- Biobase::exprs(eset1)
mat2 <- Biobase::exprs(eset2)
w1 <- base::intersect(rownames(mat1), rownames(mat2))
if(base::length(w1)==0){
message('No overlap genes between two eSet, please check and re-try!');return(FALSE);
}
if((base::length(w1)<nrow(mat1) | base::length(w1)<nrow(mat2)) & std==FALSE){
message('Original two esets contain different gene list, strongly suggest to do z transformation (set std=TRUE) across all samples before merge!');
}
if(std==TRUE){
## z-transformation
mat1 <- apply(mat1,2,do.std) # std to samples
mat2 <- apply(mat2,2,do.std)
}
rmat <- base::cbind(as.data.frame(mat1)[w1, ], as.data.frame(mat2)[w1,])
rmat <- as.matrix(rmat)
#choose1 <- apply(rmat <= quantile(rmat, probs = 0.05), 1, sum) <= ncol(rmat) * 0.90 ## low expressed genes
#rmat <- rmat[choose1, ]
phe1 <- Biobase::pData(eset1)
phe2 <- Biobase::pData(eset2)
phe1 <- as.data.frame(apply(phe1,2,clean_charVector),stringsAsFactors=F)
phe2 <- as.data.frame(apply(phe2,2,clean_charVector),stringsAsFactors=F)
if(base::length(use_col)==0){
use_col <- base::intersect(colnames(phe1),colnames(phe2))
}
rphe <- list();
if(base::length(use_col)>1)
rphe <- base::rbind(phe1[colnames(mat1), use_col], phe2[colnames(mat2), use_col])
if(base::length(use_col)==1){
rphe <- c(phe1[colnames(mat1), use_col], phe2[colnames(mat2), use_col])
rphe <- data.frame(rphe,stringsAsFactors=FALSE); colnames(rphe) <- use_col;
rownames(rphe) <- colnames(rmat)
}
if(base::length(use_col)==0){message('Warning: no intersected phenotype column!');}
if(is.null(group1)==TRUE) group1 <- 'group1'
if(is.null(group2)==TRUE) group2 <- 'group2'
rphe[[group_col_name]]<- c(rep(group1, ncol(mat1)), rep(group2, ncol(mat2)))
if (remove_batch == TRUE) {
rmat <- limma::removeBatchEffect(rmat,batch=rphe[[group_col_name]])
}
if(class(rphe)=='list'){rphe <- as.data.frame(rphe,stringsAsFactors=FALSE); rownames(rphe) <- colnames(rmat)}
reset <- generate.eset(rmat,phenotype_info = rphe, annotation_info = 'combine')
return(reset)
}
#' Reassign featureData slot of ExpressionSet and Update feature information
#'
#' \code{update_eset.feature} reassigns the featureData slot of ExpressionSet object based on user's demand. It is mainly used for gene ID conversion.
#'
#' User can pass a conversion table to \code{use_eset} for the ID conversion. A conversion table can be obtained from the original featureDta slot
#' (if one called the \code{load.exp.GEO} function and set getGPL==TRUE) or by running the \code{get_IDtransfer} function.
#' The mapping between original ID and target ID can be summerised into 4 categories.
#' 1) One-to-one, simply replaces the original ID with target ID;
#' 2) Many-to-one, the expression value for the target ID will be merged from its original ID;
#' 3) One-to-many, the expression value for the original ID will be distributed to the matched target IDs;
#' 4) Many-to-many, apply part 3) first, then part 2).
#'
#' @param use_eset ExpressionSet class object.
#' @param use_feature_info data.frame, a data frame contains feature information, it can be obtained by calling \code{fData} function.
#' @param from_feature character, orginal ID. Must be one of the column names in \code{use_feature_info} and correctly characterize the \code{use_eset}'s row names.
#' @param to_feature character, target ID. Must be one of the column names in \code{use_feature_info}.
#' @param merge_method character, the agglomeration method to be used for merging gene expression value.
#' This should be one of, "median", "mean", "max" or "min". Default is "median".
#' @param distribute_method character, the agglomeration method to be used for distributing the gene expression value.
#' This should be one of, "mean" or "equal". Default is "equal".
#'
#' @return Return an ExressionSet object with updated feature information.
#' @examples
#' mat1 <- matrix(rnorm(10000),nrow=1000,ncol=10)
#' colnames(mat1) <- paste0('Sample',1:ncol(mat1))
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' eset <- generate.eset(exp_mat=mat1)
#' test_transfer_table <- data.frame(
#' 'Gene'=c('Gene1','Gene1','Gene2','Gene3','Gene4'),
#' 'Transcript'=c('T11','T12','T2','T3','T3'))
#' new_eset <- update_eset.feature(use_eset=eset,
#' use_feature_info=test_transfer_table,
#' from_feature='Gene',
#' to_feature='Transcript',
#' merge_method='median',
#' distribute_method='equal'
#' )
#' print(Biobase::exprs(eset)[test_transfer_table$Gene,])
#' print(Biobase::exprs(new_eset))
#'
#' @export update_eset.feature
update_eset.feature <- function(use_eset=NULL,use_feature_info=NULL,from_feature=NULL,to_feature=NULL,
merge_method='median',distribute_method='equal'){
#
all_input_para <- c('use_eset','from_feature','to_feature')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('merge_method',c("median","mean","max","min"),envir=environment()),
check_option('distribute_method',c('mean','equal'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(use_feature_info)) use_feature_info <- Biobase::fData(use_eset)
n1 <- colnames(use_feature_info)
if(!from_feature %in% n1){
message(sprintf('%s not in in the colnames of use_feature_info, please re-try!',from_feature));return(use_eset)
}
if(!to_feature %in% n1){
message(sprintf('%s not in in the colnames of use_feature_info, please re-try!',to_feature));return(use_eset)
}
mat <- Biobase::exprs(use_eset)
use_feature_info <- base::unique(use_feature_info);
w1 <- which(use_feature_info[,1]!="" & use_feature_info[,2]!="" & is.na(use_feature_info[,1])==FALSE & is.na(use_feature_info[,2])==FALSE)
use_feature_info <- use_feature_info[w1,]
g1 <- rownames(mat) ## rownames for the expmat
f1 <- clean_charVector(use_feature_info[,from_feature]) ## from feature info
t1 <- clean_charVector(use_feature_info[,to_feature]) ## to feature info
w1 <- which(f1 %in% g1); f1 <- f1[w1]; t1 <- t1[w1]; ## only consider features in the rownames of expmat
if(base::length(w1)==0){
message(sprintf('Rownames of the expression matrix was not included in the %s column, please check and re-try !',from_feature))
return(use_eset)
}
message(sprintf('%d transfer pairs related with %d rows from original expression matrix will be keeped !',base::length(w1),base::length(g1)))
fc1 <- base::table(f1); tc1 <- base::table(t1); fw1 <- which(fc1>1); tw1 <- which(tc1>1); ## check duplicate records
if(base::length(fw1)>0){
message(sprintf('Original feature %s has %d items with duplicate records, will distribute the original values equal to all related items !
if do not want this, please check and retry !',from_feature,base::length(fw1)))
#return(use_eset)
w2 <- which(f1 %in% names(fw1)) ## need to distribute
w0 <- base::setdiff(1:base::length(f1),w2) ## do not need to distribute
if(distribute_method=='equal'){
v1 <- mat[f1[w2],]; ## distribute equal
}
if(distribute_method=='mean'){
v1 <- mat[f1[w2],]; ## distribute mean
tt <- as.numeric(base::table(f1[w2])[f1[w2]])
v1 <- v1/tt;
}
rownames(v1) <- paste0(f1[w2],'-',t1[w2]);
f1[w2] <- paste0(f1[w2],'-',t1[w2]); # update transfer table
mat <- base::rbind(v1,mat[f1[w0],]) # update mat table
fc1 <- base::table(f1); tc1 <- base::table(t1); fw1 <- which(fc1>1); tw1 <- which(tc1>1); ## update f1, t1 and related values
}
if(base::length(tw1)>0){
w2 <- which(t1 %in% names(tw1)) ## need to merge
w0 <- base::setdiff(1:base::length(t1),w2) ## do not need to merge
mat_new_0 <- mat[f1[w0],]; rownames(mat_new_0) <- t1[w0] ## mat do not need to merge
if(merge_method=='mean') tmp1 <- stats::aggregate(mat[f1[w2],,drop=FALSE],list(t1[w2]),function(x){base::mean(x,na.rm=TRUE)})
if(merge_method=='median') tmp1 <- stats::aggregate(mat[f1[w2],,drop=FALSE],list(t1[w2]),function(x){stats::median(x,na.rm=TRUE)})
if(merge_method=='max') tmp1 <- stats::aggregate(mat[f1[w2],,drop=FALSE],list(t1[w2]),function(x){base::max(x,na.rm=TRUE)})
if(merge_method=='min') tmp1 <- stats::aggregate(mat[f1[w2],,drop=FALSE],list(t1[w2]),function(x){base::min(x,na.rm=TRUE)})
mat_new_1 <- tmp1[,-1]; rownames(mat_new_1) <- tmp1[,1] ## mat merged
mat_new <- base::rbind(mat_new_0,mat_new_1)
}else{
mat_new <- mat[f1,]
rownames(mat_new) <- t1
}
new_eset <- generate.eset(exp_mat=mat_new, phenotype_info=Biobase::pData(use_eset), feature_info=NULL, annotation_info=annotation(use_eset))
return(new_eset)
}
#' Reassign the phenoData slot of ExpressionSet and Update phenotype information
#'
#' \code{update_eset.phenotype} reassigns the phenoData slot of ExpressionSet based on user's demand.
#' It is mainly used to modify sample names and extract interested phenotype information for further sample clustering.
#'
#' @param use_eset ExpressionSet class object.
#' @param use_phenotype_info data.frame, a dataframe contains phenotype information, can be obtained by calling \code{pData} function.
#' @param use_sample_col character, must be one of the column names in \code{use_phenotype_info}.
#' @param use_col character, the columns will be kept from \code{use_phenotype_info}.
#' 'auto', only extracting 'cluster-meaningful' sample features (e.g. it is meaningless to use 'gender' as clustering feature, if all samples are female).
#' 'GEO-auto' means it will extract the following selected columns,
#' "geo_accession", "title", "source_name_ch1", and columns ended with ":ch1". Default is "auto".
#' @return Return an ExressionSet object with updated phenotype information.
#' @examples
#' \dontrun{
#' net_eset <- load.exp.GEO(out.dir='./test',
#' GSE='GSE116028',
#' GPL='GPL6480',
#' getGPL=TRUE,
#' update=FALSE)
#' net_eset <- update_eset.phenotype(use_eset=net_eset,
#' use_phenotype_info=Biobase::pData(net_eset),
#' use_sample_col='geo_accession',
#' use_col='GEO-auto')
#' }
#' @export update_eset.phenotype
update_eset.phenotype <- function(use_eset=NULL,use_phenotype_info=NULL,use_sample_col=NULL,use_col='auto'){
#
all_input_para <- c('use_eset')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(use_eset)){
message('use_eset required, please re-try !');
return(use_eset)
}
if(is.null(use_phenotype_info)) use_phenotype_info <- Biobase::pData(use_eset)
if(is.null(use_sample_col)==FALSE){
if(!use_sample_col %in% colnames(use_phenotype_info)){
stop(sprintf('%s not in the colnames of use_phenotype_info, please re-try!',use_sample_col));#return(use_eset)
}
}
if(is.null(use_col)) use_col <- colnames(use_phenotype_info)
mat <- Biobase::exprs(use_eset)
s1 <- colnames(mat) ## all samples
if(is.null(use_sample_col)==TRUE){
p1 <- rownames(use_phenotype_info)
}else{
p1 <- use_phenotype_info[,use_sample_col] ## sample in the phenotype info
}
w1 <- which(p1 %in% s1); p1 <- p1[w1]; ## only consider samples in the colnames of expmat
if(base::length(w1)==0){
if(is.null(use_sample_col)==TRUE){
message('Colnames of the expression matrix was not included in the rownames of use_phenotype_info, please check and re-try !')
}else{
message(sprintf('Colnames of the expression matrix was not included in the %s column, please check and re-try !',use_sample_col))
}
return(use_eset)
}
message(sprintf('%d out of %d samples from the expression matrix will be keeped !',base::length(w1),base::length(s1)))
mat_new <- mat[,p1]
use_phenotype_info <- use_phenotype_info[w1,]
n1 <- colnames(use_phenotype_info)
if(use_col[1] == 'GEO-auto'){
w1 <- c('geo_accession','title','source_name_ch1',n1[grep(':ch1',n1)])
p1 <- use_phenotype_info[,w1]
colnames(p1)[4:ncol(p1)] <- gsub('(.*):ch1','\\1',colnames(p1)[4:ncol(p1)])
colnames(p1)[3] <- gsub('(.*)_ch1','\\1',colnames(p1)[3])
if(is.null(use_sample_col)==FALSE) rownames(p1) <- use_phenotype_info[,use_sample_col]
if(base::length(w1)>1) p1 <- as.data.frame(apply(p1,2,clean_charVector),stringsAsFactors=FALSE)
if(base::length(w1)==1) p1 <- as.data.frame(clean_charVector(p1),stringsAsFactors=FALSE)
new_phenotype_info <- p1;
}else{
if(use_col[1] == 'auto'){
u1 <- apply(use_phenotype_info,2,function(x)base::length(base::unique(x)))
w1 <- which(u1>=2 & u1<=nrow(use_phenotype_info)-1)
if(base::length(w1)==0){
message('No column could match the auto criteria, please check and re-try!');return(FALSE)
}
p1 <- use_phenotype_info[,w1]
if(base::length(w1)>1) p1 <- as.data.frame(apply(p1,2,clean_charVector),stringsAsFactors=FALSE)
if(base::length(w1)==1) p1 <- as.data.frame(clean_charVector(p1),stringsAsFactors=FALSE)
new_phenotype_info <- use_phenotype_info[,w1];names(new_phenotype_info) <- names(use_phenotype_info)[w1];
}else{
if(base::length(base::setdiff(use_col,n1))>0){
message(sprintf('%s not in use_phenotype_info, please re-try!',base::paste(base::setdiff(use_col,n1),collapse=';')));return(FALSE)
}
p1 <- use_phenotype_info[,use_col]
if(base::length(use_col)>1) p1 <- as.data.frame(apply(p1,2,clean_charVector),stringsAsFactors=FALSE)
if(base::length(use_col)==1) p1 <- as.data.frame(clean_charVector(p1),stringsAsFactors=FALSE)
new_phenotype_info <- p1; names(new_phenotype_info) <- use_col;
}
}
rownames(new_phenotype_info) <- rownames(use_phenotype_info)
#print(new_phenotype_info)
message(sprintf('%d out of %d sample features will be keeped !',ncol(new_phenotype_info),ncol(use_phenotype_info)))
new_eset <- generate.eset(exp_mat=mat_new, phenotype_info=new_phenotype_info, feature_info=Biobase::fData(use_eset), annotation_info=annotation(use_eset))
return(new_eset)
}
#' IQR (interquartile range) filter to extract genes from expression matrix
#'
#' \code{IQR.filter} is a function to extract genes from the expression matrix by setting threshold to their IQR value.
#' IQR (interquartile range) is a measure of statistical dispersion. It is calculated for each gene across all the samples.
#' By setting threshold value, genes with certain statistical dispersion across samples will be filtered out.
#' This step is mainly used to perform sample cluster and to prepare the input for SJAracne.
#'
#' @param exp_mat matrix, the gene expression matrix. Each row represents a gene/transcript/probe, each column represents a sample.
#' @param use_genes a vector of characters, the gene list needed to be filtered. Default is the row names of \code{exp_mat}.
#' @param thre numeric, the threshold for IQR of the genes in \code{use_genes}. Default is 0.5.
#' @param loose_gene a vector of characters, the gene list that only need to pass the \code{loose_thre}.
#' This parameter is designed for the input of possible drivers used in SJAracne. Default is NULL.
#' @param loose_thre numeric, the threshold for IQR of the genes in \code{loose_gene}. Default is 0.1.
#' @return Return a vector with logical values indicate which genes should be kept.
#' @examples
#' mat1 <- matrix(rnorm(15000),nrow=1500,ncol=10)
#' colnames(mat1) <- paste0('Sample',1:ncol(mat1))
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' choose1 <- IQR.filter(mat1,thre=0.5,
#' loose_gene=paste0('Gene',1:100))
#' @export
IQR.filter <- function(exp_mat,use_genes=rownames(exp_mat),thre = 0.5,loose_gene=NULL,loose_thre=0.1) {
#
all_input_para <- c('exp_mat','use_genes','thre')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
use_genes <- base::intersect(use_genes,rownames(exp_mat))
use_genes <- base::setdiff(use_genes,"")
exp_mat <- exp_mat[use_genes,]
iqr <- apply(exp_mat, 1, stats::IQR) ## calculate IQR for each gene
choose0 <- use_genes[iqr > quantile(iqr, loose_thre)] ## for loose_gene
choose1 <- use_genes[iqr > quantile(iqr, thre)] ## for all genes
choose2 <- base::unique(c(base::intersect(loose_gene, choose0), choose1)) ## union set
use_vec <- rep(FALSE,length.out=base::length(use_genes));names(use_vec) <- use_genes
use_vec[choose2] <- TRUE
print(base::table(use_vec))
return(use_vec)
}
#' Normalization of RNA-Seq Reads Count
#'
#' \code{RNASeqCount.normalize.scale} is a simple version to normalize the RNASeq reads count data.
#'
#' Users can also load \code{load.exp.RNASeq.demo}, and follow the \code{DESeq2} pipeline for RNASeq data processing.
#' Warning, \code{load.exp.RNASeq.demo} and \code{load.exp.RNASeq.demoSalmon} in NetBID2 may not cover all the possible scenarios.
#'
#' @param mat matrix, matrix of RNA-Seq reads data. Each row is a gene/transcript, each column is a sample.
#' @param total integer, total RNA-Seq reads count. If NULL, will use the mean of each column's summation. Default is NULL.
#' @param pseudoCount integer, the integer added to avoid "-Inf" showing up during log transformation. Default is 1.
#'
#' @return Return a numeric matrix, containing the normalized RNA-Seq reads count.
#'
#' @examples
#' mat1 <- matrix(rnbinom(10000, mu = 10, size = 1),nrow=1000,ncol=10)
#' colnames(mat1) <- paste0('Sample1',1:ncol(mat1))
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' norm_mat1 <- RNASeqCount.normalize.scale(mat1)
#' @export
RNASeqCount.normalize.scale <- function(mat,
total = NULL,
pseudoCount = 1) {
#
all_input_para <- c('mat','pseudoCount')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
d <- mat
if (!is.data.frame(d))
d <- data.frame(d)
if (!all(d > 0))
d <- d + pseudoCount
s <- apply(d, 2, sum)
m <-
ifelse(is.null(total), as.integer(base::mean(s)), as.integer(total)) ## total or mean sum
options(digits = 2 + nchar(m))
fac <- m / s
for (i in 1:base::length(s)) {
d[, i] <- d[, i] * fac[i]
#d[, i] <- round(d[, i] * fac[i], 0)
}
if (!all(d > 0))
d <- d + pseudoCount
d
}
## inner function for dist2
dist2.mod <- function (x, fun = function(a, b) base::mean(abs(a - b), na.rm = TRUE),
diagonal = 0)
{
if (!(is.numeric(diagonal) && (base::length(diagonal) == 1)))
stop("'diagonal' must be a numeric scalar.")
if (missing(fun)) {
res = apply(x, 2, function(w) base::colMeans(abs(x - w), na.rm = TRUE))
}
else {
res = matrix(diagonal, ncol = ncol(x), nrow = ncol(x))
if (ncol(x) >= 2) {
for (j in 2:ncol(x)) for (i in 1:(j - 1)) res[i,
j] = res[j, i] = fun(x[, i], x[, j])
}
}
colnames(res) = rownames(res) = colnames(x)
return(res)
}
########################### activity-related functions
## functions for activity score calculation, mean, absmean, maxmean, weighted mean ?
es <- function(z, es.method = "mean") {
if (es.method == "maxmean") {
n <- base::length(z)
m1 <- ifelse(sum(z > 0) > 0, sum(z[z > 0]) / n, 0)
m2 <- ifelse(sum(z < 0) > 0, sum(z[z < 0]) / n, 0)
if (m1 > -m2)
es <- m1
else
es <- m2
}
else if (es.method == 'absmean') {
es <- base::mean(abs(z),na.rm=TRUE)
}
else if (es.method == 'mean') {
es <- base::mean(z,na.rm=TRUE)
}
else if (es.method == 'median') {
es <- stats::median(z,na.rm=TRUE)
}
else if (es.method == 'max') {
es <- base::max(z,na.rm=TRUE)
}
else if (es.method == 'min') {
es <- base::min(z,na.rm=TRUE)
}
return(es)
}
do.std <- function(x) {
x <- x[!is.na(x)]
(x - base::mean(x,na.rm=TRUE)) / sd(x,na.rm=TRUE)
}
#' Calculate Activity Value for Each Driver
#'
#' \code{cal.Activity} calculates the activity value for each driver.
#' This function requires two inputs, the driver-to-target list object \code{target_list} and the expression matrix.
#'
#' @param target_list list, the driver-to-target list object. Either igraph_obj or target_list is necessary for this function.
#' The names of the list elements are drivers.
#' Each element is a data frame, usually contains at least three columns.
#' "target", target gene names;
#' "MI", mutual information;
#' "spearman", spearman correlation coefficient.
#' "MI" and "spearman" is necessary if es.method="weightedmean".
#' Users can call \code{get.SJAracne.network} to get this list from the network file generated by SJAracne (the second element of the return list)
#' or prepare the list object by hand but should match the data format described above.
#' @param igraph_obj igraph object, optional. Either igraph_obj or target_list is necessary for this function.
#' Users can call \code{get.SJAracne.network} to get this object from the network file generated by SJAracne (the third element of the return list),
#' or prepare the igraph network object by hand (Directed network and the edge attributes should include "weight" and "sign" if es.method="weightedmean").
#' @param cal_mat numeric matrix, the expression matrix of genes/transcripts.
#' @param es.method character, method applied to calculate the activity value. User can choose from "mean", "weightedmean", "maxmean" and "absmean".
#' Default is "weightedmean".
#' @param std logical, if TRUE, the expression matrix will be normalized by column. Default is TRUE.
#' @param memory_constrain logical, if TRUE, the calculation strategy will not use Matrix Cross Products, which is memory consuming.
#' Default is FALSE.
#' @return Return a matrix of activity values. Rows are drivers, columns are samples.
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' ac_mat <- cal.Activity(target_list=analysis.par$merge.network$target_list,
#' cal_mat=Biobase::exprs(analysis.par$cal.eset),
#' es.method='weightedmean')
#' ac_mat <- cal.Activity(igraph_obj=analysis.par$merge.network$igraph_obj,
#' cal_mat=Biobase::exprs(analysis.par$cal.eset),
#' es.method='maxmean')
#' @export
cal.Activity <- function(target_list=NULL, igraph_obj = NULL, cal_mat=NULL, es.method = 'weightedmean',std=TRUE,memory_constrain=FALSE) {
#
all_input_para <- c('cal_mat','es.method','std','memory_constrain')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('memory_constrain',c(TRUE,FALSE),envir=environment()),
check_option('std',c(TRUE,FALSE),envir=environment()),
check_option('es.method',c('mean','weightedmean','maxmean','absmean'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(target_list)==TRUE & is.null(igraph_obj)==TRUE){
message('Either target_list or igraph_obj is required, please check and re-try!');return(FALSE);
}
if(is.null(target_list)==FALSE & memory_constrain==TRUE){
ac.mat <- cal.Activity.old(target_list=target_list, cal_mat=cal_mat, es.method = es.method,std=std)
return(ac.mat)
}
if(is.null(target_list)==TRUE & memory_constrain==TRUE){
message('Only accepts target_list input when memory_constrain=TRUE, please check and re-try!');return(FALSE);
}
if(nrow(cal_mat)==0){
message('No genes in the cal_mat, please check and re-try!');return(FALSE);
}
if(std==TRUE) cal_mat <- apply(cal_mat, 2, do.std)
if(is.null(igraph_obj)==FALSE){
gr <- igraph_obj
mat1 <- get_igraph2matrix(gr,es.method=es.method)
mat2 <- get_igraph2matrix(gr,es.method='mean')
all_source <- get_gr2driver(gr)
}else{
if(is.null(target_list)==FALSE){
mat1 <- get_target_list2matrix(target_list,es.method=es.method)
mat2 <- get_target_list2matrix(target_list,es.method='mean')
all_source <- names(target_list)
}
}
##
mat1_source <- mat1[all_source,,drop=FALSE]
w1 <- base::intersect(rownames(cal_mat),colnames(mat1_source))
if(base::length(w1)==0){
message('No intersected genes found for the cal_mat and target in the network, please check and re-try!');
return(FALSE)
}
use_mat1_source <- mat1_source[,w1,drop=FALSE] ## network info
use_mat2_source <- mat2[all_source,w1,drop=FALSE] ## network binary info
## weighted mean + mean
if(es.method %in% c('weightedmean','mean')){
use_cal_mat <- cal_mat[w1,,drop=FALSE] ## expression info
out_mat <- use_mat1_source %*% use_cal_mat
out_mat <- out_mat/Matrix::rowSums(use_mat2_source) ## get mean
}
## absmean
if(es.method == 'absmean'){
use_cal_mat <- cal_mat[w1,,drop=FALSE] ## expression info
out_mat <- use_mat1_source %*% abs(use_cal_mat)
out_mat <- out_mat/Matrix::rowSums(use_mat2_source) ## get mean
}
## maxmean
if(es.method == 'maxmean'){
use_cal_mat <- cal_mat[w1,,drop=FALSE] ## expression info
use_cal_mat_pos <- use_cal_mat;use_cal_mat_pos[which(use_cal_mat_pos<0)] <- 0;
use_cal_mat_neg <- use_cal_mat;use_cal_mat_neg[which(use_cal_mat_neg>0)] <- 0;
out_mat_pos <- use_mat1_source %*% use_cal_mat_pos
out_mat_pos <- out_mat_pos/Matrix::rowSums(use_mat2_source) ## get mean
out_mat_neg <- use_mat1_source %*% use_cal_mat_neg
out_mat_neg <- out_mat_neg/Matrix::rowSums(use_mat2_source) ## get mean
out_mat_sign <- sign(abs(out_mat_pos)-abs(out_mat_neg))
out_mat_sign_pos <- out_mat_sign; out_mat_sign_pos[out_mat_sign_pos!=1] <-0;
out_mat_sign_neg <- out_mat_sign; out_mat_sign_neg[out_mat_sign_neg!= -1] <-0;
out_mat <- out_mat_pos*out_mat_sign_pos-out_mat_neg*out_mat_sign_neg
}
## median, min, max , not supported
# output
ac.mat <- as.matrix(out_mat)
w1 <- which(is.na(ac.mat[,1])==FALSE)
if(base::length(w1)==0){
message('Fail in calculating activity, please check the ID type in cal_mat and target_list and try again !')
}
ac.mat <- ac.mat[w1,,drop=FALSE]
return(ac.mat)
}
## inner functions
cal.Activity.old <- function(target_list=NULL, cal_mat=NULL, es.method = 'weightedmean',std=TRUE) {
## mean, absmean, maxmean, weightedmean
use_genes <- row.names(cal_mat)
if(base::length(use_genes)==0){
message('No genes in the cal_mat, please check and re-try!');return(FALSE);
}
all_target <- target_list
#all_target <- all_target[base::intersect(use_genes, names(all_target))] ## if the driver is not included in cal_mat but its target genes are included, will also calculate activity
ac.mat <-
matrix(NA, ncol = ncol(cal_mat), nrow = base::length(all_target)) ## generate activity matrix, each col for sample, each row for source target
#z-normalize each sample
if(std==TRUE) cal_mat <- apply(cal_mat, 2, do.std)
for (i in 1:base::length(all_target)) {
x <- names(all_target)[i]
x1 <- all_target[[x]]
x2 <- base::unique(base::intersect(rownames(x1), use_genes)) ## filter target by cal genes
x1 <- x1[x2, ] ## target info
target_num <- base::length(x2)
if (target_num == 0)
next
if (target_num == 1){
if (es.method != 'weightedmean') ac.mat[i, ] <- cal_mat[x2,] # 20230228
if (es.method == 'weightedmean') ac.mat[i, ] <- cal_mat[x2,]*x1$MI * sign(x1$spearman) # 20230228
next
}
if (es.method != 'weightedmean')
ac.mat[i, ] <- apply(cal_mat[x2,,drop=FALSE], 2, es, es.method)
if (es.method == 'weightedmean') {
weight <- x1$MI * sign(x1$spearman)
ac.mat[i, ] <- apply(cal_mat[x2,,drop=FALSE] * weight, 2, es, 'mean')
}
}
rownames(ac.mat) <- names(all_target)
colnames(ac.mat) <- colnames(cal_mat)
w1 <- which(is.na(ac.mat[,1])==FALSE)
if(base::length(w1)==0){
message('Fail in calculating activity, please check the ID type in cal_mat and target_list and try again !')
}
ac.mat <- ac.mat[w1,]
return(ac.mat)
}
get_target_list2matrix <- function(target_list=NULL,es.method = 'weightedmean') {
all_source <- names(target_list)
all_target <- base::unique(unlist(lapply(target_list,function(x)rownames(x))))
mat1 <- matrix(0,nrow=base::length(all_source),ncol=base::length(all_target))
rownames(mat1) <- all_source;
colnames(mat1) <- all_target;
for(i in all_source){
if(es.method!='weightedmean') mat1[i,rownames(target_list[[i]])] <- 1;
if(es.method=='weightedmean') mat1[i,rownames(target_list[[i]])] <- target_list[[i]]$MI*sign(target_list[[i]]$spearman);
}
return(mat1)
}
get_igraph2matrix <- function(gr=NULL,es.method = 'weightedmean'){
if(es.method=='weightedmean'){
if('weight' %in% igraph::edge_attr_names(gr) & 'sign' %in% igraph::edge_attr_names(gr)){
mat1 <- igraph::as_adjacency_matrix(gr,type='both',attr = 'weight')
mat2 <- igraph::as_adjacency_matrix(gr,type='both',attr = 'sign')
mat1 <- mat1*mat2
}else{
message('weight, sign attributes are not included in the input igraph object, weightedmean can not be used !')
return(FALSE)
}
}else{
mat1 <- igraph::as_adjacency_matrix(gr,type='both')
}
return(mat1)
}
get_gr2driver <- function(gr,mode='out'){
d1 <- igraph::degree(gr,mode=mode)
names(d1[which(d1>0)])
}
#' Clean Activity-based profile
#'
#' \code{processDriverProfile} is a helper function to pre-process Activity-based profile.
#'
#' @param Driver_profile a numeric vector, contain statistics for drivers
#' (e.g driver's target size, driver's Z-statistics)
#' @param Driver_name a character vector, contain name for drivers.
#' The length of `Driver_profile` and `Driver_name` must be equal and the order of item must match.
#' @param choose_strategy character, strategy of selection if duplicate driver name (e.g TP53_TF, TP53_SIG).
#' Choose from "min","max","absmin","absmax". Default is "min".
#' @param return_type character, strategy of return type.
#' Choose from "driver_name","gene_name", "driver_statistics", "gene_statistics".
#' If choose "*_name", only the name vector is returned.
#' If choose "*_statistics", the statistics vector is returned with character name.
#' Default is "driver_name".
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' ms_tab <- analysis.par$final_ms_tab
#' Driver_profile <- ms_tab$P.Value.G4.Vs.WNT_DA
#' Driver_name <- ms_tab$gene_label
#' res1 <- processDriverProfile(Driver_profile=Driver_profile,
#' Driver_name=Driver_name,
#' choose_strategy='min')
#' res2 <- processDriverProfile(Driver_profile=Driver_profile,
#' Driver_name=Driver_name,
#' return_type = 'gene_name',
#' choose_strategy='min')
#' Driver_profile <- ms_tab$Z.G4.Vs.WNT_DA
#' res3 <- processDriverProfile(Driver_profile=Driver_profile,
#' Driver_name=Driver_name,
#' choose_strategy='absmax',
#' return_type = 'driver_statistics')
#' res4 <- processDriverProfile(Driver_profile=Driver_profile,
#' Driver_name=Driver_name,
#' choose_strategy='absmax',
#' return_type = 'gene_statistics')
#' driver_size <- ms_tab$Size
#' res5 <- processDriverProfile(Driver_profile=driver_size,
#' Driver_name=Driver_name,
#' choose_strategy='max')
#' @export
processDriverProfile <- function(Driver_profile,Driver_name,
choose_strategy='min',
return_type='driver_name'){
all_input_para <- c('Driver_profile','Driver_name','choose_strategy')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('choose_strategy',c('min','max','absmin','absmax'),envir=environment()),
check_option('return_type',c('driver_name','gene_name','driver_statistics','gene_statistics'),envir=environment()))
if(min(check_res)==0){message('Please check and re-try!');return(FALSE)}
##
l1 <- length(Driver_profile)
l2 <- length(Driver_name)
if(l1!=l2 | l1==0 | l2==0){
message('Driver_profile and Driver_name need same length.
Please check Driver_profile, Driver_name, and re-try!');return(FALSE)
}
gene_name <- gsub('(.*)_(TF|SIG)',"\\1",Driver_name)
uni_gene_name <- unique(gene_name)
tmp1 <- lapply(uni_gene_name,function(x){
w1 <- which(gene_name == x)
x1 <- Driver_profile[w1]
x2 <- rank(x1)
x3 <- rank(abs(x1))
if(choose_strategy=='min'){w2 <- w1[which.min(x2)]}
if(choose_strategy=='max'){w2 <- w1[which.max(x2)]}
if(choose_strategy=='absmin'){w2 <- w1[which.min(x3)]}
if(choose_strategy=='absmax'){w2 <- w1[which.max(x3)]}
w2
})
remain_item <- unlist(tmp1)
if(return_type=='driver_name') new_vec <- Driver_name[remain_item]
if(return_type=='gene_name') new_vec <- gene_name[remain_item]
if(return_type=='driver_statistics'){new_vec <- Driver_profile[remain_item]; names(new_vec) <- Driver_name[remain_item]}
if(return_type=='gene_statistics'){new_vec <- Driver_profile[remain_item]; names(new_vec) <- gene_name[remain_item]}
return(new_vec)
}
#' Calculate Activity Value for Gene Sets
#'
#' \code{cal.Activity.GS} calculates activity value for each gene set, and return a numeric matrix with rows of gene sets and columns of samples.
#'
#' @param use_gs2gene list, contains elements of gene sets. Element name is gene set name, each element contains a vector of genes belong to that gene set.
#' Default is using \code{all_gs2gene[c('H','CP:BIOCARTA','CP:REACTOME','CP:KEGG_MEDICU')]}, which is loaded from \code{gs.preload}.
#' @param cal_mat numeric matrix, gene/transcript expression matrix.
#' If want to input activity matrix, need to use `processDriverProfile()` to pre-process the dataset.
#' Detailed could see demo.
#' @param es.method character, method to calculate the activity value. Users can choose from "mean", "absmean", "maxmean", "gsva", "ssgsea", "zscore" and "plage".
#' The details for using the last four options, users can check \code{gsva}. Default is "mean".
#' @param std logical, if TRUE, the expression matrix will be normalized by column. Default is TRUE.
#' @return Return an activity matrix with rows of gene sets and columns of samples.
#'
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' gs.preload(use_spe='Homo sapiens',update=FALSE)
#' use_gs2gene <- merge_gs(all_gs2gene=all_gs2gene,
#' use_gs=c('H','CP:BIOCARTA','CP:REACTOME','CP:KEGG_MEDICU','C5'))
#' exp_mat_gene <- Biobase::exprs(analysis.par$cal.eset)
#' ## each row is a gene symbol, if not, must convert ID first
#' ac_gs <- cal.Activity.GS(use_gs2gene = use_gs2gene,
#' cal_mat = exp_mat_gene)
#' ## if want to input activity-matrix
#' ac_mat <- cal.Activity(target_list=analysis.par$merge.network$target_list,
#' cal_mat=Biobase::exprs(analysis.par$cal.eset),
#' es.method='weightedmean')
#' # pre-process the activity matrix by selecting the one
#' # with larger target size for duplicate drivers
#' Driver_name <- rownames(ac_mat)
#' ms_tab <- analysis.par$final_ms_tab
#' driver_size <- ms_tab[Driver_name,]$Size
#' use_driver <- processDriverProfile(Driver_profile=driver_size,
#' Driver_name=Driver_name,
#' choose_strategy='max',
#' return_type='driver_name')
#' use_driver_gene_name <-
#' processDriverProfile(Driver_profile=driver_size,
#' Driver_name=Driver_name,
#' choose_strategy='max',
#' return_type='gene_name')
#' ac_mat_gene <- ac_mat[use_driver,]
#' rownames(ac_mat_gene) <- use_driver_gene_name
#' driver_ac_gs <- cal.Activity.GS(use_gs2gene = use_gs2gene,
#' cal_mat = ac_mat_gene)
#' @export
cal.Activity.GS <- function(use_gs2gene=all_gs2gene[c('H','CP:BIOCARTA','CP:REACTOME','CP:KEGG_MEDICU')], cal_mat=NULL, es.method = 'mean',std=TRUE) {
#
all_input_para <- c('use_gs2gene','cal_mat','es.method','std')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('std',c(TRUE,FALSE),envir=environment()),
check_option('es.method',c("mean","absmean","maxmean","gsva","ssgsea","zscore","plage"),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
while(class(use_gs2gene[[1]])=='list'){
nn <- unlist(lapply(use_gs2gene,names))
use_gs2gene <- unlist(use_gs2gene,recursive = FALSE)
names(use_gs2gene)<-nn
}
use_genes <- row.names(cal_mat)
if(es.method %in% c('gsva','ssgsea','zscore','plage')){
ac.mat <- GSVA::gsva(generate.eset(cal_mat),use_gs2gene,method=es.method)
ac.mat <- Biobase::exprs(ac.mat)
return(ac.mat)
}
ac.mat <-
matrix(NA, ncol = ncol(cal_mat), nrow = base::length(use_gs2gene)) ## generate activity matrix, each col for sample, each row for source target
#z-normalize each sample
if(std==TRUE) cal_mat <- apply(cal_mat, 2, do.std)
for (i in 1:base::length(use_gs2gene)) {
x <- names(use_gs2gene)[i]
x1 <- use_gs2gene[[x]]
x2 <- base::unique(base::intersect(x1, use_genes)) ## filter target by cal genes
target_num <- base::length(x2)
if (target_num == 0)
next
if (target_num == 1){
ac.mat[i, ] <- cal_mat[x2,]
next
}
ac.mat[i, ] <- apply(cal_mat[x2,,drop=FALSE], 2, es, es.method)
}
rownames(ac.mat) <- names(use_gs2gene)
colnames(ac.mat) <- colnames(cal_mat)
w1 <- apply(ac.mat,1,function(x)base::length(which(is.na(x)==TRUE)))
ac.mat <- ac.mat[which(w1==0),,drop=FALSE]
return(ac.mat)
}
#' Differential Expression Analysis and Differential Activity Analysis Between 2 Sample Groups Using Bayesian Inference
#'
#' \code{getDE.BID.2G} is a function performs differential gene expression analysis and differential driver activity analysis between
#' control group (parameter G0) and experimental group (parameter G1).
#'
#' @param eset ExpressionSet class object, contains gene expression data or driver activity data.
#' @param output_id_column character, the column names of Biobase::fData(eset).
#' This option is useful when the \code{eset} expression matrix is at transcript-level, and user is expecting to see the gene-level comparison statistics.
#' If NULL, rownames of the Biobase::fData(eset) will be used.
#' Default is NULL.
#' @param G1 a vector of characters, the sample names of experimental group.
#' @param G0 a vecotr of characters, the sample names of control group.
#' @param G1_name character, the name of experimental group (e.g. "Male"). Default is "G1".
#' @param G0_name character, the name of control group (e.g. "Female"). Default is "G0".
#' @param logTransformed logical, if TRUE, log tranformation of the expression value will be performed.
#' @param method character, users can choose between "MLE" and "Bayesian".
#' "MLE", the maximum likelihood estimation, will call generalized linear model(glm/glmer) to perform data regression.
#' "Bayesian", will call Bayesian generalized linear model (bayesglm) or multivariate generalized linear mixed model (MCMCglmm) to perform data regression.
#' Default is "Bayesian".
#' @param family character or family function or the result of a call to a family function.
#' This parameter is used to define the model's error distribution. See \code{?family} for details.
#' Currently, options are gaussian, poisson, binomial(for two-group sample classes)/category(for multi-group sample classes)/ordinal(for multi-group sample classes with class_ordered=TRUE).
#' If set with gaussian or poission, the response variable in the regression model will be the expression level, and the independent variable will be the sample's phenotype.
#' If set with binomial, the response variable in the regression model will be the sample phenotype, and the independent variable will be the expression level.
#' For binomial, category and ordinal input, the family will be automatically reset, based on the sample's class level and the setting of \code{class_ordered}.
#' Default is gaussian.
#' @param pooling character, users can choose from "full","no" and "partial".
#' "full", use probes as independent observations.
#' "no", use probes as independent variables in the regression model.
#' "partial", use probes as random effect in the regression model.
#' Default is "full".
#' @param verbose logical, if TRUE, sample names of both groups will be printed. Default is TRUE.
#'
#' @return
#' Return a data frame. Rows are genes/drivers, columns are "ID", "logFC", "AveExpr", "t", "P.Value", "adj.P.Val", "Z-statistics", "Ave.G1" and "Ave.G0".
#' Names of the columns may vary from different group names. Sorted by P.Value.
#'
#' @examples
#' mat <- matrix(c(0.50099,-1.2108,-1.0524,
#' 0.34881,-0.13441,0.87112,
#' 1.84579,-2.0356,-2.6025,
#' 1.62954,1.88281,1.29604),nrow=2,byrow=TRUE)
#' rownames(mat) <- c('A1','A2')
#' colnames(mat) <- c('Case-rep1','Case-rep2','Case-rep3',
#' 'Control-rep1','Control-rep2','Control-rep3')
#' tmp_eset <- generate.eset(mat,feature_info=data.frame(row.names=rownames(mat),
#' probe=rownames(mat),gene=rep('GeneX',2),
#' stringsAsFactors = FALSE))
#' res1 <- getDE.BID.2G(tmp_eset,output_id_column='probe',
#' G1=c('Case-rep1','Case-rep2','Case-rep3'),
#' G0=c('Control-rep1','Control-rep2','Control-rep3'))
#' res2 <- getDE.BID.2G(tmp_eset,output_id_column='gene',
#' G1=c('Case-rep1','Case-rep2','Case-rep3'),
#' G0=c('Control-rep1','Control-rep2','Control-rep3'))
#' res3 <- getDE.BID.2G(tmp_eset,output_id_column='gene',
#' G1=c('Case-rep1','Case-rep2','Case-rep3'),
#' G0=c('Control-rep1','Control-rep2','Control-rep3'),
#' pooling='partial')
#' \dontrun{
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' phe_info <- Biobase::pData(analysis.par$cal.eset)
#' each_subtype <- 'G4'
#' G0 <- rownames(phe_info)[which(phe_info$`subgroup`!=each_subtype)] # get sample list for G0
#' G1 <- rownames(phe_info)[which(phe_info$`subgroup`==each_subtype)] # get sample list for G1
#' DE_gene_BID <- getDE.BID.2G(eset=analysis.par$cal.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' DA_driver_BID <- getDE.BID.2G(eset=analysis.par$merge.ac.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' }
#' @export
getDE.BID.2G <-function(eset,output_id_column=NULL,G1=NULL, G0=NULL,G1_name=NULL,G0_name=NULL,method='Bayesian',family=gaussian,pooling='full',logTransformed=TRUE,verbose=TRUE){
#
all_input_para <- c('eset','G1','G0','method','family','pooling','logTransformed','verbose')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('logTransformed',c(TRUE,FALSE),envir=environment()),
check_option('verbose',c(TRUE,FALSE),envir=environment()),
check_option('method',c('Bayesian','MLE'),envir=environment()),
check_option('pooling',c('full','no','partial'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
exp_mat <- Biobase::exprs(eset)
G1 <- base::intersect(G1,colnames(exp_mat))
G0 <- base::intersect(G0,colnames(exp_mat))
if(verbose==TRUE){
print(sprintf('G1:%s', base::paste(G1, collapse = ';')))
print(sprintf('G0:%s', base::paste(G0, collapse = ';')))
}
if(base::length(G1)==0 | base::length(G0)==0){
message('Too few samples, please check the sample name of G1, G0 and samples in eset !');return(FALSE);
}
#
rn <- rownames(exp_mat)
exp_mat <- exp_mat[,c(G1,G0)]
if(is.null(dim(exp_mat))==TRUE){exp_mat <- t(as.matrix(exp_mat));rownames(exp_mat)<-rn}
if(is.null(output_id_column)==TRUE) use_id <- rownames(Biobase::fData(eset)) else use_id <- Biobase::fData(eset)[,output_id_column]
comp <- c(rep(1,length.out=base::length(G1)),rep(0,length.out=base::length(G0)))
all_id <- base::unique(use_id)
de <- lapply(all_id,function(x){
w1 <- which(use_id==x)
x1 <- exp_mat[w1,,drop=FALSE]
bid(mat=x1,use_obs_class=comp,class_order=c(0,1),family=family,method=method,
nitt=13000,burnin=3000,thin=1,pooling=pooling,class_ordered=FALSE,
logTransformed=logTransformed,std=FALSE,average.method=c('geometric'),verbose=FALSE)
})
de <- as.data.frame(do.call(base::rbind,de))
de$adj.P.Val<-p.adjust(de$P.Value,'fdr')
de$logFC<-sign(de$FC)*log2(abs(de$FC))
de$ID <- all_id
de<-de[,c('ID','logFC','AveExpr','t','P.Value','adj.P.Val','Z-statistics')]
rownames(de) <- de$ID
tT <- de;
tmp1 <- stats::aggregate(exp_mat,list(use_id),mean)
new_mat <- tmp1[,-1];rownames(new_mat) <- tmp1[,1]
tT <- tT[rownames(new_mat),,drop=FALSE]
exp_G1 <- base::rowMeans(new_mat[,G1,drop=FALSE]);
exp_G0 <- base::rowMeans(new_mat[,G0,drop=FALSE]);
tT <- base::cbind(tT,'Ave.G0'=exp_G0,'Ave.G1'=exp_G1)
if(is.null(G0_name)==FALSE) colnames(tT) <- gsub('Ave.G0',paste0('Ave.',G0_name),colnames(tT))
if(is.null(G1_name)==FALSE) colnames(tT) <- gsub('Ave.G1',paste0('Ave.',G1_name),colnames(tT))
tT <- tT[order(tT$P.Value),]
return(tT)
}
#' Combine Multiple Comparison Results from Differential Expression (DE) or Differential Activity (DA) Analysis
#'
#' \code{combineDE} combines multiple comparisons of DE or DA analysis.
#' Can combine DE with DE, DA with DA and also DE with DA if proper transfer table prepared.
#'
#' For example, there are 4 subgroups in the phenotype, G1, G2, G3 and G4. One DE analysis was performed on G1 vs. G2, and another DE was performed on G1 vs. G3.
#' If user is interested in the DE analysis between G1 vs. (G2 and G3), he can call this function to combine the two comparison results above toghether.
#' The combined P values will be taken care by \code{combinePvalVector}.
#'
#' @param DE_list list, each element in the list is one DE/DA comparison need to be combined.
#' @param DE_name a vector of characters, the DE/DA comparison names.
#' If not NULL, it must match the names of DE_list in correct order.
#' If NULL, names of the DE_list will be used.
#' Default is NULL.
#' @param transfer_tab data.frame, the ID conversion table. Users can call \code{get_IDtransfer} to get this table.
#' The purpose is to correctly mapping ID for \code{DE_list}. The column names must match \code{DE_name}.
#' If NULL, ID column of each DE comparison will be considered as the same type.
#' Default is NULL.
#' @param main_id character, a name of the element in \code{DE_list}. The ID column of that comparison will be used as the ID of the final combination.
#' If NULL, the first element name from \code{DE_list} will be used. Default is NULL.
#' @param method character, users can choose between "Stouffer" and "Fisher". Default is "Stouffer".
#' @param twosided logical, if TRUE, a two-tailed test will be performed.
#' If FALSE, a one-tailed test will be performed, and P value falls within the range of 0 to 0.5. Default is TRUE.
#' @param signed logical, if TRUE, give a sign to the P value, which indicating the direction of testing.
#' Default is TRUE.
#' @return Return a list contains the combined DE/DA analysis. Each single comparison result before combination is wrapped inside
#' (may have with some IDs filtered out, due to the combination). A data frame named "combine" inside the list is the combined analysis.
#' Rows are genes/drivers, columns are combined statistics (e.g. "logFC", "AveExpr", "t", "P.Value" etc.).
#'
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' phe_info <- Biobase::pData(analysis.par$cal.eset)
#' each_subtype <- 'G4'
#' G0 <- rownames(phe_info)[which(phe_info$`subgroup`!=each_subtype)] # get sample list for G0
#' G1 <- rownames(phe_info)[which(phe_info$`subgroup`==each_subtype)] # get sample list for G1
#' DE_gene_limma <- getDE.limma.2G(eset=analysis.par$cal.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' DA_driver_limma <- getDE.limma.2G(eset=analysis.par$merge.ac.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' DE_list <- list(DE=DE_gene_limma,DA=DA_driver_limma)
#' g1 <- gsub('(.*)_.*','\\1',DE_list$DA$ID)
#' transfer_tab <- data.frame(DE=g1,DA=DE_list$DA$ID,stringsAsFactors = FALSE)
#' res1 <- combineDE(DE_list,transfer_tab=transfer_tab,main_id='DA')
#'
#' \dontrun{
#' each_subtype <- 'G4'
#' G0 <- rownames(phe_info)[which(phe_info$`subgroup`!=each_subtype)] # get sample list for G0
#' G1 <- rownames(phe_info)[which(phe_info$`subgroup`==each_subtype)] # get sample list for G1
#' DE_gene_limma_G4 <- getDE.limma.2G(eset=analysis.par$cal.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' each_subtype <- 'SHH'
#' G0 <- rownames(phe_info)[which(phe_info$`subgroup`!=each_subtype)] # get sample list for G0
#' G1 <- rownames(phe_info)[which(phe_info$`subgroup`==each_subtype)] # get sample list for G1
#' DE_gene_limma_SHH <- getDE.limma.2G(eset=analysis.par$cal.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' DE_list <- list(G4=DE_gene_limma_G4,SHH=DE_gene_limma_SHH)
#' res2 <- combineDE(DE_list,transfer_tab=NULL)
#' }
#' @export
combineDE<-function(DE_list,DE_name=NULL,transfer_tab=NULL,main_id=NULL,method='Stouffer',twosided=TRUE,signed=TRUE){
#
all_input_para <- c('DE_list','method','twosided','signed')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('twosided',c(TRUE,FALSE),envir=environment()),
check_option('signed',c(TRUE,FALSE),envir=environment()),
check_option('method',c('Stouffer','Fisher'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
nDE<-base::length(DE_list)
if(nDE<2) stop('At least two DE outputs are required for combineDE analysis!\n')
if(is.null(DE_name)==TRUE){
DE_name <- names(DE_list)
}
if(is.null(transfer_tab)==TRUE){
transfer_tab <- lapply(DE_list,function(x)x$ID)
names(transfer_tab) <- DE_name
w1 <- names(which(base::table(unlist(transfer_tab))==nDE))
if(base::length(w1)==0){message('No intersected IDs found between DE list, please check and re-try!');return(FALSE);}
transfer_tab <- do.call(base::cbind,lapply(transfer_tab,function(x)w1))
transfer_tab <- as.data.frame(transfer_tab,stringsAsFactors=FALSE)
}
w1 <- lapply(DE_name,function(x){
which(transfer_tab[[x]] %in% DE_list[[x]]$ID)
})
w11 <- which(base::table(unlist(w1))==nDE)
if(base::length(w11)==0){
message('No intersected IDs found between DE list, please check and re-try!');return(FALSE);
}
w2 <- as.numeric(names(w11))
transfer_tab <- transfer_tab[w2,]
combine_info <- lapply(DE_name,function(x1){
DE_list[[x1]][transfer_tab[[x1]],]
})
names(combine_info) <- DE_name
dd <- do.call(base::cbind,lapply(combine_info,function(x1){
x1$P.Value*sign(x1$logFC)
}))
res1 <- t(apply(dd,1,function(x){
combinePvalVector(x,method=method,signed=signed,twosided=twosided)
}))
res1 <- as.data.frame(res1)
if(is.null(main_id)==TRUE) main_id <- DE_name[1]
res1$adj.P.Val <- p.adjust(res1$P.Value,'fdr')
res1$logFC <- base::rowMeans(do.call(base::cbind,lapply(combine_info,function(x)x$logFC)))
res1$AveExpr <- base::rowMeans(do.call(base::cbind,lapply(combine_info,function(x)x$AveExpr)))
res1 <- base::cbind(ID=transfer_tab[,main_id],transfer_tab,res1,stringsAsFactors=FALSE)
rownames(res1) <- res1$ID
combine_info$combine <- res1
return(combine_info)
}
# inner function: class_label can be obtained by get_class
get_class2design <- function(class_label){
design <- model.matrix(~0+class_label);colnames(design) <- base::unique(class_label);
rownames(design) <- names(class_label)
return(design)
#design.mat <-as.data.frame(matrix(0, nrow = base::length(class_label), ncol = base::length(base::unique(class_label))))
#rownames(design.mat) <- names(class_label) ## sample
#colnames(design.mat) <- base::unique(class_label)
#for(i in 1:base::length(class_label)){design.mat[names(class_label)[i],class_label[i]]<-1;}
#return(design.mat)
}
#' Differential Expression Analysis and Differential Activity Analysis Between 2 Sample Groups Using Limma
#'
#' \code{getDE.limma.2G} is a function performs differential gene expression analysis and differential driver activity analysis
#' between control group (parameter G0) and experimental group (parameter G1), using limma related functions.
#'
#' @param eset ExpressionSet class object, contains gene expression data or driver activity data.
#' @param G1 a vector of characters, the sample names of experimental group.
#' @param G0 a vecotr of characters, the sample names of control group.
#' @param G1_name character, the name of experimental group (e.g. "Male"). Default is "G1".
#' @param G0_name character, the name of control group (e.g. "Female"). Default is "G0".
#' @param verbose logical, if TRUE, sample names of both groups will be printed. Default is TRUE.
#' @param random_effect a vector of characters, vector or factor specifying a blocking variable.
#' Default is NULL, no random effect will be considered.
#'
#' @return
#' Return a data frame. Rows are genes/drivers, columns are "ID", "logFC", "AveExpr", "t", "P.Value", "adj.P.Val", "B", "Z-statistics", "Ave.G1" and "Ave.G0".
#' Names of the columns may vary from different group names. Sorted by P-values.
#'
#'
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' phe_info <- Biobase::pData(analysis.par$cal.eset)
#' each_subtype <- 'G4'
#' G0 <- rownames(phe_info)[which(phe_info$`subgroup`!=each_subtype)] # get sample list for G0
#' G1 <- rownames(phe_info)[which(phe_info$`subgroup`==each_subtype)] # get sample list for G1
#' DE_gene_limma <- getDE.limma.2G(eset=analysis.par$cal.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' DA_driver_limma <- getDE.limma.2G(eset=analysis.par$merge.ac.eset,
#' G1=G1,G0=G0,
#' G1_name=each_subtype,
#' G0_name='other')
#' @export
getDE.limma.2G <- function(eset=NULL, G1=NULL, G0=NULL,G1_name=NULL,G0_name=NULL,verbose=TRUE,random_effect=NULL) {
#
all_input_para <- c('eset','G1','G0','verbose')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('verbose',c(TRUE,FALSE),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
exp_mat <- Biobase::exprs(eset)
G1 <- base::intersect(G1,colnames(exp_mat))
G0 <- base::intersect(G0,colnames(exp_mat))
if(verbose==TRUE){
print(sprintf('G1:%s', base::paste(G1, collapse = ';')))
print(sprintf('G0:%s', base::paste(G0, collapse = ';')))
}
if(base::length(G1)==0 | base::length(G0)==0){
message('Too few samples, please check the sample name of G1, G0 and samples in eset !');return(FALSE);
}
#
all_samples <- colnames(Biobase::exprs(eset))
use_samples <- c(G0, G1)
phe <- as.data.frame(Biobase::pData(eset)[use_samples, ,drop=FALSE]);
rownames(phe) <- use_samples
new_eset <- generate.eset(exp_mat=Biobase::exprs(eset)[, use_samples,drop=F],phenotype_info=phe, feature_info=Biobase::fData(eset))
new_mat <- Biobase::exprs(new_eset)
##
design.mat <-as.data.frame(matrix(NA, nrow = base::length(use_samples), ncol = 1))
rownames(design.mat) <-use_samples
colnames(design.mat) <- 'group'
design.mat[base::intersect(G0, use_samples), 'group'] <- 'G0'
design.mat[base::intersect(G1, use_samples), 'group'] <- 'G1'
# design <- model.matrix( ~ group + 0, design.mat)
group <- factor(design.mat$group)
design <- model.matrix(~0+group);
colnames(design) <- levels(group); rownames(design) <- colnames(new_mat)
if(is.null(random_effect)==TRUE){
fit <- limma::lmFit(new_mat,design)
}else{
random_effect <- random_effect[colnames(Biobase::exprs(new_eset))]
corfit <- limma::duplicateCorrelation(new_eset,design,block=random_effect)
fit <- limma::lmFit(new_mat,design,block=random_effect,correlation=corfit$consensus)
}
contrasts <- limma::makeContrasts(G1-G0,levels=design)
fit2 <- limma::contrasts.fit(fit,contrasts=contrasts)
fit2 <- limma::eBayes(fit2,trend=TRUE)
#summary(decideTests(fit2, method="global"))
##
tT <- limma::topTable(fit2, adjust.method = "fdr", number = Inf,coef=1)
if(nrow(tT)==1){
rownames(tT) <- rownames(new_mat)
}
tT <- base::cbind(ID=rownames(tT),tT,stringsAsFactors=FALSE)
tT <- tT[rownames(new_mat),,drop=FALSE]
exp_G1 <- base::rowMeans(new_mat[,G1,drop=FALSE]);
exp_G0 <- base::rowMeans(new_mat[,G0,drop=FALSE]);
w1 <- which(tT$P.Value<=0);
if(base::length(w1)>0) tT$P.Value[w1] <- .Machine$double.xmin;
#z_val <- sapply(tT$P.Value*sign(tT$logFC),function(x)combinePvalVector(x,twosided = TRUE)[1])
z_val <- sapply(tT$P.Value*sign(tT$logFC),function(x)ifelse(x ==0, 0, combinePvalVector(x,twosided = TRUE)[1])) ## remove zero
if(is.null(random_effect)==TRUE){
tT <- base::cbind(tT,'Z-statistics'=z_val,'Ave.G0'=exp_G0,'Ave.G1'=exp_G1)
}else{
tT <- base::cbind(tT,'Z-statistics'=z_val,'Ave.G0'=exp_G0,'Ave.G1'=exp_G1,
'Ave.G0_RemoveRandomEffect'=fit@.Data[[1]][rownames(tT),'G0'],
'Ave.G1_RemoveRandomEffect'=fit@.Data[[1]][rownames(tT),'G1'])
}
if(is.null(G0_name)==FALSE) colnames(tT) <- gsub('Ave.G0',paste0('Ave.',G0_name),colnames(tT))
if(is.null(G1_name)==FALSE) colnames(tT) <- gsub('Ave.G1',paste0('Ave.',G1_name),colnames(tT))
tT <- tT[order(tT$P.Value, decreasing = FALSE), ]
return(tT)
}
#' Combine P Values Using Fisher's Method or Stouffer's Method
#'
#' \code{combinePvalVector} is a function to combine multiple comparison's P values using Fisher's method or Stouffer's method.
#'
#' @param pvals a vector of numerics, the P values from multiple comparison need to be combined.
#' @param method character, users can choose between "Stouffer" and "Fisher". Default is "Stouffer".
#' @param signed logical, if TRUE, will give a sign to the P value to indicate the direction of testing.
#' Default is TRUE.
#' @param twosided logical, if TRUE, P value is calculated in a one-tailed test.
#' If FALSE, P value is calculated in a two-tailed test, and it falls within the range 0 to 0.5.
#' Default is TRUE.
#' @return Return a vector contains the "Z-statistics" and "P.Value".
#' @examples
#' combinePvalVector(c(0.1,1e-3,1e-5))
#' combinePvalVector(c(0.1,1e-3,-1e-5))
#' @export
combinePvalVector <-
function(pvals,
method = 'Stouffer',
signed = TRUE,
twosided = TRUE) {
#
all_input_para <- c('pvals','method','signed','twosided')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('signed',c(TRUE,FALSE),envir=environment()),
check_option('twosided',c(TRUE,FALSE),envir=environment()),
check_option('method',c('Stouffer','Fisher'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
#remove NA pvalues
pvals <- pvals[!is.na(pvals) & !is.null(pvals)]
pvals[which(abs(pvals)<=0)] <- .Machine$double.xmin
if (sum(is.na(pvals)) >= 1) {
stat <- NA
pval <- NA
} else{
if (twosided & (sum(pvals > 1 | pvals < -1) >= 1))
stop('pvalues must between 0 and 1!\n')
if (!twosided & (sum(pvals > 0.5 | pvals < -0.5) >= 1))
stop('One-sided pvalues must between 0 and 0.5!\n')
if (!signed) {
pvals <- abs(pvals)
}
signs <- sign(pvals)
signs[signs == 0] <- 1
if (grepl('Fisher', method, ignore.case = TRUE)) {
if (twosided & signed) {
neg.pvals <- pos.pvals <- abs(pvals) / 2
pos.pvals[signs < 0] <- 1 - pos.pvals[signs < 0]
neg.pvals[signs > 0] <- 1 - neg.pvals[signs > 0]
} else{
neg.pvals <- pos.pvals <- abs(pvals)
}
pvals <-
c(1, -1) * c(
pchisq(
-2 * sum(log(as.numeric(pos.pvals))),
df = 2 * base::length(pvals),
lower.tail = FALSE
) / 2,
pchisq(
-2 * sum(log(as.numeric(neg.pvals))),
df = 2 * base::length(pvals),
lower.tail = FALSE
) / 2
)
pval <- base::min(abs(pvals))[1]
#if two pvals are equal, pick up the first one
stat <-
sign(pvals[abs(pvals) == pval])[1] * qnorm(pval, lower.tail = F)[1]
pval <- 2 * pval
}
else if (grepl('Stou', method, ignore.case = TRUE)) {
if (twosided) {
zs <- signs * qnorm(abs(pvals) / 2, lower.tail = FALSE)
stat <- sum(zs) / sqrt(base::length(zs))
pval <- 2 * pnorm(abs(stat), lower.tail = FALSE)
}
else{
zs <- signs * qnorm(abs(pvals), lower.tail = FALSE)
stat <- sum(zs) / sqrt(base::length(zs))
pval <- pnorm(abs(stat), lower.tail = FALSE)
}
}
else{
stop('Only \"Fisher\" or \"Stouffer\" method is supported!!!\n')
}
}
return(c(`Z-statistics` = stat, `P.Value` = pval))
}
#' Merge Activity Values from TF (transcription factors) ExpressionSet Object and Sig (signaling factors) ExpressionSet Object
#'
#' \code{merge_TF_SIG.AC} combines the activity value from TF (transcription factors) and Sig (signaling factors) ExpressionSet objects together,
#' and adds "_TF" or "_SIG" suffix to drivers for easier distinction.
#'
#'
#' @param TF_AC ExpressionSet object, containing the activity values for all TFs.
#' @param SIG_AC ExpressionSet object, containing the activity values for all SIGs.
#'
#' @return Return an ExpressionSet object.
#' @examples
#' if(exists('analysis.par')==TRUE) rm(analysis.par)
#' network.dir <- sprintf('%s/demo1/network/',system.file(package = "NetBID2")) # use demo
#' network.project.name <- 'project_2019-02-14' # demo project name
#' project_main_dir <- 'test/'
#' project_name <- 'test_driver'
#' analysis.par <- NetBID.analysis.dir.create(project_main_dir=project_main_dir,
#' project_name=project_name,
#' network_dir=network.dir,
#' network_project_name=network.project.name)
#' analysis.par$tf.network <- get.SJAracne.network(network_file=analysis.par$tf.network.file)
#' analysis.par$sig.network <- get.SJAracne.network(network_file=analysis.par$sig.network.file)
#' ## get eset (here for demo, use network.par$net.eset)
#' network.par <- list()
#' network.par$out.dir.DATA <- system.file('demo1','network/DATA/',package = "NetBID2")
#' NetBID.loadRData(network.par=network.par,step='exp-QC')
#' analysis.par$cal.eset <- network.par$net.eset
#' ac_mat_TF <- cal.Activity(target_list=analysis.par$tf.network$target_list,
#' cal_mat=Biobase::exprs(analysis.par$cal.eset),
#' es.method='weightedmean')
#' ac_mat_SIG <- cal.Activity(target_list=analysis.par$tf.network$target_list,
#' cal_mat=Biobase::exprs(analysis.par$cal.eset),
#' es.method='weightedmean')
#' analysis.par$ac.tf.eset <- generate.eset(exp_mat=ac_mat_TF,
#' phenotype_info=Biobase::pData(analysis.par$cal.eset))
#' analysis.par$ac.sig.eset <- generate.eset(exp_mat=ac_mat_SIG,
#' phenotype_info=Biobase::pData(analysis.par$cal.eset))
#' analysis.par$merge.ac.eset <- merge_TF_SIG.AC(TF_AC=analysis.par$ac.tf.eset,
#' SIG_AC=analysis.par$ac.sig.eset)
#' @export
merge_TF_SIG.AC <- function(TF_AC=NULL,SIG_AC=NULL){
#
all_input_para <- c('TF_AC','SIG_AC')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
mat_TF <- Biobase::exprs(TF_AC)
mat_SIG <- Biobase::exprs(SIG_AC)
funcType <- c(rep('TF',nrow(mat_TF)),rep('SIG',nrow(mat_SIG)))
rn <- c(rownames(mat_TF),rownames(mat_SIG))
rn_label <- base::paste(rn,funcType,sep='_')
mat_combine <- base::rbind(mat_TF,mat_SIG[,colnames(mat_TF)])
rownames(mat_combine) <- rn_label
eset_combine <- generate.eset(exp_mat=mat_combine,phenotype_info=Biobase::pData(TF_AC)[colnames(mat_combine),],
feature_info=NULL,annotation_info='activity in dataset')
return(eset_combine)
}
#' Merge TF (transcription factor) Network and Sig (signaling factor) Network
#'
#' \code{merge_TF_SIG.network} takes TF network and Sig network and combine them together.
#' The merged list object contains three elements, a data.frame contains all the combined network information \code{network_dat},
#' a driver-to-target list object \code{target_list}, and an igraph object of the network \code{igraph_obj}.
#'
#' @param TF_network list, the TF network created by \code{get.SJAracne.network} function.
#' @param SIG_network list, the SIG network created by \code{get.SJAracne.network} function.
#' @return
#' Return the a list containing three elements, \code{network_dat}, \code{target_list} and \code{igraph_obj}.
#' @examples
#' if(exists('analysis.par')==TRUE) rm(analysis.par)
#' network.dir <- sprintf('%s/demo1/network/',system.file(package = "NetBID2")) # use demo
#' network.project.name <- 'project_2019-02-14' # demo project name
#' project_main_dir <- 'test/'
#' project_name <- 'test_driver'
#' analysis.par <- NetBID.analysis.dir.create(project_main_dir=project_main_dir,
#' project_name=project_name,
#' network_dir=network.dir,
#' network_project_name=network.project.name)
#' analysis.par$tf.network <- get.SJAracne.network(network_file=analysis.par$tf.network.file)
#' analysis.par$sig.network <- get.SJAracne.network(network_file=analysis.par$sig.network.file)
#' analysis.par$merge.network <- merge_TF_SIG.network(TF_network=analysis.par$tf.network,
#' SIG_network=analysis.par$sig.network)
#' @export
merge_TF_SIG.network <- function(TF_network=NULL,SIG_network=NULL){
#
all_input_para <- c('TF_network','SIG_network')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
s_TF <- names(TF_network$target_list)
s_SIG <- names(SIG_network$target_list)
funcType <- c(rep('TF',base::length(s_TF)),rep('SIG',base::length(s_SIG)))
rn <- c(s_TF,s_SIG)
rn_label <- base::paste(rn,funcType,sep='_')
target_list_combine <- c(TF_network$target_list,SIG_network$target_list)
names(target_list_combine) <- rn_label
n_TF <- TF_network$network_dat
if(nrow(n_TF)>0) n_TF$source <- base::paste(n_TF$source,'TF',sep='_')
n_SIG <- SIG_network$network_dat
if(nrow(n_SIG)>0) n_SIG$source <- base::paste(n_SIG$source,'SIG',sep='_')
net_dat <- base::rbind(n_TF,n_SIG)
igraph_obj <- graph_from_data_frame(net_dat[,c('source','target')],directed=TRUE)
if('MI' %in% colnames(net_dat)) igraph_obj <- set_edge_attr(igraph_obj,'weight',index=E(igraph_obj),value=net_dat[,'MI'])
if('spearman' %in% colnames(net_dat)) igraph_obj <- set_edge_attr(igraph_obj,'sign',index=E(igraph_obj),value=sign(net_dat[,'spearman']))
return(list(network_dat=net_dat,target_list=target_list_combine,igraph_obj=igraph_obj))
}
#' Generate the Master Table for Drivers
#'
#' \code{generate.masterTable} generates a master table to show the mega information of all tested drivers.
#'
#' The master table gathers TF (transcription factor) information, Sig (signaling factor) information, all the DE (differential expression analysis)
#' and DA (differential activity analysis) from multiple comparisons. It also shows each driver's target gene size and other additional information
#' (e.g. gene biotype, chromosome name, position etc.).
#'
#' @param use_comp a vector of characters, the name of multiple comparisons. It will be used to name the columns of master table.
#' @param DE list, a list of DE comparisons, each comparison is a data.frame. The element name in the list must contain the name in \code{use_comp}.
#' @param DA list, a list of DA comparisons, each comparison is a data.frame. The element name in the list must contain the name in \code{use_comp}.
#' @param target_list list, a driver-to-target list. The names of the list elements are drivers. Each element is a data frame, usually contains three columns.
#' "target", target gene names; "MI", mutual information; "spearman", spearman correlation coefficient.
#' It is highly suggested to follow the NetBID2 pipeline, and the \code{TF_network} could be generated by \code{get_net2target_list} and \code{get.SJAracne.network}.
#' @param main_id_type character, the type of driver's ID. It comes from the attribute name in biomaRt package.
#' Such as "ensembl_gene_id", "ensembl_gene_id_version", "ensembl_transcript_id", "ensembl_transcript_id_version" or "refseq_mrna".
#' For details, user can call \code{biomaRt::listAttributes()} to display all available attributes in the selected dataset.
#' @param transfer_tab data.frame, the data frame for ID conversion. This can be obtained by calling \code{get_IDtransfer}.
#' If NULL and \code{main_id_type} is not in the column names of \code{tf_sigs}, it will use the conversion table within the function.
#' Default is NULL.
#' @param tf_sigs list, contains all the detailed information of TF and Sig. Users can call \code{db.preload} for access.
#' @param z_col character, name of the column in \code{DE} and \code{DA} contains the Z statistics. Default is "Z-statistics".
#' @param display_col character, name of the column in \code{DE} and \code{DA} need to be kept in the master table. Default is c("logFC","P.Value").
#' @param column_order_strategy character, users can choose between "type" and "comp". Default is "type".
#' If set as type, the columns will be ordered by column type; If set as comp, the columns will be ordered by comparison.
#' @return Return a data frame contains the mega information of all tested drivers.
#' The column "originalID" and "originalID_label" is the same ID as from the original dataset.
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' #analysis.par$final_ms_tab ## this is master table generated before
#' ac_mat <- cal.Activity(target_list=analysis.par$merge.network$target_list,
#' cal_mat=Biobase::exprs(analysis.par$cal.eset),es.method='weightedmean')
#' analysis.par$ac.merge.eset <- generate.eset(exp_mat=ac_mat,
#' phenotype_info=Biobase::pData(analysis.par$cal.eset))
#' phe_info <- Biobase::pData(analysis.par$cal.eset)
#' all_subgroup <- base::unique(phe_info$subgroup) ##
#' for(each_subtype in all_subgroup){
#' comp_name <- sprintf('%s.Vs.others',each_subtype) ## each comparison must give a name !!!
#' G0 <- rownames(phe_info)[which(phe_info$`subgroup`!=each_subtype)] # get sample list for G0
#' G1 <- rownames(phe_info)[which(phe_info$`subgroup`==each_subtype)] # get sample list for G1
#' DE_gene_limma <- getDE.limma.2G(eset=analysis.par$cal.eset,G1=G1,G0=G0,
#' G1_name=each_subtype,G0_name='other')
#' analysis.par$DE[[comp_name]] <- DE_gene_limma
#' DA_driver_limma <- getDE.limma.2G(eset=analysis.par$ac.merge.eset,G1=G1,G0=G0,
#' G1_name=each_subtype,G0_name='other')
#' analysis.par$DA[[comp_name]] <- DA_driver_limma
#' }
#' all_comp <- names(analysis.par$DE) ## get all comparison name for output
#' db.preload(use_level='gene',use_spe='human',update=FALSE);
#' test_ms_tab <- generate.masterTable(use_comp=all_comp,
#' DE=analysis.par$DE,
#' DA=analysis.par$DA,
#' target_list=analysis.par$merge.network$target_list,
#' tf_sigs=tf_sigs,
#' z_col='Z-statistics',
#' display_col=c('logFC','P.Value'),
#' main_id_type='external_gene_name')
#' @export
generate.masterTable <- function(use_comp=NULL,DE=NULL,DA=NULL,
target_list=NULL,main_id_type=NULL,transfer_tab=NULL,
tf_sigs=NULL,
z_col='Z-statistics',display_col=c('logFC','P.Value'),
column_order_strategy='type'){
#
all_input_para <- c('use_comp','DE','DA','target_list','main_id_type','tf_sigs','z_col','display_col','column_order_strategy')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=base::environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('column_order_strategy',c('type','comp'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(base::length(base::setdiff(use_comp,names(DE)))>0){message(sprintf('%s not included in DE, please check and re-try!',base::setdiff(use_comp,names(DE))));return(FALSE)}
if(base::length(base::setdiff(use_comp,names(DA)))>0){message(sprintf('%s not included in DA, please check and re-try!',base::setdiff(use_comp,names(DA))));return(FALSE)}
# get original ID
ori_rn <- rownames(DA[[1]]) ## with label
w1 <- grep('(.*)_TF',ori_rn); w2 <- grep('(.*)_SIG',ori_rn)
funcType <- rep(NA,length.out=base::length(ori_rn));rn <- funcType
funcType[w1] <- 'TF'; funcType[w2] <- 'SIG';
rn[w1] <- gsub('(.*)_TF',"\\1",ori_rn[w1]);rn[w2] <- gsub('(.*)_SIG',"\\1",ori_rn[w2]);
rn_label <- ori_rn
use_size <- unlist(lapply(target_list[rownames(DA[[1]])],nrow))
# id issue
current_id <- names(tf_sigs$tf)[-1]
use_info <- base::unique(base::rbind(tf_sigs$tf$info,tf_sigs$sig$info))
if(main_id_type %in% current_id){
use_info <- use_info[which(use_info[,main_id_type] %in% rn),]
}else{
if(is.null(transfer_tab)==TRUE){
transfer_tab <- get_IDtransfer(from_type=main_id_type,to_type=current_id[1],use_genes=rn,ignore_version = TRUE)
use_info <- base::merge(use_info,transfer_tab,by.x=current_id[1],by.y=current_id[1])
}else{
transfer_tab <- transfer_tab[which(transfer_tab[,main_id_type] %in% rn),]
uid <- base::intersect(colnames(transfer_tab),colnames(use_info))
if(base::length(uid)==0){message('No ID type in the transfer_tab could match ID type in tf_sigs, please check and re-try!');return(FALSE);}
uid <- uid[1]
use_info <- base::merge(use_info,transfer_tab,by.x=uid,by.y=uid)
}
use_info <- use_info[which(use_info[,main_id_type] %in% rn),]
}
use_info <- base::unique(use_info)
if(nrow(use_info)==0){message('ID issue error, please check main_id_type setting!');return(FALSE);}
# merge info
tmp1 <- stats::aggregate(use_info,list(use_info[,main_id_type]),function(x){
x1 <- x[which(x!="")]
x1 <- x1[which(is.na(x1)==FALSE)]
base::paste(sort(base::unique(x1)),collapse=';')
})
tmp1 <- tmp1[,-1]; rownames(tmp1) <- tmp1[,main_id_type]
geneSymbol <- tmp1[rn,'external_gene_name'] ## this column for function enrichment
if('external_transcript_name' %in% colnames(tmp1)){ ## this column for display
gene_label <- base::paste(tmp1[rn,'external_transcript_name'],funcType,sep = '_')
}else{
gene_label <-base::paste(tmp1[rn,'external_gene_name'],funcType,sep = '_')
}
#
#label_info <- data.frame('gene_label'=gene_label,'geneSymbol'=geneSymbol,
# 'originalID'=rn,'originalID_label'=rn_label,'funcType'=funcType,'Size'=use_size,stringsAsFactors=FALSE)
label_info <- data.frame('originalID_label'=rn_label,'originalID'=rn,'gene_label'=gene_label,'geneSymbol'=geneSymbol,
'funcType'=funcType,'Size'=use_size,stringsAsFactors=FALSE)
w1 <- which(is.na(geneSymbol)==TRUE)
label_info[w1,'geneSymbol'] <- label_info[w1,'originalID']
label_info[w1,'gene_label'] <- label_info[w1,'originalID_label']
add_info <- tmp1[rn,]
#
combine_info <- lapply(use_comp,function(x){
DA[[x]] <- DA[[x]][rn_label,,drop=F]
DE[[x]] <- as.data.frame(DE[[x]])[rn,]
avg_col <- colnames(DA[[x]])[grep('^Ave',colnames(DA[[x]]))]
uc <- c(z_col,avg_col,base::setdiff(display_col,c(z_col,avg_col))); uc <- base::intersect(uc,colnames(DA[[x]]))
DA_info <- DA[[x]][rn_label,uc,drop=F]
avg_col <- colnames(DE[[x]])[grep('^Ave',colnames(DE[[x]]))]
uc <- c(z_col,avg_col,base::setdiff(display_col,c(z_col,avg_col))); uc <- base::intersect(uc,colnames(DE[[x]]))
DE_info <- as.data.frame(DE[[x]])[rn,uc,drop=F]
colnames(DA_info) <- paste0(colnames(DA_info),'.',x,'_DA')
colnames(DE_info) <- paste0(colnames(DE_info),'.',x,'_DE')
colnames(DA_info)[1] <- paste0('Z.',x,'_DA')
colnames(DE_info)[1] <- paste0('Z.',x,'_DE')
out <- base::cbind(DA_info,DE_info,stringsAsFactors=FALSE)
rownames(out) <- rn_label
out
})
combine_info_DA <- do.call(base::cbind,lapply(combine_info,function(x)x[rn_label,grep('_DA$',colnames(x)),drop=T]))
combine_info_DE <- do.call(base::cbind,lapply(combine_info,function(x)x[rn_label,grep('_DE$',colnames(x)),drop=T]))
# re-organize the columns for combine info
if(column_order_strategy=='type' & length(use_comp)>1){
col_ord <- c('Z','AveExpr',display_col)
tmp1 <- lapply(col_ord,function(x){
x1 <- grep(sprintf('^%s\\.',x),colnames(combine_info_DA))
if(length(x1)>0) combine_info_DA[,x1] else return(NULL)
})
combine_info_DA <- do.call(base::cbind,tmp1)
tmp1 <- lapply(col_ord,function(x){
combine_info_DE[,grep(sprintf('^%s\\.',x),colnames(combine_info_DE))]
})
combine_info_DE <- do.call(base::cbind,tmp1)
}
# put them together
ms_tab <- base::cbind(label_info,combine_info_DA,combine_info_DE,add_info)
rownames(ms_tab) <- ms_tab$originalID_label
return(ms_tab)
}
#' Save the Master Table into Excel File
#'
#' \code{out2excel} is a function can save data frame as Excel File. This is mainly for the output of master table generated by \code{generate.masterTable}.
#'
#' @param all_ms_tab list or data.frame, if data.frame, it is generated by \code{generate.masterTable}.
#' If list, each list element is data.frame/master table.
#' The name of the list element will be the sheet name in the excel file.
#' @param out.xlsx character, path and file name of the output Excel file.
#' @param mark_gene list, list of marker genes. The name of the list element is the marked group name. Each element is a vector of marker genes.
#' This is optional, just to add additional information to the file.
#' @param mark_col character, the color to mark the marker genes. If NULL, will use \code{get.class.color} to get the colors.
#' @param mark_strategy character, users can choose between "color" and "add_column".
#' "Color" means the mark_gene will be marked by filling its background color;
#' "add_column" means the mark_gene will be displayed in a separate column with TRUE/FALSE, indicating whether the gene belongs to a mark group or not.
#' @param workbook_name character, name of the workbook for the output Excel. Default is "ms_tab".
#' @param only_z_sheet logical, if TRUE, will create a separate sheet only contains Z-statistics related columns from DA/DE analysis.
#' Default is FALSE.
#' @param z_column character, name of the columns contain Z-statistics. If NULL, find column names start with "Z.".
#' Default is NULL.
#' @param sig_thre numeric, threshold for the Z-statistics. Z values passed the threshold will be colored. The color scale is defined by \code{z2col}.
#' Default is 1.64.
#' @return Return a logical value. If TRUE, the Excel file has been generated successfully.
#' @examples
#' \dontrun{
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' ms_tab <- analysis.par$final_ms_tab ## this is master table generated before
#' mark_gene <- list(WNT=c('WIF1','TNC','GAD1','DKK2','EMX2'),
#' SHH=c('PDLIM3','EYA1','HHIP','ATOH1','SFRP1'),
#' Group3=c('IMPG2','GABRA5','EGFL11','NRL','MAB21L2','NPR3','MYC'),
#' Group4=c('KCNA1','EOMES','KHDRBS2','RBM24','UNC5D'))
#' mark_col <- get.class.color(names(mark_gene),
#' pre_define=c('WNT'='blue','SHH'='red',
#' 'Group3'='yellow','Group4'='green'))
#' outfile <- 'test_out.xlsx'
#' out2excel(ms_tab,out.xlsx = outfile,mark_gene,mark_col)
#' }
#' @export
out2excel <- function(all_ms_tab,out.xlsx,
mark_gene=NULL,
mark_col=NULL,
mark_strategy='color',
workbook_name='ms_tab',
only_z_sheet=FALSE,
z_column=NULL,sig_thre=1.64){
#
all_input_para <- c('all_ms_tab','out.xlsx','mark_strategy','workbook_name','only_z_sheet','sig_thre')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('only_z_sheet',c(TRUE,FALSE),envir=environment()),
check_option('mark_strategy',c('color','add_column'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
wb <- openxlsx::createWorkbook(workbook_name)
if(!'list' %in% class(all_ms_tab)){
all_ms_tab <- list('Sheet1'=as.data.frame(all_ms_tab))
}
if(only_z_sheet==TRUE){
nn <- names(all_ms_tab)
all_ms_tab <- lapply(all_ms_tab,function(x){
w1 <- grep('.*_[DE|DA]',colnames(x))
w2 <- base::setdiff(colnames(x)[w1],colnames(x)[w1][grep('^Z.*',colnames(x)[w1])])
w3 <- base::setdiff(colnames(x),w2)
list(x[,w3],x)
})
all_ms_tab <- unlist(all_ms_tab,recursive = FALSE)
nn1 <- lapply(nn,function(x){
if(x!='Sheet1'){
c(sprintf('Only_Z_%s',x),sprintf('Full_info_%s',x))
}else{
c('Only_Z','Full_info')
}
})
nn1 <- unlist(nn1)
names(all_ms_tab) <- nn1
}
if(is.null(mark_gene)==FALSE){
if(!'list' %in% class(mark_gene)){
mark_gene <- list('mark_gene'=mark_gene)
}
if(is.null(names(mark_gene))==TRUE){
message('Must give name to the mark_gene list');return(FALSE)
}
if(mark_strategy=='color' & is.null(mark_col)==TRUE){
mark_col <- get.class.color(names(mark_gene))
}
if(mark_strategy=='add_column'){
new_col_name <- paste0('is',names(mark_gene))
all_ms_tab <- lapply(all_ms_tab,function(x){
g1 <- x$geneSymbol
r1 <- do.call(base::cbind,lapply(mark_gene,function(x1){
ifelse(g1 %in% x1,'TRUE','FALSE')
}))
new_x <- base::cbind(x,r1)
colnames(new_x) <- c(colnames(x),new_col_name)
new_x
})
}
}
z_column_index <- 'defined'
if(is.null(z_column)==TRUE) z_column_index <- 'auto'
i <- 0
headerStyle <- openxlsx::createStyle(fontSize = 14, fontColour = "#FFFFFF", halign = "center",fgFill = "#4F81BD",
border="TopBottom", borderColour = "#4F81BD",wrapText=TRUE) ## style for header line
for(sheetname in names(all_ms_tab)){ ## list, each item create one sheet
i <- i +1
d <- as.data.frame(all_ms_tab[[sheetname]])
if(z_column_index=='auto') use_z_column <- colnames(d)[grep('^Z\\.',colnames(d),ignore.case = TRUE)] else use_z_column <- base::intersect(colnames(d),z_column)
openxlsx::addWorksheet(wb,sheetName=sheetname)
openxlsx::writeData(wb,sheet = i,d)
openxlsx::addStyle(wb, sheet = i, headerStyle, rows = 1, cols = 1:ncol(d), gridExpand = TRUE) ## add header style
all_c <- list()
for(z_col in use_z_column){ ## find colnames with z. (only applied to pipeline excel)
j <- which(colnames(d)==z_col)
z1 <- d[,j]; c1 <- z2col(z1,sig_thre=sig_thre)
all_c[[as.character(j)]] <- c1
}
mat_c <- do.call(base::cbind,all_c)
uni_c <- base::unique(unlist(all_c))
for(r in uni_c){
w1 <- which(mat_c==r)
nn <- nrow(mat_c)
rr <- w1%%nn+1
cc <- w1%/%nn+1
cc[which(rr==1)] <- cc[which(rr==1)]-1
rr[which(rr==1)] <- nn+1
openxlsx::addStyle(wb, sheet = i, createStyle(fgFill=r), rows =rr, cols = as.numeric(names(all_c)[cc]))
}
if(is.null(mark_gene)==FALSE){
for(k in names(mark_gene)){
openxlsx::addStyle(wb, sheet = i, openxlsx::createStyle(fgFill=mark_col[k]),
rows =which(toupper(d[,which(colnames(d)=='geneSymbol')]) %in% toupper(mark_gene[[k]]))+1,
cols = which(colnames(d)=='geneSymbol')) ## find column with geneSymbol and mark color
}
}
w1 <- which(gsub("\\s","",as.matrix(d))=='TRUE')
nn <- nrow(d)
rr <- w1%%nn+1
cc <- w1%/%nn+1
cc[which(rr==1)] <- cc[which(rr==1)]-1
rr[which(rr==1)] <- nn+1
openxlsx::addStyle(wb, sheet = i, openxlsx::createStyle(fontColour='#FF0000'), rows =rr, cols = as.numeric(cc)) ## find column with TRUE/FALSE and mark with color
}
openxlsx::saveWorkbook(wb, out.xlsx, overwrite = TRUE)
return(TRUE)
##
}
#' Load MSigDB Database into R Workspace
#'
#' \code{gs.preload} downloads data from MSigDB and stores it into two variables in R workspace, \code{all_gs2gene} and \code{all_gs2gene_info}.
#' \code{all_gs2gene} is a list object with elements of gene sets collections.
#' \code{all_gs2gene_info} is a data.frame contains the description of each gene sets.
#'
#' This is a pre-processing function for NetBID2 advanced analysis. User only need to input the species name (e.g. "Homo sapiens", "Mus musculus").
#' It will call \code{msigdbr} to download data from MSigDB and save it as RData under the \code{db/} directory with species name.
#'
#' @param use_spe character, name of interested species (e.g. "Homo sapiens", "Mus musculus").
#' Users can call \code{msigdbr_species()} to access the full list of available species names.
#' Default is "Homo sapiens".
#' @param update logical, if TRUE, the previous loaded RData will be updated. Default is FALSE.
#' @param main.dir character, the main file path of user's NetBID2 project.
#' If NULL, will be set to \code{system.file(package = "NetBID2")}. Default is NULL.
#' @param db.dir character, the file path to save the RData. Default is \code{db} directory under the \code{main.dir}, if one has a \code{main.dir}.
#'
#' @return Reture a logical value. If TRUE, MsigDB database is loaded successfully, with \code{all_gs2gene} and \code{all_gs2gene_info} created
#' in the workspace.
#'
#' @examples
#' gs.preload(use_spe='Homo sapiens',update=FALSE)
#' gs.preload(use_spe='Mus musculus',update=FALSE)
#' print(all_gs2gene_info)
#' # contain the information for all gene set collection, collection info, collection size,
#' ## sub collection,sub collection info, sub collection size
#' print(names(all_gs2gene)) # the first level of the list is the collection and sub-collection IDs
#' print(str(all_gs2gene$`CP:KEGG_MEDICUS`))
#'
#' \dontrun{
#' gs.preload(use_spe='Homo sapiens',update=TRUE)
#' }
#' @export
gs.preload <- function(use_spe = 'Homo sapiens',
update = FALSE,
main.dir = NULL,
db.dir = sprintf("%s/db/",main.dir)){
#
all_input_para <- c('use_spe','update')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
all_spe <- msigdbr::msigdbr_species()[["species_name"]]
check_res <- c(check_option('update',c(TRUE,FALSE),envir=environment()),
check_option('use_spe',all_spe,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
## only support geneSymbol (because pipeline-generated master table will contain geneSymbol column)
if(is.null(main.dir)==TRUE){
main.dir <- system.file(package = "NetBID2")
message(sprintf('main.dir not set, will use package directory: %s',main.dir))
}
if(is.null(db.dir)==TRUE){
db.dir <- sprintf("%s/db/",main.dir)
}
message(sprintf('Will use directory %s as the db.dir',db.dir))
use_spe1 <- gsub(' ','_',use_spe)
out_file <- sprintf('%s/%s_gs2gene.RData',db.dir,use_spe1)
if(file.exists(out_file)==FALSE | update==TRUE){
message('Begin generating all_gs2gene !')
all_gs_info <- msigdbr::msigdbr(species = use_spe) ## use msigdbr_species() to check possible available species
# for gs_collection
all_gs_cat <- base::unique(all_gs_info$gs_collection)
all_gs2gene_1 <- lapply(all_gs_cat,function(x){
x1 <- all_gs_info[which(all_gs_info$gs_collection==x),]
all_gs <- base::unique(x1$gs_name)
x2 <- lapply(all_gs, function(y){
base::unique(x1$gene_symbol[which(x1$gs_name==y)])
})
names(x2) <- all_gs;x2
})
names(all_gs2gene_1) <- all_gs_cat
all_gs_subcat <- base::setdiff(base::unique(all_gs_info$gs_subcollection),"")
all_gs2gene_2 <- lapply(all_gs_subcat,function(x){
x1 <- all_gs_info[which(all_gs_info$gs_subcollection==x),]
all_gs <- base::unique(x1$gs_name)
x2 <- lapply(all_gs, function(y){
base::unique(x1$gene_symbol[which(x1$gs_name==y)])
})
names(x2) <- all_gs;x2
})
names(all_gs2gene_2) <- all_gs_subcat
all_gs2gene <- c(all_gs2gene_1,all_gs2gene_2)
#
all_gs2gene <- all_gs2gene[sort(names(all_gs2gene))]
gs_size <- unlist(lapply(all_gs2gene,length))
# info for cat
info_cat <- c('C1'='positional gene sets', 'C2'='curated gene sets', 'C3'='regulatory target gene sets', 'C4'='computational gene sets', 'C5'='ontology gene sets',
'C6'='oncogenic signature gene sets', 'C7'='immunologic signature gene sets', 'C8'='cell type signature gene sets', 'H'='hallmark gene sets')
info_subcat <- c('CGP'='Chemical and Genetic Perturbations', 'CP'='Canonical Pathways', 'CP:BIOCARTA'='BioCarta subset', 'CP:KEGG_MEDICUS'='KEGG MEDICUS subset', 'CP:KEGG_LEGACY'='KEGG LEGACY subset', 'CP:REACTOME'='Reactome subset', 'CP:WIKIPATHWAYS'='WikiPahtways subset', 'CP:PID'='PID subset',
'MIR:MIRDB'='microRNA targets predicted by miRDB v6.0', 'MIR:MIR_LEGACY'='microRNA targets predicted by motif matching', 'TFT:GTRD'='Transcription Factor Targets mined from GTRD', 'TFT:TFT_LEGACY'='Transcription Factor Targets predicted by motif matching',
'3CA'='Computational gene sets mined from Curated Cancer Cell Atlas (3CA) metaprograms', 'CGN'='Computational gene sets defined by expression neighborhoods centered on 380 cancer-associated genes', 'CM'='Computational gene sets compiled from cancer modules significantly changed in a variety of cancer conditions',
'GO:BP'='GO Biological Process', 'GO:MF'='GO Molecular Function', 'GO:CC'='GO Cellular Component', 'HPO'='Human Phenotype Ontology', 'VAX'='Immune signatures curated from vaccine response studies', 'IMMUNESIGDB'='Immune signatures collected from ImmuneSigDB')
cat_rel <- base::unique(as.data.frame(all_gs_info[,c('gs_collection','gs_subcollection')]))
all_gs2gene_info <- data.frame(cat_rel[,1],info_cat[cat_rel[,1]],gs_size[cat_rel[,1]],cat_rel[,2],info_subcat[cat_rel[,2]],gs_size[cat_rel[,2]],stringsAsFactors = FALSE)
colnames(all_gs2gene_info) <- c('Collection','Collection_Info','Collection_Size','Subcollection','Subcollection_Info','Subcollection_Size')
all_gs2gene_info <- all_gs2gene_info[order(all_gs2gene_info[,1]),]
all_gs2gene_info[,c(1,2,4,5)] <- as.data.frame(apply(all_gs2gene_info[,c(1,2,4,5)],2,function(x){x[which(is.na(x)==TRUE)] <- "";x}),stringsAsFactors=FALSE)
save(all_gs2gene,all_gs2gene_info,file=out_file)
}
load(out_file,.GlobalEnv)
message('all_gs2gene loaded, you could see all_gs2gene_info to check the details !')
return(TRUE)
}
######################################################### visualization functions
## simple functions to get info
#' Create a vector of each sample's selected phenotye descriptive information.
#'
#' \code{get_obs_label} creates a vector of each sample's selected phenotype descriptive information.
#' This is a helper function for data visualization.
#'
#' @param phe_info data.frame, the phenotype data of the samples.
#' It is a data frame that can store any number of descriptive columns (covariates) for each sample row.
#' To get the phenotype data, using the accessor function \code{pData}.
#' @param use_col a vector of numerics or characters.
#' Users can select the interested descriptive column(s) by calling index or name of the column(s).
#' @param collapse character, an optional character string to separate the results when the length
#' of \code{use_col} is more than 1. Not NA_character. Default is "|".
#'
#' @return
#' Return a vector of selected phenotype descriptive information (covariates) for each sample.
#' Vector name is the sample name.
#' @examples
#' analysis.par <- list()
#' analysis.par$out.dir.DATA <- system.file('demo1','driver/DATA/',package = "NetBID2")
#' NetBID.loadRData(analysis.par=analysis.par,step='ms-tab')
#' phe_info <- Biobase::pData(analysis.par$cal.eset)
#' use_obs_class <- get_obs_label(phe_info = phe_info,'subgroup')
#' print(use_obs_class)
#' \dontrun{
#'}
#' @export
get_obs_label <- function(phe_info,use_col,collapse='|'){
#
all_input_para <- c('phe_info','use_col','collapse')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
w1 <- base::setdiff(use_col,colnames(phe_info))
if(length(w1)>0){
message(sprintf('%s not in the colnames(phe_info),please check and re-try!',paste(w1,collapse=';')));return(FALSE);
}
obs_label<-phe_info[,use_col];
if(base::length(use_col)>1){
obs_label<-apply(obs_label,1,function(x)base::paste(x,collapse=collapse))
}
names(obs_label) <- rownames(phe_info);
obs_label
}
#' Get interested phenotype groups from pData slot of the ExpressionSet object.
#'
#' \code{get_int_group} is a function to extract interested phenotype groups from the ExpressionSet object
#' with 'cluster-meaningful' sample features.
#'
#' @param eset an ExpressionSet object.
#' @return Return a vector of phenotype groups which could be used for sample cluster analysis.
#'
#' @examples
#' network.par <- list()
#' network.par$out.dir.DATA <- system.file('demo1','network/DATA/',package = "NetBID2")
#' NetBID.loadRData(network.par=network.par,step='exp-QC')
#' intgroups <- get_int_group(network.par$net.eset)
#' @export
get_int_group <- function(eset){
all_input_para <- c('eset')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
phe <- Biobase::pData(eset)
feature_len <- apply(phe,2,function(x)base::length(base::unique(x)))
intgroup <- colnames(phe)[which(feature_len>1 & feature_len<nrow(phe))]
return(intgroup)
}
#' Get Score to Measure Similarity Between Observed Classification and Predicted Classification.
#'
#' \code{get_clustComp} calculates a score to measure the similarity between two classifications.
#'
#' @param pred_label a vector of characters, the predicted classification labels.
#' @param obs_label a vector of characters, the observed classification labels.
#' @param strategy character, the method applied to calculate the score.
#' Users can choose "ARI (adjusted rand index)", "NMI (normalized mutual information)" or "Jaccard".
#' Default is "ARI".
#' @return Return a score for the measurement of similarity.
#' @examples
#' obs_label <- c('A','A','A','B','B','C','D')
#' pred_label <- c(1,1,1,1,2,2,2)
#' get_clustComp(pred_label,obs_label)
#' @export
get_clustComp <- function(pred_label, obs_label,strategy='ARI') {
#
all_input_para <- c('pred_label','obs_label','strategy')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('strategy',c('ARI','NMI','Jaccard'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(names(pred_label))==TRUE & is.null(names(obs_label))==FALSE){
message('The names of pred_label will use the order of obs_label')
names(pred_label) <- names(obs_label)
}
if(is.null(names(obs_label))==TRUE & is.null(names(pred_label))==FALSE){
message('The names of obs_label will use the order of pred_label')
names(obs_label) <- names(pred_label)
}
if(is.null(names(obs_label))==TRUE & is.null(names(pred_label))==TRUE){
message('Assume pred_label and obs_label have the same order!')
names(obs_label) <- as.character(1:base::length(obs_label))
names(pred_label) <- names(obs_label)
}
if(strategy=='Jaccard') res1 <- get_jac(pred_label, obs_label) else res1 <- clustComp(pred_label, obs_label)[[strategy]]
return(res1)
}
# get jaccard accuracy
get_jac <- function(pred_label, obs_label) {
jac1 <- c()
for (i in base::unique(pred_label)) {
jac_index <- c()
x1 <- names(pred_label)[which(pred_label == i)]
for (j in base::unique(obs_label)) {
x2 <- names(obs_label)[which(obs_label == j)]
jac_index <-
c(jac_index, base::length(base::intersect(x1, x2)) / base::length(union(x1, x2)))
}
jac1 <- c(jac1, base::max(jac_index) * base::length(x1))
}
jac1 <- sum(jac1) / base::length(pred_label)
return(jac1)
}
#' Visualize Each Sample's Observed Label vs. Predicted Label in Table
#'
#' \code{draw.clustComp} draws a table to show each sample's observed label vs. its predicted label.
#' Each row represents an observed label (e.g. one subgroup of disease), each column represents the predicted label created by classification algorithm (e.g K-means).
#'
#' The table provides more details about the side-by-side PCA biplot created by \code{draw.emb.kmeans}.
#' The purpose is to find if any abnormal sample (outlier) exists. The darker the table cell is,
#' the more samples are gathered in the corresponding label.
#'
#' @param pred_label a vector of characters, the predicted labels created by classification (e.g K-means).
#' @param obs_label a vector of characters, the observed labels annotated by phenotype data.
#' @param strategy character, method to quantify the similarity between predicted labels vs. observed labels.
#' Users can choose from "ARI (adjusted rand index)", "NMI (normalized mutual information)" and "Jaccard".
#' Default is "ARI".
#' @param use_col logical, If TRUE, the table will be colored. The more sample gathered in one table cell, the darker shade it has.
#' Default is TRUE.
#' @param low_K integer, a threshold of sample number to be shown in a single cell.
#' If too many samples gathered in a single table cell, it will be challenging for eyes.
#' By setting the value of this threshold, if the number of samples gathered in one table cell exceeded the threshold,
#' only the number will be shown. Otherwise, all samples' names will be listed.
#' Default is 5.
#' @param highlight_clust a vector of characters, the predicted label need to be highlighted in the figure.
#' @param main character, an overall title for the plot.
#' @param clust_cex numeric, text size for the predicted label (column names). Default is 1.
#' @param outlier_cex numeric, text size for the observed label (row names). Default is 0.3.
#' @return Return a matrix of integers and a table for visualization. Rows are predicted label, columns are observed label.
#' Integer is the number of samples gathered in the corresponding label.
#' @examples
#' network.par <- list()
#' network.par$out.dir.DATA <- system.file('demo1','network/DATA/',package = "NetBID2")
#' NetBID.loadRData(network.par=network.par,step='exp-QC')
#' mat <- Biobase::exprs(network.par$net.eset)
#' phe <- Biobase::pData(network.par$net.eset)
#' intgroup <- 'subgroup'
#' pred_label <- draw.emb.kmeans(mat=mat,all_k = NULL,
#' obs_label=get_obs_label(phe,intgroup),
#' kmeans_strategy='consensus')
#' draw.clustComp(pred_label,get_obs_label(phe,intgroup),outlier_cex=1,low_K=2,use_col=TRUE)
#' draw.clustComp(pred_label,get_obs_label(phe,intgroup),outlier_cex=1,low_K=2,use_col=FALSE)
#' @export
draw.clustComp <- function(pred_label, obs_label,strategy='ARI',
use_col=TRUE,low_K=5,
highlight_clust=NULL,
main=NULL,clust_cex=1,outlier_cex=0.3) {
#
all_input_para <- c('pred_label','obs_label','strategy','use_col','low_K','clust_cex','outlier_cex')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
check_res <- c(check_option('use_col',c(TRUE,FALSE),envir=environment()),
check_option('strategy',c('ARI','NMI','Jaccard'),envir=environment()))
if(base::min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
if(is.null(names(pred_label))==TRUE & is.null(names(obs_label))==FALSE){
message('The names of pred_label will use the order of obs_label')
names(pred_label) <- names(obs_label)
}
if(is.null(names(obs_label))==TRUE & is.null(names(pred_label))==FALSE){
message('The names of obs_label will use the order of pred_label')
names(obs_label) <- names(pred_label)
}
if(is.null(names(obs_label))==TRUE & is.null(names(pred_label))==TRUE){
message('Assume pred_label and obs_label have the same order!')
names(obs_label) <- as.character(1:base::length(obs_label))
names(pred_label) <- names(obs_label)
}
nn <- names(pred_label)
k1 <- get_clustComp(pred_label,obs_label,strategy=strategy)
if(is.null(main)==TRUE){
mm <- sprintf('%s:%s',strategy,format(k1,digits=4),
format(k1,digits=4))
}else{
mm <- main
}
t1 <- base::table(list(pred_label[nn],obs_label[nn]))
graphics::layout(1)
textWidth <- base::max(strwidthMod(colnames(t1),units='inch',cex=clust_cex))+par.char2inch()[1]*1.5
par(mai=c(0.5,textWidth,1,1))
if(use_col==TRUE){
graphics::image(t1,col=c('white',grDevices::colorRampPalette(brewer.pal(8,'Reds'))(base::length(base::unique(as.numeric(t1)))-1)),bty='n',xaxt='n',yaxt='n',
main=mm)
}else{
graphics::image(t1,bty='n',xaxt='n',yaxt='n',
main=mm,col='white')
}
pp <- par()$usr
graphics::rect(xleft=pp[1],xright=pp[2],ybottom=pp[3],ytop=pp[4])
xx <- base::seq(pp[1],pp[2],length.out = nrow(t1)+1)
yy <- base::seq(pp[3],pp[4],length.out = ncol(t1)+1)
xxx <- (xx[1:(base::length(xx)-1)]+xx[2:base::length(xx)])/2
yyy <- (yy[1:(base::length(yy)-1)]+yy[2:base::length(yy)])/2
graphics::abline(h=yy);graphics::abline(v=xx)
graphics::text(pp[1]-par.char2pos()[1],yyy,colnames(t1),adj=1,xpd=TRUE,col=ifelse(colnames(t1) %in% highlight_clust,2,1),cex=clust_cex)
graphics::text(xxx,pp[4],rownames(t1),srt=0,xpd=TRUE,
col=ifelse(rownames(t1) %in% highlight_clust,2,1),cex=clust_cex,pos=3)
for(i in 1:nrow(t1)){
for(j in 1:ncol(t1)){
v1 <- t1[i,j]
if(v1==0) next
if(v1>low_K){
graphics::text(xxx[i],yyy[j],v1,cex=clust_cex)
}else{
v2 <- names(obs_label)[which(pred_label==rownames(t1)[i] & obs_label==colnames(t1)[j])]
v2 <- base::paste(v2,collapse='\n')
graphics::text(xxx[i],yyy[j],v2,cex=outlier_cex)
}
}
}
return(t1)
}
#' Set Color Scale for Z Statistics Value
#'
#' \code{z2col} is a helper function in \code{out2excel}. It defines the color scale of the Z statistics value.
#'
#' @param x a vector of numerics, a vector of Z statistics.
#' @param n_len integer, number of unique colors. Default is 60.
#' @param sig_thre numeric, the threshold for significance (absolute value of Z statistics). Z values failed to pass the threshold will be colored "white".
#' @param col_min_thre numeric, the lower threshold for the color bar value. Default is 0.01.
#' @param col_max_thre numeric, the upper threshold for the color bar value. Default is 3.
#' @param blue_col a vector of characters, the blue colors used to show the negative Z values. Default is brewer.pal(9,'Set1')[2].
#' @param red_col a vector of characters, the red colors used to show positive Z values. Default is brewer.pal(9,'Set1')[1].
#' @return Return a vector of color codes.
#' @examples
#' t1 <- sort(rnorm(mean=0,sd=2,n=100))
#' graphics::image(as.matrix(t1),col=z2col(t1))
#' @export
z2col <- function(x,n_len=60,sig_thre=0.01,col_min_thre=0.01,col_max_thre=3,
blue_col=brewer.pal(9,'Set1')[2],
red_col=brewer.pal(9,'Set1')[1]){
#
tmp_x <- setdiff(x,c(Inf,-Inf))
if(length(tmp_x)==0) return(ifelse(x>0,'red','blue'))
all_input_para <- c('x','n_len','sig_thre','col_min_thre','col_max_thre','blue_col','red_col')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
## create vector for z-score, can change sig threshold
x[which(is.na(x)==TRUE)] <- 0
x[which(x==Inf)]<- base::max(x[which(x!=Inf)])+1
x[which(x==-Inf)]<- base::min(x[which(x!=-Inf)])-1
if(col_min_thre<0) col_min_thre<-0.01
if(col_max_thre<0) col_max_thre<-3
c2 <- grDevices::colorRampPalette(c(blue_col,'white',red_col))(n_len)
r1 <- 1.05*base::max(abs(x)) ## -r1~r1
if(r1 < col_max_thre){
r1 <- col_max_thre
}
if(col_min_thre>r1){
r2 <- seq(-r1,r1,length.out=n_len+1)
}else{
r21 <- seq(-r1,-col_min_thre,length.out=n_len/2)
r22 <- base::seq(col_min_thre,r1,length.out=n_len/2)
r2 <- c(r21,r22)
}
x1 <- cut(x,r2)
names(c2) <- levels(x1)
x2 <- c2[x1]
x2[which(abs(x)<sig_thre)] <- 'white'
x2
}
#' Create Color Codes for a Vector of Characters
#'
#' \code{get.class.color} creates a vector of color codes for the input character vector. This is a helper function to assign nice looking colors for better visualization.
#'
#' @param x a vector of characters, names or labels.
#' @param use_color a vector of color codes, colors to be assigned to each member of \code{x}. Default is brewer.pal(9, 'Set1').
#' @param pre_define a vector of characters, pre-defined color codes for a certain input (e.g. c("blue", "red") with names c("A", "B")). Default is NULL.
#'
#' @return Return a vector of color codes, with input character vector as names.
#' @examples
#' get.class.color(c('ClassA','ClassB','ClassC','ClassA','ClassC','ClassC'))
#' get.class.color(c('ClassA','ClassB','ClassC','SHH','WNT','Group3','Group4'))
#' get.class.color(c('ClassA','ClassB','ClassC','SHH','WNT','Group3','Group4'),
#' use_color=brewer.pal(8, 'Set1'))
#'
#' pre_define <- c('blue', 'red', 'yellow', 'green','yellow', 'green')
#' ## pre-defined colors for MB
#' names(pre_define) <- c('WNT', 'SHH', 'Group3', 'Group4','GroupC', 'GroupD')
#' ##pre-defined color name for MB
#' get.class.color(c('ClassA','ClassB','ClassC','SHH','WNT','Group3','Group4'),
#' pre_define=pre_define)
#'
#' \dontrun{
#'}
#' @export
get.class.color <- function(x,use_color=NULL,pre_define=NULL) {
#
all_input_para <- c('x')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(min(check_res)==0){message('Please check and re-try!');return(FALSE)}
#
x <- clean_charVector(x);
#
if(is.null(pre_define)==FALSE & is.null(names(pre_define))==TRUE){
message('No class name for the color vector, please check and re-try !');return(FALSE);
}
x <- clean_charVector(x);
if(is.null(use_color)==TRUE){
use_color <- brewer.pal(9, 'Set1')
}
if (base::length(base::intersect(x, names(pre_define))) == 0) {
w1 <- base::length(base::unique(x))
if(w1 < length(use_color)){
cc2 <- use_color[1:w1]
}else{
cc2 <- grDevices::colorRampPalette(use_color)(base::length(base::unique(x)))
}
names(cc2) <- base::unique(x)
cc2 <- cc2[x]
} else{
x1 <- base::unique(x)
x2 <- base::setdiff(x1, names(pre_define))
cc1 <- NULL
w1 <- base::length(x2)
if (w1 > 0) {
if(w1 < length(use_color)){
cc1 <- use_color[1:w1]
}else{
cc1 <- grDevices::colorRampPalette(use_color)(w1)
}
names(cc1) <- x2
}
cc2 <- c(pre_define, cc1)
cc2 <- cc2[x]
}
return(cc2)
}
## get color box text,inner function ## refer from web
# https://stackoverflow.com/questions/45366243/text-labels-with-background-colour-in-r
boxtext <- function(x, y, labels = NA, col.text = NULL, col.bg = NA,
border.bg = NA, adj = NULL, pos = NULL, offset = 0.5,
padding = c(0.5, 0.5), cex = 1, font = graphics::par('font')){
## The Character expansion factro to be used:
theCex <- graphics::par('cex')*cex
## Is y provided:
if (missing(y)) y <- x
## Recycle coords if necessary:
if (base::length(x) != base::length(y)){
lx <- base::length(x)
ly <- base::length(y)
if (lx > ly){
y <- rep(y, ceiling(lx/ly))[1:lx]
} else {
x <- rep(x, ceiling(ly/lx))[1:ly]
}
}
## Width and height of text
textHeight <- graphics::strheight(labels, cex = theCex, font = font)
textWidth <- graphics::strwidth(labels, cex = theCex, font = font)
## Width of one character:
charWidth <- graphics::strwidth("e", cex = theCex, font = font)
## Is 'adj' of length 1 or 2?
if (!is.null(adj)){
if (base::length(adj == 1)){
adj <- c(adj[1], 0.5)
}
} else {
adj <- c(0.5, 0.5)
}
## Is 'pos' specified?
if (!is.null(pos)){
if (pos == 1){
adj <- c(0.5, 1)
offsetVec <- c(0, -offset*charWidth)
} else if (pos == 2){
adj <- c(1, 0.5)
offsetVec <- c(-offset*charWidth, 0)
} else if (pos == 3){
adj <- c(0.5, 0)
offsetVec <- c(0, offset*charWidth)
} else if (pos == 4){
adj <- c(0, 0.5)
offsetVec <- c(offset*charWidth, 0)
} else {
stop('Invalid argument pos')
}
} else {
offsetVec <- c(0, 0)
}
## Padding for boxes:
if (base::length(padding) == 1){
padding <- c(padding[1], padding[1])
}
## Midpoints for text:
xMid <- x + (-adj[1] + 1/2)*textWidth + offsetVec[1]
yMid <- y + (-adj[2] + 1/2)*textHeight + offsetVec[2]
## Draw rectangles:
rectWidth <- textWidth + 2*padding[1]*charWidth
rectHeight <- textHeight + 2*padding[2]*charWidth
graphics::rect(xleft = xMid - rectWidth/2,ybottom = yMid - rectHeight/2,
xright = xMid + rectWidth/2,ytop = yMid + rectHeight/2,
col = col.bg, border = border.bg,xpd=TRUE)
## Place the text:
graphics::text(xMid, yMid, labels, col = col.text, cex = theCex, font = font,adj = c(0.5, 0.5),xpd=TRUE)
## Return value:
if (base::length(xMid) == 1){
invisible(c(xMid - rectWidth/2, xMid + rectWidth/2, yMid - rectHeight/2,yMid + rectHeight/2))
} else {
invisible(base::cbind(xMid - rectWidth/2, xMid + rectWidth/2, yMid - rectHeight/2,yMid + rectHeight/2))
}
}
#' Visualize Sample Clustering Result in 2D Plot
#'
#' \code{draw.2D} creats a 2D plot to visualize the sample clustering result.
#'
#' @param X a vector of numerics, the x coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the first component.
#' @param Y a vector of numerics, the y coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the second component.
#' @param class_label a vector of characters, labels or categories of samples. The vector name should be sample names.
#' @param xlab character, the label for x-axis. Default is "PC1".
#' @param ylab character, the label for y-axis. Default is "PC2".
#' @param legend_cex numeric, giving the amount by which the text of legend should be magnified relative to the default. Default is 0.8.
#' @param main character, an overall title for the plot. Default is "".
#' @param point_cex numeric, giving the amount by which the size of the data points should be magnified relative to the default. Default is 1.
#' @param use_color a vector of color codes, colors to be assigned to each member of display label. Default is brewer.pal(9, 'Set1').
#' @param pre_define a vector of characters, pre-defined color codes for a certain input (e.g. c("blue", "red") with names c("A", "B")). Default is NULL.
#'
#' @return Return a logical value. If TRUE, the plot has been created successfully.
#' @examples
#' mat1 <- matrix(rnorm(2000,mean=0,sd=1),nrow=100,ncol=20)
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' colnames(mat1) <- paste0('Sample',1:ncol(mat1))
#' pc <- stats::prcomp(t(mat1))$x
#' pred_label <- kmeans(pc,centers=4)$cluster ## this can use other cluster results
#' draw.2D(X=pc[,1],Y=pc[,2],class_label=pred_label)
#' @export
draw.2D <- function(X,Y,class_label,xlab='PC1',ylab='PC2',legend_cex=0.8,main="",point_cex=1,use_color=NULL,pre_define=NULL){
#
all_input_para <- c('X','Y','class_label','xlab','ylab','legend_cex','main','point_cex')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(min(check_res)==0){message('Please check and re-try!');return(FALSE)}
class_label <- clean_charVector(class_label)
#
if(base::length(X)!=base::length(Y)){
message('Input two dimension vector with different length, please check and re-try !');return(FALSE);
}
if(base::length(X)!=base::length(class_label)){
message('Input dimension vector has different length with the class_label, please check and re-try !');return(FALSE);
}
par(mai = c(1, 1, 1, 0.5+base::max(strwidthMod(class_label,units='inch',cex=legend_cex,ori=FALSE,mod=FALSE))))
cls_cc <- get.class.color(class_label,use_color=use_color,pre_define=pre_define) ## get color for each label
graphics::plot(Y ~ X,pch = 16,cex = point_cex,col = cls_cc,main=main,xlab=xlab,ylab=ylab)
graphics::legend(par()$usr[2],par()$usr[4],sort(base::unique(class_label)),fill = cls_cc[sort(base::unique(class_label))],
horiz = FALSE,xpd = TRUE,border = NA,bty = 'n',cex=legend_cex)
return(TRUE)
}
#' Visualize Sample Clustering Result in 2D Plot with interactive mode
#'
#' \code{draw.2D.interactive} creats a 2D plot to visualize the sample clustering result with interactive mode realized by plotly.
#'
#' @param X a vector of numerics, the x coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the first component.
#' @param Y a vector of numerics, the y coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the second component.
#' @param sample_label a vector of characters, name of samples to be displayed on the figure.
#' @param color_label a vector of characters, labels used to define the point color.
#' @param shape_label a vector of characters, labels used to define the point shape.
#' @param xlab character, the label for x-axis. Default is "PC1".
#' @param ylab character, the label for y-axis. Default is "PC2".
#' @param main character, an overall title for the plot. Default is "".
#' @param point_cex numeric, giving the amount by which the size of the data points should be magnified relative to the default. Default is 1.
#' @param use_color a vector of color codes, colors to be assigned to each member of display label. Default is brewer.pal(9, 'Set1').
#' @param pre_define a vector of characters, pre-defined color codes for a certain input (e.g. c("blue", "red") with names c("A", "B")). Default is NULL.
#'
#' @return Return the plotly class object for interactive visualization.
#' @examples
#' mat1 <- matrix(rnorm(2000,mean=0,sd=1),nrow=100,ncol=20)
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' colnames(mat1) <- paste0('Sample',1:ncol(mat1))
#' pc <- stats::prcomp(t(mat1))$x
#' pred_label <- kmeans(pc,centers=4)$cluster ## this can use other cluster results
#' draw.2D.interactive(X=pc[,1],Y=pc[,2],
#' sample_label=rownames(pc),
#' color_label=pred_label,
#' pre_define = c('1'='blue','2'='red','3'='yellow','4'='green'))
#' draw.2D.interactive(X=pc[,1],Y=pc[,2],
#' sample_label=rownames(pc),
#' shape_label=pred_label)
#' @export
draw.2D.interactive <- function(X,Y,sample_label=NULL,color_label=NULL,shape_label=NULL,
xlab='PC1',ylab='PC2',main="",point_cex=1,
use_color=NULL,pre_define=NULL){
if(!'plotly' %in% rownames(installed.packages())){
message('plotly not installed!');return(FALSE);
}
#
all_input_para <- c('X','Y','sample_label','xlab','ylab','main','point_cex')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(min(check_res)==0){message('Please check and re-try!');return(FALSE)}
sample_label <- clean_charVector(sample_label)
if(is.null(color_label)==FALSE){
color_label <- clean_charVector(color_label)
}
if(is.null(shape_label)==FALSE){
shape_label <- clean_charVector(shape_label)
}
#
if(base::length(X)!=base::length(Y)){
message('Input two dimension vector with different length, please check and re-try !');return(FALSE);
}
if(base::length(X)!=base::length(sample_label)){
message('Input dimension vector has different length with the sample_label, please check and re-try !');return(FALSE);
}
if(is.null(shape_label)==TRUE & is.null(color_label)==TRUE){
message('Either color_label or shape_label is required!');return(FALSE);
}
if(is.null(shape_label)==FALSE & is.null(color_label)==FALSE){
if(base::length(X)!=base::length(color_label)){
message('Input dimension vector has different length with the color_label, please check and re-try !');return(FALSE);
}
color_label.factor <- as.factor(color_label)
cls_cc <- get.class.color(levels(color_label.factor),use_color=use_color,pre_define=pre_define) ## get color for each label
if(base::length(X)!=base::length(shape_label)){
message('Input dimension vector has different length with the shape_label, please check and re-try !');return(FALSE);
}
shape_label.factor <- as.factor(shape_label)
data <- data.frame(X=X,Y=Y,color_label=color_label.factor,shape_label=shape_label.factor);
display_text <- paste0(sample_label,':',color_label,':',shape_label)
p <- plotly::plot_ly(data = data, x = ~X, y = ~Y,
marker = list(size = point_cex*12),type='scatter',color=~color_label,colors=cls_cc,
hoverinfo='text',text=display_text,
mode='markers',symbol=~shape_label) %>%
plotly::layout(title = main,
xaxis = list(zeroline = FALSE,title=list(text=xlab)),#20240327
yaxis = list(zeroline = FALSE,title=list(text=ylab),
showlegend=TRUE)
)
}
if(is.null(shape_label)==TRUE & is.null(color_label)==FALSE){
if(base::length(X)!=base::length(color_label)){
message('Input dimension vector has different length with the color_label, please check and re-try !');return(FALSE);
}
color_label.factor <- as.factor(color_label)
cls_cc <- get.class.color(levels(color_label.factor),use_color=use_color,pre_define=pre_define) ## get color for each label
data <- data.frame(X=X,Y=Y,color_label=color_label.factor);
display_text <- paste0(sample_label,':',color_label)
p <- plotly::plot_ly(data = data, x = ~X, y = ~Y,
marker = list(size = point_cex*12),type='scatter',color=~color_label,colors=cls_cc,
hoverinfo='text',text=display_text,
mode='markers') %>%
plotly::layout(title = main,
yaxis = list(zeroline = FALSE,title=list(text=xlab)),
xaxis = list(zeroline = FALSE,title=list(text=ylab),showlegend=TRUE)
)
}
if(is.null(shape_label)==FALSE & is.null(color_label)==TRUE){
if(base::length(X)!=base::length(shape_label)){
message('Input dimension vector has different length with the shape_label, please check and re-try !');return(FALSE);
}
shape_label.factor <- as.factor(shape_label)
data <- data.frame(X=X,Y=Y,shape_label=shape_label.factor);
display_text <- paste0(sample_label,':',shape_label)
p <- plotly::plot_ly(data = data, x = ~X, y = ~Y,
marker = list(size = point_cex*12),type='scatter',color = I('black'),
hoverinfo='text',text=display_text,
mode='markers',symbol=~shape_label) %>%
plotly::layout(title = main,
yaxis = list(zeroline = FALSE,title=list(text=xlab)),
xaxis = list(zeroline = FALSE,title=list(text=ylab))
)
}
return(p)
}
#' Visualize Sample Clustering Result in 2D Plot with Sample Names
#'
#' \code{draw.2D.text} creates a 2D plot with sample names labeled, to visualize the sample clustering result.
#'
#' @param X a vector of numerics, the x coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the first component.
#' @param Y a vector of numerics, the y coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the second component.
#' @param class_label a vector of characters, labels or categories of samples. The vector name should be sample names.
#' @param xlab character, the label for x-axis. Default is "PC1".
#' @param ylab character, the label for y-axis. Default is "PC2".
#' @param legend_cex numeric, giving the amount by which the text of legend should be magnified relative to the default. Default is 0.8.
#' @param main character, an overall title for the plot. Default is "".
#' @param point_cex numeric, giving the amount by which the size of the data points should be magnified relative to the default. Default is 1.
#' @param class_text a vector of characters, the user-defined sample names to label each data points in the plot.
#' If NULL, will use the names of \code{class_label}. Default is NULL.
#' @param text_cex numeric, giving the amount by which the text of \code{class_text} should be magnified relative to the default. Default is NULL.
#' @param use_color a vector of color codes, colors to be assigned to each member of display label. Default is brewer.pal(9, 'Set1').
#' @param pre_define a vector of characters, pre-defined color codes for a certain input (e.g. c("blue", "red") with names c("A", "B")). Default is NULL.
#'
#' @return Return a logical value. If TRUE, the plot has been created successfully.
#' @examples
#' mat1 <- matrix(rnorm(2000,mean=0,sd=1),nrow=100,ncol=20)
#' rownames(mat1) <- paste0('Gene',1:nrow(mat1))
#' colnames(mat1) <- paste0('Sample',1:ncol(mat1))
#' pc <- stats::prcomp(t(mat1))$x
#' pred_label <- kmeans(pc,centers=4)$cluster ## this can use other cluster results
#' draw.2D.text(X=pc[,1],Y=pc[,2],class_label=pred_label,
#' point_cex=5,text_cex=0.5)
#' @export
draw.2D.text <- function(X,Y,class_label,class_text=NULL,xlab='PC1',ylab='PC2',legend_cex=0.8,main="",
point_cex=1,text_cex=NULL,use_color=NULL,pre_define=NULL){
#
all_input_para <- c('X','Y','class_label','xlab','ylab','legend_cex','main','point_cex')
check_res <- sapply(all_input_para,function(x)check_para(x,envir=environment()))
if(min(check_res)==0){message('Please check and re-try!');return(FALSE)}
class_label <- clean_charVector(class_label)
#
if(base::length(X)!=base::length(Y)){
message('Input two dimension vector with different length, please check and re-try !');return(FALSE);
}
if(base::length(X)!=base::length(class_label)){
message('Input dimension vector has different length with the class_label, please check and re-try !');return(FALSE);
}
par(mai = c(1, 1, 1, 0.5+base::max(strwidthMod(class_label,units='inch',cex=legend_cex,ori=FALSE,mod=FALSE))))
if(is.null(class_text)==TRUE){
class_text <- names(class_label)
}
cc <- 10/base::length(class_label)
if(cc<0.05) cc<-0.05
if(cc>1) cc<-1
if(is.null(text_cex)==FALSE) cc <- text_cex
cls_cc <- get.class.color(class_label,use_color=use_color,pre_define=pre_define) ## get color for each label
graphics::plot(Y ~ X,pch = 16,cex = point_cex,col = cls_cc,main=main,xlab=xlab,ylab=ylab)
graphics::text(x=X,y=Y,labels=class_text,cex=cc,xpd=TRUE,adj=0.5)
#print(cc);print(str(nn))
graphics::legend(par()$usr[2],par()$usr[4],sort(base::unique(class_label)),fill = cls_cc[sort(base::unique(class_label))],
horiz = FALSE,xpd = TRUE,border = NA,bty = 'n',cex=legend_cex)
return(TRUE)
}
#' Visualize Sample Clustering Result in 3D Plot
#'
#' \code{draw.3D} creates a 3D plot to visualize the sample clustering result.
#' @param X a vector of numerics, the x coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the first component.
#' @param Y a vector of numerics, the y coordinates of points in the plot. If user would like to create a PCA biplot, this parameter should be the

pipeline_functions.R at commit 5defa45, under Apache-2.0 · at the source

Overview

Authors: Winson S Ho^1,2, Isha Mondal^1,2, Jingjing Liu³, Raymond Sun^1,2, Jiawei Huo⁴, Chao Gao^1,2, Oishika Das^1,2, Daren Tieu^1,2, Jingqi Sun⁴, Hanchen Lin⁴, Peng Zhang⁴, Jiyang Yu³, Rongze Olivia Lu^1,2

ORCID iDs: Jiawei Huo, Jiyang Yu

Institutions: University of California, San Francisco (United States); UCSF Helen Diller Family Comprehensive Cancer Center (United States); St. Jude Children's Research Hospital (United States); Northwestern University (United States); Robert H. Lurie Comprehensive Cancer Center of Northwestern University

Journal: The Journal of clinical investigation, volume 136, issue 13, article e196753

Dates: received 24 July 2025; accepted 15 April 2026; published online 23 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1172/jci196753 · PMID 42085538 · PMCID PMC13318113 · OpenAlex W7160296515

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), mouse (organism), other condition (population), cellular / molecular (subfield)

Methods: Statistics, fMRI & imaging

Keywords: Immunology, Oncology, Cancer immunotherapy, Cellular senescence, Phosphoprotein phosphatases

MeSH: Cellular Senescence*, Cerebellar Neoplasms*, Medulloblastoma*, Neoplasm Proteins*, Neoplasms, Experimental*, Protein Phosphatase 2*, Animals, CD8-Positive T-Lymphocytes, Cell Line, Tumor, Humans, Mice, Mice, Knockout, Piperazines, Protein Phosphatase 2C (* major topic)

Topic: Protein Tyrosine Phosphatases (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: NINDS NIH HHS (R01 NS126501)

Citations: not cited yet (Europe PMC); 46 references in the paper

Abstract

Medulloblastoma (MB) is the most common malignant pediatric brain tumor. Current therapies are associated with substantial morbidity, and prognosis remains poor in high-risk subgroups, particularly those with TP53 mutations or relapsed disease. Cellular senescence is a tumor-suppressive program implicated in MB, but its role in antitumor immunity remains incompletely understood. We found that protein phosphatase 2A (PP2A) regulated immunogenic senescence in MB. Genetic ablation of the PP2A catalytic subunit PP2Ac or depletion of the regulatory subunit PP2A-B56α induced senescence in MB models. PP2Ac-deficient senescent cells exhibited increased MHC class I expression and enhanced immunogenicity. In syngeneic orthotopic models, PP2Ac loss prolonged survival in an immune- and CD8+ T cell–dependent manner. Analysis of patient datasets showed that senescence-associated gene signatures correlated with improved survival. Single-cell transcriptomic analysis further revealed that senescent MB cells were heterogeneous and that reduced PP2A activity was associated with an immunogenic senescence state. Because the PP2A inhibitor LB-100 has limited potency and off-target effects, we developed a lipid nanoparticle (LNP) platform to deliver siRNA targeting PPP2CA. LNP–small-interfering PP2Ac efficiently silenced PP2Ac in vitro and, when delivered locally in vivo, prolonged survival in a CD8+ T cell–dependent manner. Together, these findings identify PP2A as a regulator of immunogenic senescence in MB and support PP2Ac targeting as a therapeutic strategy.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

jyyulab/NetBID

License: Apache-2.0

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: 5defa454d600b94f5dd6d1f9f4428f99759a6821, 23 July 2026

Languages: R (7), JavaScript (3), Shell (2)

Size: 273 files, 12 scripts

Software Heritage: not archived

Found in: the text, “RNA-Seq.”

Holds: README, license file, environment (DESCRIPTION, Dockerfile, Dockerfile.cgc), continuous integration, documentation, 2 notebooks

Not found: CITATION.cff, tests

Tools: igraph (2 files), Plotly (2 files), ComplexHeatmap (1 file), DESeq2 (1 file), limma (1 file), lme4 (1 file), reshape2 (1 file), UMAP (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

14 files

IBC_CCDI/cloudAppNetBID.R, R, 132 lines
R/pipeline_functions.R, R, 3,700 lines, 1 match
build_docker_image.sh, Shell, 1 line
demo_scripts/analysis_and_plot_demo1.R, R, 413 lines
demo_scripts/pipeline_analysis_demo1.R, R, 170 lines
demo_scripts/pipeline_network_demo1.R, R, 151 lines
docs/_site/assets/js/dark-mode-preview.js, JavaScript, 23 lines
docs/_site/assets/js/just-the-docs.js, JavaScript, 295 lines
docs/_site/assets/js/vendor/lunr.min.js, JavaScript, 6 lines
inst/Rmd/eset_QC.Rmd, R, 173 lines
inst/Rmd/net_QC.Rmd, R, 199 lines
push_docker_image.sh, Shell, 1 line
LICENSE, License, 201 lines
README.md, Text, 74 lines

At the source: github.com/jyyulab/NetBID

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
12 scripts, each with its path and the digest of its content;
1 match between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

geo:GSE85217, at NCBI GEO; found in “Data availability”

Data availability

Publicly available datasets analyzed in this study include GEO GSE85217 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE85217) and GSE155446 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE155446). Sequencing data generated in this study are available in the GEO under accession GSE302307 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE302307). All other data generated in this study are provided in the Supporting Data Values file or are available from the corresponding author upon reasonable request. Supporting data values underlying all figures are provided in an XLS file.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 5 keywords, 14 MeSH terms, 1 funder, 46 references.

Cite

This paper

Ho, W. S., Mondal, I., Liu, J., Sun, R., Huo, J., Gao, C., Das, O., Tieu, D., Sun, J., Lin, H., Zhang, P., Yu, J., & Lu, R. O. (2026). Protein phosphatase 2A regulates senescence and immunogenicity in medulloblastoma models. The Journal of clinical investigation, 136(13), e196753. https://doi.org/10.1172/jci196753

BibTeX

@article{ho2026protein, author = {Ho, Winson S and Mondal, Isha and Liu, Jingjing and Sun, Raymond and Huo, Jiawei and Gao, Chao and Das, Oishika and Tieu, Daren and Sun, Jingqi and Lin, Hanchen and Zhang, Peng and Yu, Jiyang and Lu, Rongze Olivia}, title = {{Protein phosphatase 2A regulates senescence and immunogenicity in medulloblastoma models}}, journal = {The Journal of clinical investigation}, year = {2026}, month = apr, volume = {136}, number = {13}, pages = {e196753}, publisher = {American Society for Clinical Investigation}, issn = {0021-9738}, doi = {10.1172/jci196753}, url = {https://doi.org/10.1172/jci196753}, pmid = {42085538}, pmcid = {PMC13318113} }

RIS

TY - JOUR AU - Ho, Winson S AU - Mondal, Isha AU - Liu, Jingjing AU - Sun, Raymond AU - Huo, Jiawei AU - Gao, Chao AU - Das, Oishika AU - Tieu, Daren AU - Sun, Jingqi AU - Lin, Hanchen AU - Zhang, Peng AU - Yu, Jiyang AU - Lu, Rongze Olivia TI - Protein phosphatase 2A regulates senescence and immunogenicity in medulloblastoma models T2 - The Journal of clinical investigation J2 - J Clin Invest PY - 2026 DA - 2026/04/23 VL - 136 IS - 13 SP - e196753 SN - 0021-9738 PB - American Society for Clinical Investigation DO - 10.1172/jci196753 UR - https://doi.org/10.1172/jci196753 LA - en ER -

CSL-JSON

{ "id": "10.1172/jci196753", "type": "article-journal", "title": "Protein phosphatase 2A regulates senescence and immunogenicity in medulloblastoma models", "container-title": "The Journal of clinical investigation", "author": [ { "family": "Ho", "given": "Winson S" }, { "family": "Mondal", "given": "Isha" }, { "family": "Liu", "given": "Jingjing" }, { "family": "Sun", "given": "Raymond" }, { "family": "Huo", "given": "Jiawei" }, { "family": "Gao", "given": "Chao" }, { "family": "Das", "given": "Oishika" }, { "family": "Tieu", "given": "Daren" }, { "family": "Sun", "given": "Jingqi" }, { "family": "Lin", "given": "Hanchen" }, { "family": "Zhang", "given": "Peng" }, { "family": "Yu", "given": "Jiyang" }, { "family": "Lu", "given": "Rongze Olivia" } ], "container-title-short": "J Clin Invest", "volume": "136", "issue": "13", "page": "e196753", "DOI": "10.1172/jci196753", "PMID": "42085538", "PMCID": "PMC13318113", "ISSN": "0021-9738", "publisher": "American Society for Clinical Investigation", "URL": "https://doi.org/10.1172/jci196753", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 23 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.