R/complex_functions.R

In YAPSA: Yet Another Package for Signature Analysis

# Copyright © 2014-2019  The YAPSA package contributors
# This file is part of the YAPSA package. The YAPSA package is licenced under
# GPL-3

#' Construct a VRanges Object from a data frame
#'
#' In this package, big data frames are generated from cohort wide vcf-like
#' files. This function constructs a VRanges object from such a data frame by
#' using \code{\link{makeGRangesFromDataFrame}} from the package
#' \code{\link{GenomicRanges}}
#'
#' @param in_df A big dataframe constructed from a vcf-like file of a whole
#'   cohort. The first columns are those of a standard vcf file, followed by an
#'   arbitrary number of custom or user defined columns. One of these can carry
#'   a PID (patient or sample identifyier) and one can carry subgroup
#'   information.
#' @param in_keep.extra.columns in_seqinfo Argument passed on to
#'   \code{\link{makeGRangesFromDataFrame}}
#' @param in_seqinfo A seqInfo object, referring to the reference genome used.
#'   Argument passed on to \code{\link{makeGRangesFromDataFrame}}
#' @param in_seqnames.field Indicates the name of the column in which the
#'   chromosome is encoded
#' @param in_start.field Indicates the name of the column in which the start
#'   coordinate is encoded
#' @param in_end.field Indicates the name of the column in which the end
#'   coordinate is encoded
#' @param in_PID.field Indicates the name of the column in which the PID
#'   (patient or sample identifier) is encoded
#' @param in_subgroup.field Indicates the name of the column in which the
#'   subgroup information is encoded
#' @param in_strand.field Indicates the name of the column in which the
#'   strandedness is encoded
#' @param verbose_flag Verbose if 1
#'
#' @return The constructed VRanges object
#'
#' @examples
#'  data(lymphoma_test)
#'  temp_vr <- makeVRangesFromDataFrame(lymphoma_test_df,
#'                                      in_seqnames.field="CHROM",
#'                                      in_subgroup.field="SUBGROUP",
#'                                      verbose_flag=1)
#'
#' @seealso \code{\link{makeGRangesFromDataFrame}}
#'
#' @import GenomicRanges
#' @import VariantAnnotation
#' @importFrom GenomeInfoDb seqlengths seqlengths<-
#' @export
#' 
makeVRangesFromDataFrame <- function(in_df,in_keep.extra.columns=TRUE,
                                    in_seqinfo=NULL,
                                    in_seqnames.field="X.CHROM",
                                    in_start.field="POS",
                                    in_end.field="POS",in_PID.field="PID",
                                    in_subgroup.field="subgroup",
                                    in_strand.field="strand",
                                    verbose_flag=1) {
    out_vr <- NULL
    name_list <- tolower(names(in_df))
    match_list <- c("ref","alt",tolower(in_seqnames.field),
                    tolower(in_start.field),tolower(in_end.field))
    if(length(which(match_list %in% name_list)) == length(match_list)) {
        my_gr <- GenomicRanges::makeGRangesFromDataFrame(in_df,
                                keep.extra.columns=in_keep.extra.columns,
                                seqinfo=in_seqinfo,
                                seqnames.field=in_seqnames.field,
                                start.field=in_start.field,
                                end.field=in_end.field,
                                strand.field=in_strand.field)
        if(!("+" %in% unique(strand(my_gr))) & 
                !("-" %in% unique(strand(my_gr)))){
                strand(my_gr) <- "+"
                if(verbose_flag==1){
                    cat("YAPSA:::makeVRangesFromDataFrame::warning:strand ",
                    "information missing, set to \"+\".\n")
        }
    }    
    out_vr <- VRanges(seqnames=seqnames(my_gr),ranges=ranges(my_gr),
                        ref=my_gr$REF,alt=my_gr$ALT)
    if(tolower(in_PID.field) %in% name_list) {
        if(verbose_flag==1){
            cat("YAPSA:::makeVRangesFromDataFrame::in_PID.field found. ",
                "Retrieving PID information.\n")
        }
        column_ind <- 
            min(which(tolower(names(mcols(my_gr)))==tolower(in_PID.field)))
        out_vr$PID <- mcols(my_gr)[,column_ind]
        sampleNames(out_vr) <- mcols(my_gr)[,column_ind]
        } else if("pid" %in% name_list) {
            if(verbose_flag==1){
                cat(paste0(
                "YAPSA:::makeVRangesFromDataFrame::warning:in_PID.field ",
                "not a valid column name,",
                " but default is valid. Retrieving PID information.\n"))
        }
        out_vr$PID <- mcols(my_gr)[,"PID"]
        sampleNames(out_vr) <- mcols(my_gr)[,column_ind]
    } else {
        if(verbose_flag==1){
            cat("YAPSA:::makeVRangesFromDataFrame::warning:PID information ",
            "missing. Filling up with dummy entries.\n");}
            out_vr$PID <- "dummy_PID"
            sampleNames(out_vr) <- "dummyPID"
    }
    if(tolower(in_subgroup.field) %in% name_list) {
        if(verbose_flag==1){
            cat("YAPSA:::makeVRangesFromDataFrame::in_subgroup.field found. ",
                "Retrieving subgroup information.\n")
        }
    column_ind <- 
        min(which(tolower(names(mcols(my_gr)))==tolower(in_subgroup.field)))
        #out_vr$Type <- mcols(my_gr)[,column_ind]
        mcols(out_vr)[,in_subgroup.field] <- mcols(my_gr)[,column_ind]
    } else if("subgroup" %in% name_list) {
        if(verbose_flag==1){
            cat(paste0("YAPSA:::makeVRangesFromDataFrame::warning:",
                "in_subgroup.field not a valid column name, ",
                "but default is valid. Retrieving subgroup information.\n"))
    }
    #out_vr$Type <- my_gr$subgroup
    mcols(out_vr)[,in_subgroup.field] <- my_gr$subgroup
    } else {
        if(verbose_flag==1){
            cat("YAPSA:::makeVRangesFromDataFrame::warning:",
                "subgroup information missing. ",
                "Filling up with dummy entries.\n");}
        #out_vr$Type <- "dummy_subgroup"
        mcols(out_vr)[,in_subgroup.field] <- "dummy_subgroup"
    }
    seqlengths(out_vr) <- seqlengths(my_gr)    
    } else {
        if(verbose_flag==1){
            cat("YAPSA:::makeVRangesFromDataFrame::
                error:mismatch in column names, ","return NULL.\n")
        }
    }
    return(out_vr)
}


#' Create a Mutational Catalogue from a VRanges Object
#'
#' This function creates a mutational catalogue from a VRanges Object by first
#' calling \code{\link[SomaticSignatures]{mutationContext}} to establish the
#' motif context of the variants in the input VRanges and then calling
#' \code{\link[SomaticSignatures]{motifMatrix}} to build the mutational
#' catalogue \code{V}.
#'
#' @param in_vr A VRanges object constructed from a vcf-like file of a whole
#'   cohort. The first columns are those of a standard vcf file, followed by an
#'   arbitrary number of custom or used defined columns. One of these can carry
#'   a PID (patient or sample identifyier) and one can carry subgroup
#'   information.
#' @param in_refGenome The reference genome handed over to
#'   \code{\link[SomaticSignatures]{mutationContext}} and used to extract the
#'   motif context of the variants in \code{in_vr}.
#' @param in_wordLength The size of the motifs to be extracted by
#'   \code{\link[SomaticSignatures]{mutationContext}}
#' @param in_PID.field Indicates the name of the column in which the PID
#'   (patient or sample identifier) is encoded
#' @param in_verbose Verbose if \code{in_verbose=1}
#' @param in_rownames Optional parameter to specify rownames of the mutational
#'   catalogue \code{V} i.e. the names of the features.
#' @param adapt_rownames Rownames of the output \code{matrix} will be adapted if
#'   \code{adapt_rownames=1}
#'
#' @return A list with entries \code{matrix}, \code{frame}, \itemize{ \item
#'   \code{matrix}: The mutational catalogue \code{V} \item \code{frame}:
#'   Additional and meta information on rownames (features), colnames (PIDs) and
#'   subgroup attribution. }
#'
#' @examples
#'  library(BSgenome.Hsapiens.UCSC.hg19)
#'  data(lymphoma_test)
#'  data(sigs)
#'  word_length <- 3
#'  temp_vr <- makeVRangesFromDataFrame(
#'    lymphoma_test_df,in_seqnames.field="CHROM",
#'    in_subgroup.field="SUBGROUP",verbose_flag=1)
#'  temp_list <- create_mutation_catalogue_from_VR(
#'    temp_vr,in_refGenome=BSgenome.Hsapiens.UCSC.hg19,
#'    in_wordLength=word_length,in_PID.field="PID",
#'    in_verbose=1)
#'  dim(temp_list$matrix)
#'  head(temp_list$matrix)
#'  test_list <- split(lymphoma_test_df,f=lymphoma_test_df$PID)
#'  other_list <- list()
#'  for(i in seq_len(length(test_list))){
#'    other_list[[i]] <- test_list[[i]][c(1:80),]
#'  }
#'  other_df <- do.call(rbind,other_list)
#'  other_vr <- makeVRangesFromDataFrame(
#'    other_df,in_seqnames.field="CHROM",
#'    in_subgroup.field="SUBGROUP",verbose_flag=1)
#'  other_list <- create_mutation_catalogue_from_VR(
#'    other_vr,in_refGenome=BSgenome.Hsapiens.UCSC.hg19,
#'    in_wordLength=word_length,in_PID.field="PID",
#'    in_verbose=1,in_rownames=rownames(AlexCosmicValid_sig_df))
#'  dim(other_list$matrix)
#'  head(other_list$matrix)
#'
#' @seealso \code{\link[SomaticSignatures]{mutationContext}}
#' @seealso \code{\link[SomaticSignatures]{motifMatrix}}
#'
#' @importFrom SomaticSignatures mutationContext
#' @importFrom SomaticSignatures motifMatrix
#' @export
#' 
create_mutation_catalogue_from_VR <- function(in_vr,in_refGenome,in_wordLength,
                                            in_PID.field="PID",in_verbose=0,
                                            in_rownames=c(),adapt_rownames=1){
    # adapt levels of in_PID.field
    PID_index <- which(names(mcols(in_vr))==in_PID.field)
    names(mcols(in_vr))[PID_index] <- "PID"
    in_vr$PID <- factor(as.character(in_vr$PID))
    # now start real work of the function
    if ( in_verbose==1 ) {
        cat("YAPSA:::create_mutation_catalogue_from_VR::Extracting mutation ",
        "context for the SNVs from the reference genome...\n");
    }
    my_motifs <- mutationContext(in_vr, in_refGenome, 
                                    k=in_wordLength, unify=TRUE)
    if ( in_verbose==1 ) {
        cat("YAPSA:::create_mutation_catalogue_from_VR::Building mutational ",
        "catalogue...\n");
    }
    ## cave: mutationContextMatrix is the old version of the SomaticSignatures 
    ##       package motifMatrix is the new version, but return a list instead 
    ##       of a matrix when the data is sparse
    #my_matrix = mutationContextMatrix(my_motifs, group=in_PID.field,
    #                                  normalize = FALSE)
    my_matrix = motifMatrix(my_motifs, group=in_PID.field,normalize = FALSE)
    ## check if unallowed nucleotides in matrix
    deselect_ind <- grep("[NX]",rownames(my_matrix))
    if (length(deselect_ind) > 0) {
        my_matrix <- my_matrix[-deselect_ind]
    }
    # account for sparsity of the matrix
    if(length(in_rownames) > 0) {
        if(length(in_rownames) < dim(my_matrix)[1]) {
            if(in_verbose==1) {
                cat("YAPSA:::create_mutation_catalogue_from_VR::
                Warning: supplied ", "rownames too short\n");
            }
        } else if (!is.character(in_rownames)) {
            if(in_verbose==1) {
                cat("YAPSA:::create_mutation_catalogue_from_VR::
                    Warning: supplied ", "rownames not character format\n");
            }
        } else {
            my_rownames <- transform_rownames_MATLAB_to_R(in_rownames,
                                                        in_wordLength)
            temp_matrix <- matrix(0,nrow=length(my_rownames),
                                    ncol=dim(my_matrix)[2])
            match_ind <- match(rownames(my_matrix),my_rownames)
            if( (length(match_ind) < dim(my_matrix)[1]) |
                (any(is.na(match_ind)))){
                if(in_verbose==1) {
                    cat("YAPSA:::create_mutation_catalogue_from_VR::
                        Warning: my_matrix"," contains rownames which are not
                        in supplied rownames.\n");
                }
            } else {
                temp_colnames <- colnames(my_matrix)
                temp_matrix[match_ind,] <- my_matrix
                my_matrix <- temp_matrix
                rownames(my_matrix) <- my_rownames
                colnames(my_matrix) <- temp_colnames
            }
        }
    }
    # string manipulations for export to matlab
    if ( in_verbose==1 ) {
        cat("YAPSA:::create_mutation_catalogue_from_VR::
            Performing string ", "manipulations for output to matlab...\n");
    }
    temp_frame <- data.frame(strsplit(rownames(my_matrix)," "))
    colnames(temp_frame) <- rownames(my_matrix)
    my_frame <- as.data.frame(t(temp_frame[,]))
    colnames(my_frame) <- c("exchange","context")
    my_frame$types <- 
        paste0(substr(my_frame$exchange,1,1),">",substr(my_frame$exchange,2,2))
    ## adapt for word length here
    fraction_length <- floor(in_wordLength/2)
    my_frame$subtypes <- 
        paste0(substr(my_frame$context,1,fraction_length),
                substr(my_frame$exchange,1,1),
                substr(my_frame$context,
                in_wordLength-fraction_length+1,in_wordLength))  
    my_frame$newnames <- paste0(my_frame$types," ",my_frame$subtypes)
    if(adapt_rownames==1) {
        rownames(my_matrix) <- my_frame$newnames
    }
    return(list(matrix=my_matrix,frame=my_frame))
}


#' Create a Mutational Catalogue from a data frame
#'
#' This function creates a mutational catalogue from a data frame. It is a
#' wrapper function for \code{\link{create_mutation_catalogue_from_VR}}: it
#' first creates a VRanges object from the data frame by
#' \code{\link{makeVRangesFromDataFrame}} and then passes this object on to the
#' above mentioned custom function.
#'
#' @param this_df A data frame constructed from a vcf-like file of a whole
#'   cohort. The first columns are those of a standard vcf file, followed by an
#'   arbitrary number of custom or used defined columns. One of these can carry
#'   a PID (patient or sample identifyier) and one can carry subgroup
#'   information.
#' @param this_refGenome_Seqinfo A seqInfo object, referring to the reference
#'   genome used. Argument passed on to \code{\link{makeGRangesFromDataFrame}}
#'   and thus indirectly to \code{\link{makeGRangesFromDataFrame}}.
#' @param this_seqnames.field Indicates the name of the column in which the
#'   chromosome is encoded
#' @param this_start.field Indicates the name of the column in which the start
#'   coordinate is encoded
#' @param this_end.field Indicates the name of the column in which the end
#'   coordinate is encoded
#' @param this_PID.field Indicates the name of the column in which the PID
#'   (patient or sample identifier) is encoded
#' @param this_subgroup.field Indicates the name of the column in which the
#'   subgroup information is encoded
#' @param this_refGenome The reference genome handed over to
#'   \code{\link{create_mutation_catalogue_from_VR}} and indirectly to
#'   \code{\link[SomaticSignatures]{mutationContext}} and used to extract the
#'   motif context of the variants in \code{in_vr}.
#' @param this_wordLength The size of the motifs to be extracted by
#'   \code{\link[SomaticSignatures]{mutationContext}}
#' @param this_verbose Verbose if \code{this_verbose=1}
#' @param this_rownames Optional parameter to specify rownames of the mutational
#'   catalogue \code{V} i.e. the names of the features.
#' @param this_adapt_rownames Rownames of the output \code{matrix} will be
#'   adapted if \code{this_adapt_rownames=1}
#'
#' @return A list with entries \code{matrix} and \code{frame} obtained from
#'   \code{\link{create_mutation_catalogue_from_VR}}: \itemize{ \item
#'   \code{matrix}: The mutational catalogue \code{V} \item \code{frame}:
#'   Additional and meta information on rownames (features), colnames (PIDs) and
#'   subgroup attribution. }
#'
#' @examples
#'  library(BSgenome.Hsapiens.UCSC.hg19)
#'  data(lymphoma_test)
#'  word_length <- 3
#'  temp_list <- create_mutation_catalogue_from_df(
#'    lymphoma_test_df,this_seqnames.field = "CHROM",
#'    this_start.field = "POS",this_end.field = "POS",
#'    this_PID.field = "PID",this_subgroup.field = "SUBGROUP",
#'    this_refGenome = BSgenome.Hsapiens.UCSC.hg19,
#'    this_wordLength = word_length)
#'    dim(temp_list$matrix)
#'    head(temp_list$matrix)
#'
#' @seealso \code{\link{makeVRangesFromDataFrame}}
#' @seealso \code{\link{create_mutation_catalogue_from_VR}}
#'
#' @export
#' 
create_mutation_catalogue_from_df <- function(this_df,
                                        this_refGenome_Seqinfo=NULL,
                                        this_seqnames.field="X.CHROM",
                                        this_start.field="POS",
                                        this_end.field="POS",
                                        this_PID.field="PID",
                                        this_subgroup.field="subgroup",
                                        this_refGenome,this_wordLength,
                                        this_verbose=1,
                                        this_rownames=c(),
                                        this_adapt_rownames=1) {
    if(is.null(this_refGenome_Seqinfo)){
        this_refGenome_Seqinfo <- seqinfo(this_refGenome)
    }
    my_vr <- makeVRangesFromDataFrame(this_df,
                                    in_keep.extra.columns=TRUE,
                                    in_seqinfo=this_refGenome_Seqinfo,
                                    in_seqnames.field=this_seqnames.field,
                                    in_start.field=this_start.field,
                                    in_end.field=this_end.field,
                                    in_PID.field=this_PID.field,
                                    in_subgroup.field=this_subgroup.field,
                                    verbose_flag=this_verbose)
    my_result_list <- create_mutation_catalogue_from_VR(my_vr,
                                                this_refGenome,
                                                this_wordLength,
                                                this_PID.field,
                                                this_verbose,
                                                this_rownames,
                                                this_adapt_rownames)
    return(my_result_list)
}


#' Change rownames from one naming convention to another
#'
#' Rownames or names of the features used differ between the different contexts
#' a signature analysis is carried out in. The function
#' \code{transform_rownames_R_to_MATLAB} changes from the convention used in
#' the YAPSA pacakge to the one used by Alexandrov et al. in the MATLAB
#' framework.
#'
#' @param in_rownames 
#'  Character vector of input rownames
#' @param wordLength 
#'  Size of the considered motif context
#'
#' @return A character vector of the translated rownames.
#' 
#' @examples
#' NULL
#' 
#' @export
#'
transform_rownames_R_to_MATLAB <- function(in_rownames,wordLength=3) {
    out_rownames <- in_rownames
    temp_frame <- data.frame(strsplit(in_rownames," "))
    colnames(temp_frame) <- in_rownames
    my_frame <- as.data.frame(t(temp_frame[,]))
    colnames(my_frame) <- c("exchange","context")
    my_frame$types <- paste0(substr(my_frame$exchange,1,1),">",
                                substr(my_frame$exchange,2,2))
    ## adapt for word length here
    fraction_length <- floor(wordLength/2)
    my_frame$subtypes <- paste0(
        substr(my_frame$context,1,fraction_length),
        substr(my_frame$exchange,1,1),
        substr(my_frame$context,wordLength-fraction_length+1,wordLength))  
    out_rownames <- paste0(my_frame$types," ",my_frame$subtypes)
    return(out_rownames)
}


#' Change rownames from one naming convention to another
#'
#' The function \code{transform_rownames_MATLAB_to_R} changes from the 
#' convention used in Alexandrov et al. in the MATLAB framework to the one used
#' by the YAPSA pacakge.
#' 
#' @export
#' @rdname transform_rownames_R_to_MATLAB
#'
transform_rownames_MATLAB_to_R <- function(in_rownames,wordLength=3) {
    out_rownames <- in_rownames
    temp_frame <- data.frame(strsplit(in_rownames," "))
    colnames(temp_frame) <- in_rownames
    my_frame <- as.data.frame(t(temp_frame[,]))
    colnames(my_frame) <- c("exchange","context")
    my_frame$types <- paste0(substr(my_frame$exchange,1,1),
                                substr(my_frame$exchange,3,3))
    ## adapt for word length here
    fraction_length <- floor(wordLength/2)
    my_frame$subtypes <- paste0(
        substr(my_frame$context,1,fraction_length),".",
        substr(my_frame$context,wordLength-fraction_length+1,wordLength))  
    out_rownames <- paste0(my_frame$types," ",my_frame$subtypes)
    return(out_rownames)
}


#' Change rownames from one naming convention to another
#'
#' The function \code{transform_rownames_MATLAB_to_R} changes from the 
#' convention used in stored mutational catalogues by Alexandrov et al. to the 
#' one used by the YAPSA pacakge.
#'
#' @export
#' @rdname transform_rownames_R_to_MATLAB
#'
transform_rownames_nature_to_R <- function(in_rownames,wordLength=3) {
    out_rownames <- in_rownames
    ## adapt for word length here
    fraction_length <- floor(wordLength/2)
    string_length <- nchar(as.character(in_rownames[1]))
    my_frame <- repeat_df(in_value = "a",in_rows = length(in_rownames),
                            in_cols = 2)
    colnames(my_frame) <- c("types","subtypes")
    my_frame$subtypes <- paste0(
        substr(in_rownames,1,fraction_length),
        substr(in_rownames,fraction_length+2,fraction_length+2),
        substr(in_rownames,string_length-fraction_length+1,string_length))  
    my_frame$types <- substr(in_rownames,fraction_length+2,fraction_length+4)
    out_rownames <- paste0(my_frame$types," ",my_frame$subtypes)
    return(out_rownames)
}


#' Change rownames from one naming convention to another
#'
#' The function \code{transform_rownames_YAPSA_to_deconstructSigs} changes from
#' the convention used in the YAPSA package to the one used by the 
#' deconstructSigs package.
#'
#' @export
#' @rdname transform_rownames_R_to_MATLAB
#'
transform_rownames_YAPSA_to_deconstructSigs <- function(in_rownames,
                                                        wordLength=3) {
    out_rownames <- in_rownames
    ## adapt for word length here
    fraction_length <- floor(wordLength/2)
    string_length <- nchar(as.character(in_rownames[1]))
    my_frame <- repeat_df(in_value = "a",in_rows = length(in_rownames),
                            in_cols = 2)
    colnames(my_frame) <- c("types","subtypes")
    my_frame$subtypes <- lapply(in_rownames,function(x) strsplit(x," ")[[1]][2])
    my_frame$types <- lapply(in_rownames,function(x) strsplit(x," ")[[1]][1])
    out_rownames <- paste0(
        substr(my_frame$subtypes,1,fraction_length),"[",
        my_frame$types,"]",
        substr(my_frame$subtypes,fraction_length+2,2*fraction_length+1))
    return(out_rownames)
}


#' Change rownames from one naming convention to another
#'
#' The function \code{transform_rownames_YAPSA_to_deconstructSigs} changes from
#' the convention used in the deconstructSigs package to the one used by the 
#' YAPSA pacakge.
#'
#' @export
#' @rdname transform_rownames_R_to_MATLAB
#'
transform_rownames_deconstructSigs_to_YAPSA <- function(in_rownames,
                                                        wordLength=3) {
    out_rownames <- in_rownames
    ## adapt for word length here
    fraction_length <- floor(wordLength/2)
    string_length <- nchar(as.character(in_rownames[1]))
    my_frame <- repeat_df(in_value = "a",in_rows = length(in_rownames),
                            in_cols = 2)
    colnames(my_frame) <- c("types","subtypes")
    my_frame$subtypes <- paste0(
        substr(in_rownames,1,fraction_length),
        substr(in_rownames,2+fraction_length,2+fraction_length),
        substr(in_rownames,6+fraction_length,5+2*fraction_length))
    my_frame$types <- substr(in_rownames,2+fraction_length,4+fraction_length)
    out_rownames <- paste0(my_frame$types," ",my_frame$subtypes)  
    return(out_rownames)
}


#' Normalize Somatic Motifs with different rownames
#'
#' This is a wrapper function to
#' \code{\link[SomaticSignatures]{normalizeMotifs}}. The rownames are first
#' transformed to fit the convention of the
#' \code{\link[SomaticSignatures]{SomaticSignatures}} package and then passed on
#' to the above mentioned function.
#'
#' @param in_matrix,in_norms Arguments to
#'   \code{\link[SomaticSignatures]{normalizeMotifs}}
#' @param adjust_counts Whether to rescale the counts after adaption or not.
#'   Default is true.
#'
#' @return The matrix returned by
#' \code{\link[SomaticSignatures]{normalizeMotifs}}, but with rownames
#' transformed back to the convention of the input
#'
#' @examples
#'  NULL
#'
#' @importFrom SomaticSignatures normalizeMotifs
#' @export
#' 
normalizeMotifs_otherRownames <- function(in_matrix,
                                            in_norms,
                                            adjust_counts=TRUE) {
    my_matrix <- in_matrix
    rownames(my_matrix) <- transform_rownames_MATLAB_to_R(rownames(in_matrix))
    out_matrix <- SomaticSignatures::normalizeMotifs(my_matrix,in_norms)
    rownames(out_matrix) <- rownames(in_matrix)
    if(adjust_counts){
        total_counts <- colSums(in_matrix)
        intermediate_counts <- colSums(out_matrix)
        mult_factors <- total_counts/intermediate_counts
        out_matrix <- as.data.frame(t(t(out_matrix)*mult_factors))
    }
    return(out_matrix)
}


adjust_number_of_columns_in_list_of_catalogues <- 
    function(in_all_list, in_merged_results_list) {
    out_list <- list()
    number_of_strata <- length(in_all_list$results_lists_list)
    reference_matrix <- in_merged_results_list$matrix
    zero_df <- as.data.frame(reference_matrix-reference_matrix)
    colname_vector <- colnames(reference_matrix)
    for (i in seq_len(number_of_strata)) {
        temp_df <- as.data.frame(in_all_list$results_lists_list[[i]]$matrix)
        match_ind <- match(names(temp_df),colname_vector)
        new_df <- zero_df
        new_df[,match_ind] <- temp_df
        out_list[[i]] <- new_df
    }
    return(out_list)
}


save_mutation_catalogue <- function(in_table,in_matrix,in_frame,
                                    in_subgroup_flag,in_meta_info,in_subDir,
                                    in_verbose) {
    type_df <- aggregate(data=in_table,subgroup~PID, function(l) return(l[1]))
    if (dim(type_df)[1] != length(colnames(in_matrix))) {
        cat("YAPSA:::save_mutation_catalogue::Warning: 
            mismatch between subgroup",
            " information and colnames(in_matrix).\n");
    }
    subgroup_vector <- rep("dummy",dim(in_matrix)[2])
    if (in_subgroup_flag==1) {
        subgroup_vector <- 
            as.character(in_meta_info$Diagnosis[match(colnames(in_matrix),
                                                in_meta_info$PID)])
    }
    if (in_subgroup_flag==2) {
        subgroup_vector <- 
            as.character(type_df$subgroup[match(colnames(in_matrix),
                                                type_df$PID)])
    }
    # store to text files
    if ( in_verbose==1 ) {cat("Store to specified files...\n");}
    save(in_table,in_matrix,file=file.path(in_subDir,"MotifMatrix.RData"))
    write.table(in_matrix,sep="\t",file=file.path(in_subDir,"MotifMatrix.txt"),
                col.names=FALSE,row.names=FALSE)
    write(c("PIDs",colnames(in_matrix)),sep="\t",
        file=file.path(in_subDir,"MotifMatrix_colnames.txt"))
    write(c("subgroup",subgroup_vector),sep="\t",
        file=file.path(in_subDir,"MotifMatrix_subgroups.txt"))
    write(c("types",in_frame$types),sep="\t",
        file=file.path(in_subDir,"MotifMatrix_types.txt"))
    write(c("subtypes",in_frame$subtypes),sep="\t",
        file=file.path(in_subDir,"MotifMatrix_subtypes.txt"))
}


stratify_vcf_like_df <- function(in_table,in_column_name,in_verbose=0) {
    out_table_list <- list()
    out_name_list <- list()
    stratification_index <- which(names(in_table)==in_column_name)
    name_list <- unique(in_table[,stratification_index])
    index_range <- sort(name_list)
    counter <- 0
    if(in_verbose == 1) {
        cat("\nYAPSA:::stratify_vcf_like_df::Stratify...\n")
    }
    for (i in index_range) {
        counter <- counter + 1
        temp_ind <- which(in_table[,stratification_index]==i)
        temp_table <- in_table[temp_ind,]
        temp_name <- i
        if(in_verbose == 1) {
            cat("YAPSA:::stratify_vcf_like_df:: temp_name =",temp_name,
                "; number of variants:",dim(temp_table)[1],"; number of PIDs:",
                length(unique(temp_table$PID)),"\n")
        }
        out_table_list[[counter]] <- temp_table
        out_name_list[[counter]] <- temp_name
    }    
    return(list(table_list=out_table_list,name_list=out_name_list))
}


save_stratified_vcf_like_df <- function(in_table_list,in_name_list,
                                        in_target_dir,in_verbose) {
    number_of_strata <- length(in_table_list)
    subDir <- file.path(in_target_dir, "vcf_like")
    strataDir <- file.path(in_target_dir, "strata")
    dir.create(in_target_dir)
    dir.create(subDir)
    dir.create(strataDir)
    if(in_verbose == 1) {
        cat("\nYAPSA:::save_stratified_vcf_like_df::Write to file...\n")
    }
    for (i in seq_len(number_of_strata)) {
        temp_table <- in_table_list[[i]]
        temp_name <- in_name_list[[i]]
        if(in_verbose == 1) {
            cat("YAPSA:::save_stratified_vcf_like_df:: 
                temp_name =",temp_name,"; number of variants:",
                dim(temp_table)[1],"; number of PIDs:",
                length(unique(temp_table$PID)),"\n")
        }
        file_name <- paste0(subDir,"/somatic_SNVs_",temp_name,".tsv")
        names(temp_table) <- sub("X.CHROM","\\#CHROM",names(temp_table))
        write.table(temp_table, file=file_name,sep="\t",
                    quote =FALSE, row.names=FALSE)
        temp_dir <- file.path(in_target_dir, "strata", temp_name)
        dir.create(temp_dir)
    }    
    return(1)
}


stratify_and_create_mutational_catalogue <- 
    function(our_table, our_column_name, vcf_target_dir, strata_target_dir,
            our_refGenome_Seqinfo, our_seqnames.field="X.CHROM",
            our_start.field="POS", our_end.field="POS", our_PID.field="PID",
            our_subgroup.field="subgroup", our_refGenome, our_wordLength,
            in_verbose=1, our_rownames=c()) {
    results_lists_list <- list()
    results_list <- stratify_vcf_like_df(our_table,our_column_name,in_verbose)
    table_list <- results_list$table_list
    name_list <- results_list$name_list
    if(!is.null(vcf_target_dir)) {
        status <- save_stratified_vcf_like_df(table_list,name_list,
                                            vcf_target_dir,in_verbose)    
    }
    number_of_strata <- length(table_list)
    if(in_verbose == 1) {
        cat("\nYAPSA:::stratify_and_create_mutational_catalogue::Create ",
            "mutational catalogues...\n")
    }
    for (i in seq_len(number_of_strata)) {
        if(in_verbose == 1) {
            cat("\nYAPSA:::stratify_and_create_mutational_catalogue::
                Processing ","stratum: ",i,"\n")
    }
        temp_table <- table_list[[i]]
        temp_name <- name_list[[i]]
        temp_results_list <- create_mutation_catalogue_from_df(
            temp_table,
            our_refGenome_Seqinfo,
            this_seqnames.field=our_seqnames.field,
            this_start.field=our_start.field,
            this_end.field=our_end.field,
            this_PID.field=our_PID.field,
            this_subgroup.field=our_subgroup.field,
            our_refGenome,our_wordLength,in_verbose,
            our_rownames)
        if(!is.null(strata_target_dir)) {
            temp_matrix <- temp_results_list$matrix
            temp_frame <- temp_results_list$frame
            MotifMatrix_dir <- file.path(strata_target_dir,temp_name)
            save_mutation_catalogue(temp_table,temp_matrix,temp_frame,
                                    2,NULL,MotifMatrix_dir,in_verbose)
        }
        results_lists_list[[i]] <- temp_results_list
    }    
    return(list(table_list=table_list,name_list=name_list,
                results_lists_list=results_lists_list))  
}


#' Wrapper for cut
#'
#' In this wrapper function for the known \code{\link{cut}} function, the
#' \code{breaks} vector need not be supplied directly, instead, for every break,
#' an interval is supplied and the function optimizes the choice of the
#' breakpoint by chosing a local minimum of the distribution.
#'
#' @param in_vector Vector of numerical continuously distributed input
#' @param in_outlier_cutoffs Interval specifyinf the upper and lower bounds of
#'   the range to be considered
#' @param in_cutoff_ranges_list List if intervals in which the cutoffs for
#'   \code{\link{cut}} have to be optimized.
#' @param in_labels Labels assigned to the strata or factors returned
#' @param in_name String specifying the name of the quantity analyzed (and
#'   plotted on the x-axis of the figure to be created).
#' @param output_path Path where the figure produced by the density function
#'   should be stored if non-NULL.
#'
#' @return A list with entries \code{category_vector}, and \code{density_plot}
#'   and \code{cutoffs} \itemize{ \item \code{category_vector}: Factor vector of
#'   the categories or strata, of the same length as \code{in_vector} \item
#'   \code{density_plot}: Density plot produced by the density function and
#'   indication of the chosen cutoffs. \item \code{cutoffs}: Vector of the
#'   computed optimal cutoffs }
#'
#' @examples
#'  data(lymphoma_test)
#'  lymphoma_test_df$random_norm <- rnorm(dim(lymphoma_test_df)[1])
#'  temp_list <- cut_breaks_as_intervals(
#'    lymphoma_test_df$random_norm,
#'    in_outlier_cutoffs=c(-4,4),
#'    in_cutoff_ranges_list=list(c(-2.5,-1.5),c(0.5,1.5)),
#'    in_labels=c("small","intermediate","big"))
#'   temp_list$density_plot
#'
#' @seealso \code{\link{cut}}
#' @seealso \code{\link[stats]{density}}
#'
#' @export
#' 
cut_breaks_as_intervals <- function(in_vector,
                                    in_outlier_cutoffs=c(0,3000),
                                    in_cutoff_ranges_list=list(c(60,69),
                                                                c(25,32)),
                                    in_labels=c("late","intermediate","early"),
                                    in_name="",
                                    output_path=NULL) {
    ## prepare
    .e <- environment()
    finite_ind <- which(is.finite(in_vector))
    finite_vector <- in_vector[finite_ind]
    choice_ind <- which(finite_vector<=in_outlier_cutoffs[2] &
                        finite_vector>=in_outlier_cutoffs[1])
    my_density <- density(finite_vector[choice_ind])
    my_density_df <- data.frame(x=my_density$x,y=my_density$y)
    ## find appropriate cutoffs in the ranges provided (refine!!)
    number_of_cutoffs <- length(in_cutoff_ranges_list)
    my_cutoff_vector <- rep(0,number_of_cutoffs)
    for(j in seq_len(number_of_cutoffs)) {
        my_cutoff_interval <- in_cutoff_ranges_list[[j]]
        my_range_ind <- which(my_density_df$x<my_cutoff_interval[2] & 
                                my_density_df$x>my_cutoff_interval[1])
        my_min_ind <- which.min(my_density_df$y[my_range_ind])
        my_cutoff_vector[j] <- my_density_df$x[my_range_ind[my_min_ind]]
    }
    ## make density plot
    p <- ggplot(my_density_df,aes_string(x="x",y="y"), environment = .e) + 
        geom_line() + labs(x=in_name, y="density") + 
        geom_vline(xintercept=my_cutoff_vector,colour="red")
    if(!is.null(output_path)) {
        png(output_path,width=400,height=400); p; dev.off()    
    }
    my_cutoff_vector <- c(in_outlier_cutoffs[1],my_cutoff_vector,
                            in_outlier_cutoffs[2])
    ## discretize into categories
    cut_vector <- cut(in_vector,my_cutoff_vector,right=FALSE,labels=in_labels)
    cut_vector <- as.character(cut_vector)
    cut_vector[which(is.na(cut_vector))] <- "undetermined"
    cut_vector <- factor(cut_vector)
    return(list(category_vector=cut_vector,density_plot=p,
                cutoffs=my_cutoff_vector))
}


#' Return those PIDs which have an extreme pattern for signature exposure
#'
#' For all signatures found in a project, this function returns the sample
#' identifiers (PIDs) with extremely high or extremely low exposures of the
#' respective signatures.
#'
#' @param in_exposures_df Data frame with the signature exposures
#' @param in_quantile Quantile for the amount of extreme PIDs to be selected.
#'
#' @return A data frame with 4 rows per signature (high PIDs, high exposures,
#'   low PIDs, low exposures); the number of columns depends on the quantile
#'   chosen.
#'
#' @examples
#'  data(lymphoma_cohort_LCD_results)
#'  get_extreme_PIDs(lymphoma_Nature2013_COSMIC_cutoff_exposures_df,0.05)
#'
#' @import ggplot2
#' @export
#' 
get_extreme_PIDs <- function(in_exposures_df,in_quantile=0.03) {
    my_number_of_sigs <- dim(in_exposures_df)[1]
    my_number_of_PIDs <- dim(in_exposures_df)[2]
    my_number_of_extremes <- ceiling(in_quantile*my_number_of_PIDs)
    extreme_PID_df <- repeat_df(0,my_number_of_extremes,my_number_of_sigs*4)
    for(i in seq_len(my_number_of_sigs)) {
        this_sig <- rownames(in_exposures_df)[i]
        my_exposure_vector <- in_exposures_df[i,]
        my_breaks <- as.numeric(quantile(my_exposure_vector,
                        probs=c(in_quantile,(1-in_quantile))))
        extreme_low_ind <- which(my_exposure_vector<=my_breaks[1])
        extreme_high_ind <- which(my_exposure_vector>=my_breaks[2])
        extreme_PID_df[,4*i-3] <- names(
            my_exposure_vector[extreme_low_ind[seq_len
                                                (my_number_of_extremes)]])
        names(extreme_PID_df)[4*i-3] <- paste0(this_sig,"_low_PIDs")
        extreme_PID_df[,4*i-2] <- as.numeric(
            my_exposure_vector[extreme_low_ind[seq_len
                                                (my_number_of_extremes)]])
        names(extreme_PID_df)[4*i-2] <- paste0(this_sig,"_low_exp")
        extreme_PID_df[,4*i-1] <- names(
            my_exposure_vector[extreme_high_ind[seq_len
                                                (my_number_of_extremes)]])
        names(extreme_PID_df)[4*i-1] <- paste0(this_sig,"_high_PIDs")
        extreme_PID_df[,4*i] <- as.numeric(
            my_exposure_vector[extreme_high_ind[seq_len
                                                (my_number_of_extremes)]])
        names(extreme_PID_df)[4*i] <- paste0(this_sig,"_high_exp")
    }
    out_df <- data.frame(t(extreme_PID_df))
    names(out_df) <- paste0("PID_",seq_len(dim(out_df)[2]))
    return(out_df)
}


#' Test if mutated PIDs are enriched in signatures
#'
#' For all signatures found in a project, this function tests whether PIDs
#' having mutations in a specified list of genes of interest have significantly
#' higher exposures.
#'
#' @param in_gene_list List with genes of interest
#' @param in_exposure_df Data frame with the signature exposures
#' @param in_mut_table Data frame or table of mutations (derived from
#'   vcf-format)
#' @param in_gene.field Name of the column in which the gene names are to be
#'   looked up
#' @param in_p_cutoff Significance threshold
#'
#' @return A list with entries \code{pvals}, \code{exposure_df},
#'   \code{number_of_mutated}, \itemize{ \item \code{pvals}: p-values of the
#'   t-tests performed on mutated vs. unmutated PIDs \item \code{exposure_df}:
#'   Transposed input exposures data frame with additional annotations for
#'   mutation status \item \code{number_of_mutated}: Number of PIDs carrying a
#'   mutation }
#'
#' @examples
#'  NULL
#'
#' @export
#' 
test_gene_list_in_exposures <- function(in_gene_list,
                                        in_exposure_df,
                                        in_mut_table,
                                        in_gene.field="GENE_short",
                                        in_p_cutoff=0.05) {
    gene_index <- which(names(in_mut_table)==in_gene.field)
    number_of_signatures <- dim(in_exposure_df)[1]
    choice_ind <- which(in_mut_table[,gene_index] %in% in_gene_list)
    #choice_table <- in_mut_table[choice_ind,c(1:5,16:18,20,24,26,29)]
    choice_table <- in_mut_table[choice_ind,]
    this_PIDs <- unique(choice_table$PID)
    temp_ind <- match(this_PIDs,names(in_exposure_df))
    temp_ind_ind <- which(temp_ind>0)
    cat("test_gene_list_in_exposures: number of PIDs in overlap ",
        "between exposures, gene list and mutation table: ",
        length(temp_ind_ind),"\n")
    this_exposure_df <- as.data.frame(t(in_exposure_df))
    this_exposure_df$this_status <- "none"
    this_exposure_df$this_col <- "green"
    this_exposure_df$this_status[temp_ind[temp_ind_ind]] <- "mut"
    this_exposure_df$this_col[temp_ind[temp_ind_ind]] <- "red"
    ## now check for differences in signature counts
    significance_counter <- 0
    out_t_test_vector <- c()
    for (i in seq_len(number_of_signatures)) {
        temp_test_dat <- t.test(
            this_exposure_df[which(this_exposure_df$this_status=="none"),i],
            this_exposure_df[which(this_exposure_df$this_status=="mut"),i],
            alternative="less")
        out_t_test_vector[i] <- temp_test_dat$p.value
        if(out_t_test_vector[i] < in_p_cutoff & !is.na(out_t_test_vector[i])) {
            significance_counter <- significance_counter + 1
        }
    }
    cat("test_gene_list_in_exposures: Number of significancies: ",
        significance_counter,"\n")
    return(list(pvals=out_t_test_vector,
                exposure_df=this_exposure_df,
                number_of_mutated=length(temp_ind_ind)))
}


#' Compare one mutational catalogue to reference mutational catalogues
#'
#' Compare one mutational catalogue (e.g. of one index patient) to a list of
#' reference mutational catalogues (e.g. from the initial Alexandrov
#' puplication) by cosine similarities
#'
#' @param in_index_df Data frame containing the mutational catalogue of interest
#' @param in_comparison_list List of data frames (ideally named) containing the
#'   reference mutational catalogues
#'
#' @return A similarity dataframe
#'
#' @examples
#'  NULL
#'
#' @export
#' 
compare_to_catalogues <- function(in_index_df,in_comparison_list){
    sim_vec_list <- lapply(in_comparison_list,FUN=function(x){
        temp_list <- compare_sets(in_index_df,x)
        return(c(as.numeric(1-as.matrix(temp_list$distance))))
    })
    sim_df_list <- list()
    for(i in seq_along(sim_vec_list)) {
        my_entity <- names(sim_vec_list)[i]
        temp_df <- data.frame(entity=my_entity,sim=sim_vec_list[[i]])
        sim_df_list[[i]] <- temp_df
    }
    sim_df <- do.call(rbind,sim_df_list)
    return(sim_df)
}


#' Extract statistical measures for entity comparison
#'
#' Compare one mutational catalogue (e.g. of one index patient) to a list of
#' reference mutational catalogues (e.g. from the initial Alexandrov
#' puplication) by cosine similarities
#'
#' @param in_sim_df A similarity data frame as extracted by
#'   \code{\link{compare_to_catalogues}}
#'
#' @return A dataframe containing statistical measures, prepared for bar plot
#'
#' @examples
#'  NULL
#'
#' @export
#' 
compute_comparison_stat_df <- function(in_sim_df){
    mean_df <- aggregate(sim~entity,data=in_sim_df,FUN=mean)
    median_df <- aggregate(sim~entity,data=in_sim_df,FUN=median)
    sd_df <- aggregate(sim~entity,data=in_sim_df,FUN=sd)
    stat_df <- mean_df
    stat_df$median_sim <- median_df$sim
    stat_df$sd_sim <- sd_df$sim
    names(stat_df)[2] <- "mean_sim"
    return(stat_df)
}


#' Aggregate exposures by category
#'
#' If a valid category (i.e. it matches to a category specified in
#' in_sig_ind_df) is supplied, then the exposures are aggregated over this
#' category.
#'
#' @param in_exposures_df Input data frame of exposures.
#' @param in_sig_ind_df Input data frame of meta information on the signatures.
#'   Has to match the signatures in \code{in_exposures_df}
#' @param in_category Category to be aggregated over
#'
#' @return A list with entries: \itemize{ \item \code{exposures}: The exposures
#'   \code{H}, a numeric data frame with \code{l} rows and \code{m} columns,
#'   \code{l} being the number of aggregated signatures and \code{m} being the
#'   number of samples \item \code{norm_exposures}: The normalized exposures
#'   \code{H}, a numeric data frame with \code{l} rows and \code{m} columns,
#'   \code{l} being the number of aggregated signatures and \code{m} being the
#'   number of samples \item \code{out_sig_ind_df}: Data frame of the type
#'   \code{signature_indices_df}, i.e. indicating name, function and
#'   meta-information of the aggregated signatures.. }
#'
#' @seealso \code{\link{LCD_complex_cutoff}}
#'
#' @examples
#'  NULL
#'
#' @export
#' 
aggregate_exposures_by_category <- function(
    in_exposures_df,in_sig_ind_df,in_category){
    if(any(in_category==names(in_sig_ind_df))) {
        ## process exposures
        temp_exposures <- cbind(
            in_exposures_df,in_sig_ind_df[,in_category,drop=FALSE])
        names(temp_exposures) <- gsub(in_category,"cat",names(temp_exposures))
        temp_aggregate_exposures <- aggregate(.~cat,data=temp_exposures,FUN=sum)
        rownames(temp_aggregate_exposures) <- temp_aggregate_exposures$cat
        temp_aggregate_exposures$cat <- NULL
        ## normalize exposures
        norm_aggregate_exposures <- 
            normalize_df_per_dim(temp_aggregate_exposures, 2)
        ## process sig_ind_df
        names(in_sig_ind_df) <- gsub(in_category,"cat",names(in_sig_ind_df))
        out_sig_ind_df <- aggregate(.~cat,data=in_sig_ind_df,FUN=head,1)
        out_sig_ind_df$sig <- out_sig_ind_df$cat
        out_sig_ind_df$cat <- NULL
        out_sig_ind_df$colour <- as.character(out_sig_ind_df$colour)
        return(list(exposures=temp_aggregate_exposures,
                    norm_exposures=norm_aggregate_exposures,
                    out_sig_ind_df=out_sig_ind_df))
    } else {
        return(NULL)
    }
}