The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
Code for reproducing the results shown in the manuscript.
### Loading functions. ###
inst <- rownames(utils::installed.packages())
cran <- c("devtools","R.utils","Matrix","glmnet","pROC","BiocManager","ashr")
# "googledrive", "httpuv"
if(!all(cran %in% inst)){
for(i in seq_along(cran)){
if(!cran[i] %in% inst){
install.packages(cran[i])
}
}
}
bioc <- c("edgeR","TCGAbiolinks")
if(!all(bioc %in% inst)){
#source("http://bioconductor.org/biocLite.R")
for(i in seq_along(bioc)){
if(!bioc[i] %in% inst){
#biocLite(bioc[i])
BiocManager::install(bioc[i])
}
}
}
#if(!"ashr" %in% inst){
# devtools::install_github("stephens999/ashr")
#}
user <- Sys.getenv("USERNAME")
path <- file.path("C:","Users",user,"Desktop","palasso")
if(user=="arra"){path <- "C:/Users/arra/Desktop/MATHS/palasso_desktop"}
if(user==""){path <- "/virdir/Scratch/arauschenberger/palasso"}
setwd(path)
folders <- c("data","results")
invisible(sapply(folders,function(x) if(!dir.exists(x)){dir.create(x)}))
if(user!="arra"){
devtools::install_github("kkdey/CorShrink") # ref="a9f6ba0"
devtools::install_github("rauschenberger/palasso") # ref="4a995a2"
}
if(FALSE){
# The functions <<save>>, <<file.exists>> and <<file.remove>> access the hard disk, but also try to access googledrive.
save <- function(object,file){
base::save(object,file=file)
tryCatch(expr=googledrive::drive_upload(media=file,path=file),
error=function(x) NULL)
#Sys.sleep(0.5)
}
file.exists <- function(file){
offline <- base::file.exists(file)
online <- FALSE
if(!offline){
d <- googledrive::as_dribble(x=file)
online <- tryCatch(expr=googledrive::some_files(d),
error=function(x) FALSE)
#Sys.sleep(0.5)
}
return(offline|online)
}
file.remove <- function(file){
base::file.remove(file)
tryCatch(expr=googledrive::drive_trash(file=file),
error=function(x) NULL)
#Sys.sleep(0.5)
}
# The function <<update>> moves files from the research cloud to the remote drive. Given both paths, it first verifies whether the folders SIM, GDC and CCA are available, and then copies all missing files from the research cloud to the remote drive.
# from: path to the origin
# to: path to the destination
update <- function(from,to){
dir <- c("SIM","GDC","CCA")
if(any(!dir.exists(file.path(from,dir)))){stop("Invalid.")}
if(any(!dir.exists(file.path(to,dir)))){stop("Invalid.")}
pb <- utils::txtProgressBar(min=0,max=1,style=3)
for(i in seq_along(dir)){
files0 <- dir(file.path(from,dir[i]))
files1 <- dir(file.path(to,dir[i]))
names <- files0[!files0 %in% files1]
for(j in seq_along(names)){
utils::setTxtProgressBar(pb=pb,value=(i-1)/3+(j*i)/(3*length(names)))
file.copy(from=file.path(from,dir[i],names[j]),
to=file.path(to,dir[i],names[j]),
copy.date=TRUE)
}
utils::setTxtProgressBar(pb=pb,value=i/3)
}
}
# update(from="results",to="//tsclient/N/palasso/results")
}Last data download started on 2019-03-26 (R version 3.5.3).
### Downloading "Isoform Expression Quantification". ###
#rm(list=ls())
#<<functions>>
directory <- file.path(path,"data")
setwd(directory)
# Retrieving cancer types:
project <- TCGAbiolinks::getGDCprojects()$id
project <- project[grepl(x=project,pattern="TCGA")]
# Downloading isoform expression data:
y <- X <- list()
for(i in seq_along(project)){
query <- TCGAbiolinks::GDCquery(project=project[i],
data.category="Transcriptome Profiling",
data.type="Isoform Expression Quantification")
TCGAbiolinks::GDCdownload(query,method="api",directory=directory)
trace(TCGAbiolinks:::readTranscriptomeProfiling,tracer=quote(ignore.case<-TRUE))
X[[i]] <- TCGAbiolinks::GDCprepare(query,directory=directory)
X[[i]][,c("miRNA_ID","reads_per_million_miRNA_mapped",
"cross-mapped","miRNA_region")] <- NULL
y[[i]] <- rep(project[i],times=length(unique(X[[i]]$barcode)))
}
save(list=c("y","X"),file=file.path(path,"data","isoform_raw.RData"))
load(file.path(path,"data","isoform_raw.RData"),verbose=TRUE)
# Merging isoform expression data:
Xs <- do.call(what=rbind,args=X) # sparse matrix
y <- do.call(what="c",args=y)
# Transform to matrix
Xs$isoform_coords <- gsub(pattern="hg38:chr",replacement="",x=Xs$isoform_coords)
samples <- unique(Xs$barcode)
covariates <- unique(Xs$isoform_coords)
row <- match(Xs$barcode,samples)
col <- match(Xs$isoform_coords,covariates)
X <- Matrix::sparseMatrix(i=row,j=col,x=Xs$read_count,dimnames=list(samples,covariates))
# Order by molecular location
split <- strsplit(x=colnames(X),split=":|-")
chr <- sapply(split,function(x) x[[1]])
pos <- sapply(split,function(x) x[[2]])
order <- order(chr,pos)
X <- X[,order]
if(FALSE){ # testing
i <- sample(seq_len(nrow(Xs)),size=1)
Xs$read_count[i]
X[Xs$barcode[i],Xs$isoform_coords[i]]
}
save(list=c("y","X"),file=file.path(path,"data","isoform_all.RData"))### Downloading "miRNA Expression Quantification". ###
#rm(list=ls())
#<<functions>>
directory <- file.path(path,"data")
setwd(directory)
# Downloading data.
project <- TCGAbiolinks::getGDCprojects()$id
project <- project[grepl(x=project,pattern="TCGA")]
y <- X <- list()
for(i in seq_along(project)){
query <- TCGAbiolinks::GDCquery(project=project[i],
data.category="Transcriptome Profiling",
data.type="miRNA Expression Quantification")
TCGAbiolinks::GDCdownload(query,method="api",directory=directory)
trace(TCGAbiolinks:::readTranscriptomeProfiling,tracer=quote(ignore.case<-TRUE))
data <- TCGAbiolinks::GDCprepare(query,directory=directory)
X[[i]] <- t(data[,c(seq(from=2,to=ncol(data),by=3))])
y[[i]] <- rep(project[i],times=nrow(X[[i]]))
}
save(list=c("y","X"),file=file.path(path,"data","miRNA_raw.RData"))
load(file.path(path,"data","miRNA_raw.RData"))
X <- do.call(what=rbind,args=X)
y <- do.call(what="c",args=y)
rownames(X) <- gsub(pattern="read_count_",replacement="",x=rownames(X))
save(list=c("y","X"),file=file.path(path,"data","miRNA_all.RData"))### Downloading "Gene Expression Quantification". ###
#rm(list=ls())
#<<functions>>
directory <- file.path(path,"data")
setwd(directory)
# Retrieving cancer types:
project <- TCGAbiolinks::getGDCprojects()$id
project <- project[grepl(x=project,pattern="TCGA")]
# Downloading data:
memory.limit(size=16000) # Activate virtual memory in system control!
y <- X <- list()
for(i in seq_along(project)){
query <- TCGAbiolinks::GDCquery(project=project[i],
data.category="Transcriptome Profiling",
data.type="Gene Expression Quantification",
workflow.type="HTSeq - Counts"); gc()
TCGAbiolinks::GDCdownload(query=query,method="api",directory=directory); gc()
trace(TCGAbiolinks:::readTranscriptomeProfiling,tracer=quote(ignore.case<-TRUE)); gc()
X[[i]] <- TCGAbiolinks::GDCprepare(query,directory=directory); gc()
y[[i]] <- rep(project[i],times=ncol(X[[i]])); gc()
}
save(list=c("y","X"),file=file.path(path,"data","gene_raw.RData"))
load(file.path(path,"data","gene_raw.RData"))
genes <- SummarizedExperiment::rowData(X[[1]])
mart <- biomaRt::useMart("ensembl",dataset="hsapiens_gene_ensembl") #
char <- biomaRt::getBM(attributes=c("ensembl_gene_id","chromosome_name","transcript_start","gene_biotype"),filters=c("biotype","chromosome_name"),values=list("protein_coding",c(1:22,"X")),mart=mart)
select <- genes$ensembl_gene_id[genes$ensembl_gene_id %in% char$ensembl_gene_id]
X <- lapply(X,function(x) t(SummarizedExperiment::assays(x)$"HTSeq - Counts"[select,]))
X <- do.call(what=rbind,args=X)
y <- do.call(what="c",args=y)
save(list=c("y","X"),file=file.path(path,"data","gene_all.RData"))### Downloading "Copy Number Variation". ###
#rm(list=ls())
#<<functions>>
directory <- file.path(path,"data")
setwd(directory)
project <- TCGAbiolinks::getGDCprojects()$id
project <- project[grepl(x=project,pattern="TCGA")]
y <- X <- list()
for(i in seq_along(project)){
query <- TCGAbiolinks::GDCquery(project=project[i],
data.category="Copy Number Variation",
data.type="Masked Copy Number Segment")
TCGAbiolinks::GDCdownload(query=query,method="api",directory=directory)
trace(TCGAbiolinks:::readTranscriptomeProfiling,tracer=quote(ignore.case<-TRUE))
X[[i]] <- TCGAbiolinks::GDCprepare(query,directory=directory)
y[[i]] <- rep(project[i],times=length(unique(X[[i]]$Sample))) # correct?
}
save(list=c("y","X"),file=file.path(path,"data","CNV_raw.RData"))
load(file.path(path,"data","CNV_raw.RData"),verbose=TRUE)
# Merging CNV data:
Xs <- do.call(what=rbind,args=X) # sparse matrix
y <- do.call(what="c",args=y)
#table(Xs$Sample)
# Prepare cutting.
cut <- list()
cut$chr <- c(1:22,"X")
cut$start <- sapply(cut$chr,function(x) min(Xs$Start[Xs$Chromosome==x]))
cut$end <- sapply(cut$chr,function(x) max(Xs$End[Xs$Chromosome==x]))
cut$length <- cut$end-cut$start
cut$dist <- sum(cut$length)/10000
cut$num <- round(cut$length/cut$dist)
# Create covariates.
cov <- list()
cov$p <- sum(cut$num)
cov$chromosome <- unlist(sapply(cut$chr,function(i) rep(i,times=cut$num[i])))
cov$location <- unlist(sapply(cut$chr,function(i) round(seq(from=cut$start[i],to=cut$end[i],length.out=cut$num[i]))))
cov$name <- paste0(cov$chromosome,":",cov$location)
# Create indices for each covariate.
index <- rep(list(integer()),times=cov$p)
pb <- utils::txtProgressBar(min=0,max=cov$p,style=3)
for(j in seq_len(cov$p)){
utils::setTxtProgressBar(pb=pb,value=j)
index[[j]] <- which((Xs$Chromosome==cov$chromosome[j]) &
(Xs$Start<=cov$location[j]) & (cov$location[j]<=Xs$End)) # consider <
}
# Expand indices to matrix.
X <- matrix(0,nrow=length(unique(Xs$Sample)),ncol=cov$p,
dimnames=list(unique(Xs$Sample),cov$name))
for(j in seq_along(index)){
mean <- Xs$Segment_Mean[index[[j]]]
i <- Xs$Sample[index[[j]]]
X[i,j] <- mean
}
if(FALSE){ # test
sample <- sample(rownames(X),size=1)
covariate <- sample(colnames(X),size=1)
split <- strsplit(covariate,split=":")[[1]]
a <- X[sample,covariate]
b <- Xs$Segment_Mean[(Xs$Sample==sample) & (Xs$Chromosome==split[1]) & (Xs$Start<=as.numeric(split[2])) & (as.numeric(split[2]) < Xs$End)]
all(a==b)
}
save(list=c("y","X","index"),file=file.path(path,"data","CNV_all.RData"))### Extracting samples of interest. ###
#rm(list=ls())
#<<functions>>
type <- c("isoform","miRNA","CNV","gene")
for(i in seq_along(type)){
cat(type[i],"\n")
load(file.path(path,"data",paste0(type[i],"_all.RData")),verbose=TRUE)
# TCGA barcode
barcode <- rownames(X)
code <- sapply(barcode,function(x) strsplit(x,split="-"))
code <- as.data.frame(do.call(what=rbind,args=code))
colnames(code) <- c("project","TSS","participant","sample_vial",
"portion_analyte","plate","center")
code$sample <- substr(code$sample_vial,start=1,stop=2)
code$vial <- substr(code$sample_vial,start=3,stop=3)
code$portion <- substr(code$portion_analyte,start=1,stop=2)
code$analyte <- substr(code$portion_analyte,start=3,stop=3)
code$sample_vial <- code$portion_analyte <- NULL
# solid tumour (except blood for LAML)
solid <- (code$sample=="01" | (y=="TCGA-LAML" & code$sample=="03"))
X <- X[solid,]
y <- y[solid]
# unique samples
unique <- !duplicated(substr(rownames(X),start=1,stop=12))
X <- X[unique,]
y <- y[unique]
save(list=c("y","X"),file=file.path(path,"data",paste0(type[i],"_sub.RData")))
}
# isoform: n=9'794, p=194'595, k=32
# miRNA: n=9'794, p=1'881, k=32
# gene: n=9'830, p=19'602, k=33
# CNV: n=10'578, p=10'000, k=33
## Understanding barcodes:
# overview: https://wiki.nci.nih.gov/display/TCGA/TCGA+barcode
# details: https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables
## Understanding replicate samples:
# http://gdac.broadinstitute.org/runs/sampleReports/latest/READ_Replicate_Samples.htmlLast data analysis started on 2019-04-12 (R version 3.5.3).
### Analysing the TCGA data. ###
#rm(list=ls())
#<<functions>>
for(type in c("miRNA","isoform","CNV","gene")){
library(Matrix)
for(k in c(207,sample(528))){
set.seed(k); cat(k," ")
if(type %in% c("isoform","miRNA") & k > 496){next}
if(type %in% c("CNV","gene") & k > 528){next}
# searching for missing cancer-cancer combinations
rm(list=setdiff(ls(),c("k","type","path","save","file.exists","file.remove"))); gc()
file0 <- paste0("results/",type,"_start_",k,".RData")
file1 <- paste0("results/",type,"_loss_",k,".RData")
if(file.exists(file0)||file.exists(file1)){next}
save(object=k,file=file0)
load(paste0("data/",type,"_sub.RData"))
# indicating the cancer-cancer combination
cancer <- substring(text=unique(y),first=6)
comb <- utils::combn(x=cancer,m=2)
select <- paste0("TCGA-",comb[,k])
y <- ifelse(y==select[1],1,ifelse(y==select[2],0,NA))
rm(cancer,select)
# removing other cancer types
cond <- !is.na(y)
y <- y[cond]
X <- X[cond,]
rm(cond)
# pre-processing
if(type %in% c("isoform","miRNA")){
x <- palasso:::.prepare(X,cutoff="zero")
} else if(type=="gene"){
x <- palasso:::.prepare(X,cutoff="knee")
} else if(type=="CNV"){
x <- list(X=X,Z=sign(X))
x <- lapply(x,function(x) scale(x))
attributes(x)$info <- data.frame(n=nrow(X),p=ncol(X),prop=mean(x$Z==0))
}
rm(X)
# cross-validation
loss <- tryCatch(expr=palasso:::.predict(y=y,X=x,nfolds.int=10),error=function(e) palasso:::.predict(y=y,X=x,nfolds.int=10))
# information
loss$info <- cbind(k=k,
y0=comb[2,k],
y1=comb[1,k],
n0=sum(y==0),
n1=sum(y==1),
attributes(x)$info,
loss$info)
# refit
object <- palasso::palasso(y=y,X=x,nfolds=10,family="binomial",standard=TRUE,elastic=TRUE,shrink=TRUE)
model <- c(names(object),"elastic",
paste0("paired.",c("adaptive","standard","combined")))
for(max in c(10,50,Inf)){
temp <- list()
temp$nzero <- data.frame(model=model,x=NA,z=NA)
for(i in seq_along(model)){
coef <- palasso:::coef.palasso(object=object,model=model[i],max=max)
temp$nzero$x[i] <- sum(coef$x!=0)
temp$nzero$z[i] <- sum(coef$z!=0)
}
temp$select <- palasso:::subset.palasso(x=object,max=max,
model="paired.adaptive")$palasso$select
temp$weights <- palasso:::weights.palasso(object=object,max=max,
model="paired.adaptive")
temp$coef <- palasso:::coef.palasso(object=object,max=max,
model="paired.adaptive")
loss[[paste0("fit",max)]] <- temp
}
save(object=loss,file=file1)
file.remove(file0)
}
index <- sum(grepl(dir(),pattern="sessionInfo"))
sink(paste0("sessionInfo",index+1,".txt"))
date()
utils::sessionInfo()
devtools::session_info()
sink()
}# The function <<collect>> loads all files from PATH including PATTERN in the file name, loads OBJECT into a list, and executes a function call.
#<<functions>>
# path: folder
# pattern: character, or NULL (all files)
# object: character vector, or NULL (all objects)
# what: function call
collect <- function(path=getwd(),pattern="",object=NULL,what="rbind"){
OBJECT = object
# identify files
files <- dir(path)
files <- files[grepl(x=files,pattern=pattern)]
files <- files[grepl(x=files,pattern=".RData")]
number <- gsub(pattern=paste0(pattern,"|.RData"),replacement="",x=files)
files <- files[order(as.numeric(number))] # trial1
names <- gsub(pattern=".RData",replacement="",x=files) # trial1
if(length(files)==0){stop("Invalid datasets.")}
# load data
all <- list()
for(i in seq_along(files)){
x <- load(file.path(path,files[i]))
x <- eval(parse(text=x))
if(is.null(OBJECT)){
all[[i]] <- x
names(all)[i] <- names[i]
} else {
for(j in seq_along(OBJECT)){
all[[OBJECT[j]]][[i]] <- x[[OBJECT[j]]]
names(all[[OBJECT[j]]])[i] <- names[i]
}
}
}
# fuse data
if(is.null(OBJECT)){
all <- do.call(what=what,args=all)
} else {
all <- lapply(all,function(x) do.call(what=what,args=x))
}
return(all)
}
LOSS <- list()
type <- c("gene","isoform","miRNA","CNV")
#type <- "miRNA"
for(i in seq_along(type)){
LOSS[[type[i]]] <- collect(path="results",
pattern=paste0(type[i],"_loss_"),
object=c("deviance","auc","class","info",
paste0("fit",c(10,50,Inf))))
}
for(i in seq_along(LOSS)){
for(j in 1:3){
colnames(LOSS[[i]][[j]])[colnames(LOSS[[i]][[j]])=="paired.adaptive"] <- "paired"
}
}
#type <- "gene"
#a <- LOSS[[type]]$deviance[rownames(LOSS[[type]]$deviance)=="10","paired"]
#b <- LOSS[[type]]$deviance[rownames(LOSS[[type]]$deviance)=="10","elastic"]
#mean(a<b)### Testing for significant differences. ###
#rm(list=ls())
#<<functions>>
#<<collect>>
row <- c("gene","isoform","miRNA","CNV")
col <- c("10","Inf")
lay <- c("standard_x","standard_z","standard_xz",
"adaptive_x","adaptive_z","adaptive_xz",
"elastic") # added "elastic"
M <- array(NA,dim=c(length(row),length(col),length(lay)),dimnames=list(row,col,lay))
for(i in seq_along(row)){
loss <- LOSS[[row[i]]][c("info","deviance")]
y0 <- as.character(loss$info$y0)
y1 <- as.character(loss$info$y1)
cancer <- sort(unique(c(y0,y1)))
Z <- palasso:::.design(x=cancer)
for(j in seq_along(col)){
# differences
cond <- rownames(loss$deviance)==col[j]
for(k in seq_along(lay)){
fill <- loss$deviance[cond,lay[k]] - loss$deviance[cond,"paired"]
X <- matrix(NA,nrow=length(cancer),ncol=length(cancer),
dimnames=list(cancer,cancer))
X[cbind(y0,y1)] <- X[cbind(y1,y0)] <- fill
X[lower.tri(X)] <- NA
# p-values
pvalue <- rep(NA,times=max(Z))
for(l in seq_len(max(Z))){
x <- as.numeric(X[Z==l])
if(col[j]=="10"){
alternative <- "greater" # Never use "two.sided"!
}
if(col[j]=="Inf"){
alternative <- "less" # Never use "two.sided"!
}
pvalue[l] <- stats::wilcox.test(x=x,alternative=alternative,
exact=FALSE)$p.value
}
# Simes
M[i,j,k] <- palasso:::.combine(pvalue,method="simes")
}
}
}
# Table SIG: significance
constraint <- "10"
table <- format(M[,constraint,1:6],digits=1,scientific=FALSE)
for(i in seq_len(nrow(table))){
for(j in seq_len(ncol(table))){
if(M[i,constraint,j]>=0.05){
table[i,j] <- paste0("{\\color{gray}{",table[i,j],"}}")
}
}
}
one <- c("","\\text{standard}","","","\\text{adaptive}","")
two <- paste0("\\text{",c("x","z","xz","x","z","xz"),"}")
rownames(table) <- paste0("\\text{",rownames(table),"}")
table <- rbind(one,two,table,deparse.level=0)
rownames(table)[1] <- "~"
xtable <- xtable::xtable(table,align=c("r","|","c","c","c","|","c","c","c"))
xtable::print.xtable(xtable,type="latex",include.colnames=FALSE,sanitize.text.function=identity)### Comparison with elastic net. ###
#rm(list=ls())
#<<functions>>
#<<collect>>
row <- c("gene","isoform","miRNA","CNV")
col <- c("10","50","Inf")
better <- worse <- less <- matrix(NA,nrow=length(row),ncol=length(col),
dimnames=list(row,col))
for(i in seq_along(row)){
for(j in seq_along(col)){
# proportion of improvements (cross-validation)
cond <- rownames(LOSS[[row[i]]]$deviance)==col[j]
loss <- LOSS[[row[i]]]$deviance[cond,c("paired","elastic")]
better[i,j] <- round(mean(loss[,"paired"]<loss[,"elastic"]),digits=2)
worse[i,j] <- round(mean(loss[,"paired"]>loss[,"elastic"]),digits=2)
# average difference in nzero (refitted models)
df_paired <- apply(LOSS[[row[i]]][[paste0("fit",col[j])]],1,function(x) sum(x$nzero[x$nzero[,"model"]=="paired.adaptive",c("x","z")]))
df_elastic <- apply(LOSS[[row[i]]][[paste0("fit",col[j])]],1,function(x) sum(x$nzero[x$nzero[,"model"]=="elastic95",c("x","z")]))
df_diff <- df_elastic-df_paired
less[i,j] <- round(mean(df_diff),digits=2)
#graphics::hist(df_diff,main=paste(row[i],col[j]),xlim=c(-1,1)*max(abs(df_diff)))
}
}
better
worse
less### Analysing the refitted models. ###
#rm(list=ls())
#<<functions>>
#<<collect>>
# Table SEL: selected model
nzero <- paste0("fit",c(5,10,Inf))
model <- c(paste0("standard_",c("x","z","xz")),
paste0("adaptive_",c("x","z","xz")),
"between_xz","within_xz")
type <- c("gene","isoform","miRNA","CNV")
table <- array(NA,dim=c(length(nzero),length(model),length(type)),
dimnames=list(nzero,model,type))
for(i in seq_along(nzero)){
for(j in seq_along(model)){
for(k in seq_along(type)){
sub <- LOSS[[type[k]]][[nzero[i]]]
table[i,j,k] <- sum(sub[,"select"]==model[j])
}
}
}
colSums(table["fit10",,]) # CHECK WHETHER COMPLETE!
table <- round(prop.table(table["fit10",,],margin=2),digits=2)
table <- t(table)
table <- table[,apply(table,2,function(x) any(x!=0))]
rownames(table) <- paste0("\\text{",rownames(table),"}")
xtable <- xtable::xtable(table,align=c("r","|","c","c","c","c"))
xtable::print.xtable(xtable,type="latex",include.colnames=FALSE,sanitize.text.function=identity)
# selected weights and covariates
type <- c("gene","isoform","miRNA","CNV")
group <- c("x","z")
model <- c(paste0("standard_",c("x","z","xz")),
paste0("adaptive_",c("x","z","xz")),
"paired.adaptive","elastic") # added "elastic" , paste0("elastic",c(100,75,50,25))
weights10 <- weightsInf <- matrix(NA,nrow=length(group),ncol=length(type),
dimnames=list(group,type))
coef10 <- coefInf <- array(NA,dim=c(length(group),length(type),length(model)),
dimnames=list(group,type,model))
for(i in seq_along(group)){
for(j in seq_along(type)){
weights10[,j] <- rowMeans(sapply(LOSS[[type[j]]]$fit10[,"weights"],colMeans))
weightsInf[,j] <- rowMeans(sapply(LOSS[[type[j]]]$fitInf[,"weights"],colMeans))
for(k in seq_along(model)){
coef10[i,j,k] <- mean(sapply(LOSS[[type[j]]]$fit10[,"nzero"], function(x) sum(x[x$model==model[k],group[i]])))
coefInf[i,j,k] <- mean(sapply(LOSS[[type[j]]]$fitInf[,"nzero"], function(x) sum(x[x$model==model[k],group[i]])))
}
}
}
# coef10["x",,]+coef10["z",,]
# with sparsity constraint
round(prop.table(weights10,margin=2),2)
round(prop.table(coef10[,,"paired.adaptive"],margin=2),2)
round(colSums(coef10[,,"paired.adaptive"]),2)
# natural sparsity
round(prop.table(weightsInf,margin=2),2)
round(prop.table(coefInf[,,"paired.adaptive"],margin=2),2)
round(colSums(coefInf[,,"paired.adaptive"]),2)
round(colSums(coefInf[,,"elastic"])/colSums(coefInf[,,"paired.adaptive"]),1) # multiple nzero of elastic net and paired lasso
# Table NSC: number of non-zero coefficients
table <- round(coefInf["x",,]+coefInf["z",,])
colnames(table)[7] <- "paired"
one <- c("","\\text{standard}","","","\\text{adaptive}","","\\text{paired}","\\text{elastic}")
two <- paste0("\\text{",c("x","z","xz","x","z","xz","xz","xz"),"}")
rownames(table) <- paste0("\\text{",rownames(table),"}")
table <- rbind(one,two,table,deparse.level=0)
rownames(table)[1] <- "~"
xtable <- xtable::xtable(table,align=c("r","|","c","c","c","|","c","c","c","|","c","|","c"))
xtable::print.xtable(xtable,type="latex",include.colnames=FALSE,sanitize.text.function=identity)### FIGURE CSW ###
#rm(list=ls())
#<<functions>>
set.seed(1)
overfit <- TRUE
# simulate
n <- 10
cx <- stats::rbeta(n=n,shape1=0.9,shape2=1)
cz <- stats::rbeta(n=n,shape1=0.4,shape2=0.9)
# collection
x <- list()
y <- list()
# adaptive weights (X only)
x[[1]] <- rep(1,times=n)
if(overfit){x[[1]] <- cx}
y[[1]] <- rep(0,times=n)
# adaptive weights (Z only)
x[[2]] <- rep(0,times=n)
y[[2]] <- rep(1,times=n)
if(overfit){y[[2]] <- cz}
# adaptive weights (X and Z)
x[[3]] <- y[[3]] <- rep(0.5,times=n)
if(overfit){x[[3]] <- cx}
if(overfit){y[[3]] <- cz}
# within-pair weights
x[[4]] <- cx^2/(cx+cz)
y[[4]] <- cz^2/(cx+cz)
# visualisation
graphics::par(mfrow=c(1,4),mar=c(4.5,0.5,0.5,0.5),oma=c(0,2,0,0))
for(i in seq_len(4)){
palasso:::plot_pairs(x=x[[i]],y=y[[i]],lwd=4)
if(i==1){
graphics::mtext(text="covariate pair",side=2,line=1)
}
}### FIGURE DIA ###
#rm(list=ls())
#<<functions>>
ellipse <- function(x,y,a=0.2,b=0.25,border=NA){
n <- max(c(length(x),length(y)))
if(length(x)==1){x <- rep(x,times=n)}
if(length(y)==1){y <- rep(y,times=n)}
if(length(border)==1){border <- rep(border,times=n)}
for(i in seq_len(n)){
angle <- seq(from=0,to=2*pi,length=100)
xs <- x[i] + a * cos(angle)
ys <- y[i] + b * sin(angle)
graphics::polygon(x=xs,y=ys,col=grey(0.9),border=border[i])
}
}
cancer <- c("ACC","BLCA","BRCA","UVM")
number <- c(80,409,1078,80)
col <- grDevices::rgb(red=0,green=0,blue=sample(seq(from=75,to=255,length.out=length(number))),maxColorValue=255)
lwd <- log(2*number)-2
lwd = pmax(5*number/max(number),1)
x1 <- 1.3; x2 <- 2; x3 <- 3
set.seed(1)
# first layer (bubble)
graphics::plot.new()
graphics::par(mar=c(0,0,0,0))
graphics::plot.window(xlim=c(1.1,3.3),ylim=c(-1.6,1.6))
ellipse(x=x1,y=0)
graphics::text(x=x1,y=0,labels="TCGA",font=2,adj=c(0.5,0))
graphics::text(x=x1,y=0,labels="n=9794",adj=c(0.5,1.2),cex=0.9)
# second layer (colon)
y1 <- seq(from=1,to=-1,length.out=length(cancer)+1)
graphics::text(x=x2,y=y1[length(cancer)],labels="...",font=2,srt=90)
y1 <- y1[-length(cancer)]
# first-second layer (connect)
graphics::segments(x0=x1+0.2,y0=0,x1=x2-0.2,y1=y1,lwd=lwd,col=col)
# second layer (bubble)
ellipse(x=x2,y=y1,a=0.2,b=0.2)
graphics::text(x=x2,y=y1,labels=cancer,font=2,col=col,adj=c(0.5,0))
graphics::text(x=x2,y=y1,labels=paste0("n=",number,""),adj=c(0.5,1.2),cex=0.9)
comb <- utils::combn(x=seq_along(y1),m=2)
# third layer (colon)
y2 <- seq(from=1.5,to=-1.5,length.out=ncol(comb)+1)
graphics::text(x=x3,y=y2[ncol(comb)],labels="...",font=2,srt=90)
y2 <- y2[-ncol(comb)]
# second-third layer (connect)
graphics::segments(x0=2.2,y0=y1[comb[1,]],x1=2.7,y1=y2,lwd=lwd[comb[1,]],col=col[comb[1,]])
graphics::segments(x0=2.2,y0=y1[comb[2,]],x1=2.7,y1=y2,lwd=lwd[comb[2,]],col=col[comb[2,]])
# third layer (bubble)
ellipse(x=x3,y=y2,a=0.3,b=0.22)
graphics::text(x=x3,y=y2,labels=paste0(cancer[comb[1,]]," "),
font=2,col=col[comb[1,]],adj=c(1,0))
graphics::text(x=x3,y=y2,labels=paste0(" ",cancer[comb[2,]]),
font=2,col=col[comb[2,]],adj=c(0,0))
graphics::text(x=x3,y=y2,labels=":",font=2,adj=c(0.5,0))
labels <- apply(comb,2,function(x) sum(number[x]))
labels <- paste0("n=",labels,"")
graphics::text(x=x3,y=y2,labels=labels,adj=c(0.5,1.2),cex=0.9)### FIGURE CLA ###
#rm(list=ls())
#<<functions>>
graphics::par(mfrow=c(4,2),oma=c(0,0,0,0),mar=c(2.1,3.5,0.5,0.5))
for(type in c("gene","isoform","miRNA","CNV")){
loss <- LOSS[[type]][c("deviance","auc","class")]
choice <- "paired"
loss <- lapply(loss,function(x) x[,c(paste0("standard_",c("x","z","xz")),paste0("adaptive_",c("x","z","xz")),choice)])
for(constraint in c("10")){ # c("5","10","Inf")
# change
sub <- lapply(loss,function(x) x[rownames(x)==constraint,])
palasso:::plot_score(sub$deviance,choice=choice)
change <- sub$deviance[,7]-sub$deviance[,-7]
palasso:::plot_box(change,ylab="change",zero=TRUE,choice=NA)
# info
info <- list()
info$select <- names(which.min(apply(sub$deviance,2,median)[-7]))
info$DEV_paired <- median(sub$deviance[,choice])
info$DEV_select <- median(sub$deviance[,info$select])
info$improve <- mean(sub$deviance[,info$select]>sub$deviance[,choice])
info$AUC_paired <- median(sub$auc[,choice])
info$CLASS_paired <- median(sub$class[,choice])
print(as.data.frame(info)) # important
}
}### FIGURE DEC ###
#rm(list=ls())
#<<functions>>
#graphics::par(oma=c(1.0,1.0,0,0),mar=c(1.5,3.0,0.5,0.5),mfrow=c(1,1))
graphics::par(oma=c(1.0,1.0,0,0),mar=c(1.5,3.0,0.5,0.5),mfrow=c(2,2))
for(type in c("gene","isoform","miRNA","CNV")){
models <- c(paste0("standard_",c("x","z","xz")),
paste0("adaptive_",c("x","z","xz")),"paired")
constraint <- c("3","4","5","10","15","20","25","50","Inf")
loss <- LOSS[[type]]["deviance"]
loss <- lapply(loss,function(x) x[,models])
table <- matrix(NA,nrow=length(constraint),ncol=length(models),
dimnames=list(constraint,models))
for(i in seq_along(constraint)){
sub <- lapply(loss,function(x) x[rownames(x)==constraint[i],])
table[i,] <- apply(sub$deviance,2,median)
}
# table <- log(table)
graphics::plot.new()
graphics::plot.window(xlim=c(1,length(constraint)),ylim=range(table))
graphics::box()
constraint[constraint=="Inf"] <- "n"
graphics::axis(side=2)
graphics::axis(side=1,at=seq_along(constraint),labels=constraint,tick=FALSE,line=-1)
for(k in c(1,2)){
for(i in seq_along(models)){
lty <- ifelse(i%in%c(1,2,3),3,ifelse(i%in%c(4,5,6),2,1))
col <- ifelse(i==7,"#00007F","#FF3535")
pch <- ifelse(i%in%c(1,4),"x",ifelse(i%in%c(2,5),"z",1))
if(k==1){
graphics::lines(table[,i],col=col,lty=lty,lwd=2)
graphics::points(table[,i],col="white",pch=16,cex=1.2)
} else {
graphics::points(table[,i],col=col,pch=1,font=2)
}
}
}
}
graphics::title(ylab="deviance",line=0.0,outer=TRUE)
graphics::title(xlab="sparsity constraint",ylab="deviance",line=0.0,outer=TRUE)### FIGURE CNV ###
#rm(list=ls())
#<<functions>>
graphics::par(oma=c(0,0,0,0),mar=c(2.1,3.5,0.5,0.5))
graphics::layout(matrix(c(1,1,2,2),nrow=1))
loss <- LOSS[[type]][c("deviance","auc","class")]
loss <- lapply(loss,function(x) x[rownames(x)=="10",])
model <- c(paste0("standard_",c("x","z","xz")),
paste0("adaptive_",c("x","z","xz")))
diff <- loss$auc[,"paired"]-loss$auc[,model]
palasso:::plot_box(diff,zero=TRUE,invert=TRUE,ylab="change")
diff <- loss$class[,"paired"]-loss$class[,model]
palasso:::plot_box(diff,zero=TRUE,ylab="change")### FIGURE MAP ###
#rm(list=ls())
#<<functions>>
loss <- LOSS[["CNV"]][c("info","auc")]
cancer <- sort(unique(c(levels(loss$info$y0),levels(loss$info$y1))))
X <- matrix(NA,nrow=length(cancer),ncol=length(cancer),dimnames=list(cancer,cancer))
#Z <- palasso:::.design(x=cancer)
y0 <- as.character(loss$info$y0)
y1 <- as.character(loss$info$y1)
X[cbind(y0,y1)] <- X[cbind(y1,y0)] <- loss$auc[rownames(loss$auc)=="10","paired"]
graphics::par(mar=c(0.5,3.0,3.0,0.5))
dimnames(X) <- lapply(dimnames(X),function(x) paste0(" ",x," "))
palasso:::plot_table(X=X,margin=-1,labels=FALSE,las=2,cex=0.7)
#sort(rowMeans(X,na.rm=TRUE),decreasing=TRUE)[1:2] # keep!### FIGURE COM ###
# 32 cancer types for isoform and miRNA
# 33 cancer types for gene and CNV
#rm(list=ls())
#<<functions>>
for(type in c("miRNA")){
cancer <- sort(unique(as.character(unlist(LOSS[[type]]$info[,c("y0","y1")]))))
n <- length(cancer)
z <- as.numeric(palasso:::.design(x=n))
x <- rep(seq_len(n),each=n)
y <- rep(seq(from=n,to=1,by=-1),times=n)
pch <- z
pch[pch==0] <- NA
pex <- c(".","O","*","+","o","-","'","x")
# colour
base <- grDevices::colorRampPalette(colors=c('darkblue','blue','red','darkred'))(n)
col <- rep(NA,times=length(z))
col[z==0] <- "white"
for(i in seq_len(n)){
col[z==i] <- base[i]
}
graphics::par(mfrow=c(1,1),mar=c(0,0,2,2))
graphics::plot.new()
graphics::plot.window(xlim=c(1,n),ylim=c(1,n))
graphics::points(x=x[pch<=25],y=y[pch<=25],
pch=pch[pch<=25],col=col[pch<=25],cex=0.9)
graphics::points(x=x[pch>25],y=y[pch>25],
pch=pex[(pch-25)[pch>25]],col=col[pch>25],cex=0.9)
graphics::segments(x0=0,x1=n+1,y0=n+1)
graphics::segments(x0=n+1,y0=n+1,y1=0)
graphics::segments(x0=0,x1=n+1,y0=n+1,y1=0,lty=2)
graphics::mtext(text=cancer,side=3,at=1:n,las=2,cex=0.7)
graphics::mtext(text=cancer,side=4,at=n:1,las=2,cex=0.7)
}These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.