Sparse regression with paired covariates

Code for reproducing the results shown in the manuscript.

Initialisation

### 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")
    
}

Preparation

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.html

Analysis

Last 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)

Figures

### 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)

}