R/test_isoform_switches.R

In kvittingseerup/IsoformSwitchAnalyzeR: Identify, Annotate and Visualize Isoform Switches with Functional Consequences from both short- and long-read RNA-seq data

### Test via satuRn
isoformSwitchTestSatuRn <- function(
    ### Core arguments
    switchAnalyzeRlist,
    alpha = 0.05,
    dIFcutoff = 0.1,

    ### Advanced arguments
    reduceToSwitchingGenes = TRUE,
    reduceFurtherToGenesWithConsequencePotential = TRUE,
    onlySigIsoforms = FALSE,
    keepIsoformInAllConditions = TRUE,
    diagplots = TRUE,
    showProgress = TRUE,
    quiet = FALSE
){
    # preamble
    if (TRUE) {
        if (class(switchAnalyzeRlist) != "switchAnalyzeRlist") {
            stop(paste("The object supplied to 'switchAnalyzeRlist'",
                       "must be a 'switchAnalyzeRlist'"))
        }
        if (switchAnalyzeRlist$sourceId == "preDefinedSwitches") {
            stop(paste("The switchAnalyzeRlist is made from pre-defined isoform switches",
                       "which means it is made without defining conditions (as it should be).",
                       "\nThis also means it cannot be used to test for new isoform switches -",
                       "if that is your intention you need to create a new",
                       "switchAnalyzeRlist with the importRdata() function and try again.",
                       sep = " "))
        }
        if (!is.logical(reduceToSwitchingGenes)) {
            stop(paste("The argument supplied to 'reduceToSwitchingGenes'",
                       "must be an a logic"))
        }
        if (dIFcutoff < 0 | dIFcutoff > 1) {
            stop("The dIFcutoff must be in the interval [0,1].")
        }
        if (reduceToSwitchingGenes) {
            if (alpha < 0 | alpha > 1) {
                stop("The alpha parameter must be between 0 and 1 ([0,1]).")
            }
            if (alpha > 0.05) {
                warning(paste("Most journals and scientists consider an alpha",
                              "larger than 0.05 untrustworthy. We therefore recommend",
                              "using alpha values smaller than or equal to 0.05"))
            }
        }
        if (is.null(switchAnalyzeRlist$isoformCountMatrix)) {
            stop(paste("An isoform replicate count matrix is nessecary for using",
                       "the satuRn approach - please reinitalize the",
                       "swithcAnalyzeRlist object with one of the import*()",
                       "functions and try again."))
        }
        if (any(switchAnalyzeRlist$conditions$nrReplicates ==
                1)) {
            stop(paste("A statistical test cannot be performed without replicates.",
                       "Please remove all conditions with only 1 replicate and try again.",
                       "\nThis can either be done before creating the switchAnalyzeRlist",
                       "in the first place or with the subsetSwitchAnalyzeRlist() function",
                       sep = " "))
        }
        comparisonsToMake <- unique(switchAnalyzeRlist$isoformFeatures[, c(
            'condition_1', 'condition_2'
        )])
        nConditions <- length(unique(switchAnalyzeRlist$designMatrix$condition))
        oneWayData <- nConditions <= 2

        if(all( switchAnalyzeRlist$conditions$nrReplicates <= 5)) {
            warning(paste(
                'You seem to have few replicates. We therfore recomend you use isoformSwitchTestDEXSeq() instead.'
            ))
        }
    }

    # Create SummarizedExperiment
    if(TRUE) {
        if (!quiet) {
            message("Step 1 of 4: Creating SummarizedExperiment data object...")
        }
        localDesign <- switchAnalyzeRlist$designMatrix

        ### Convert group of interest to factors
        localDesign$condition <- factor(localDesign$condition, levels=unique(localDesign$condition))

        if( ncol(localDesign) > 2 ) {
            for(i in 3:ncol(localDesign) ) { # i <- 4
                if( class(localDesign[,i]) %in% c('numeric', 'integer') ) {
                    ### if there are at least two samples per "condition"
                    if( uniqueLength( localDesign[,i] ) * 2 <= length(localDesign[,i]) ) {
                        localDesign[,i] <- factor(localDesign[,i])
                    }
                }
            }
        }

        ### Make formula for model
        localFormula <- '~ 0 + condition'
        if (ncol(localDesign) > 2) {
            localFormula <- paste(
                localFormula,
                '+',
                paste(
                    colnames(localDesign)[3:ncol(localDesign)],
                    collapse = ' + '
                ),
                sep=' '
            )
        }
        localFormula <- as.formula(localFormula)

        ### Make model (design matrix)
        localModel <- model.matrix(localFormula, data = localDesign)
        indexToModify <- 1:nConditions
        colnames(localModel)[indexToModify] <- gsub(
            pattern =  '^condition',
            replacement =  '',
            x =  colnames(localModel)[indexToModify]
        )

        ### Subset count matrix (after prefiltering)
        localCount <- switchAnalyzeRlist$isoformCountMatrix[which(
            switchAnalyzeRlist$isoformCountMatrix$isoform_id %in%
                switchAnalyzeRlist$isoformFeatures$isoform_id
        ),]
        rownames(localCount) <- localCount$isoform_id

        ### add gene id
        localCount$gene_id <- switchAnalyzeRlist$isoformFeatures$gene_id[match(
            localCount$isoform_id, switchAnalyzeRlist$isoformFeatures$isoform_id
        )]

        ### extract isoform and gene annotation
        txInfo <- localCount[,c('isoform_id', 'gene_id')]

        ### extract counts
        localCount <- localCount[,setdiff(colnames(localCount), c('isoform_id', 'gene_id'))]

        ### this column always exist if the input was prepared using `importRdata`
        rownames(localDesign) <- localDesign$sampleID

        ### Make summarizedExperiment object
        sumExp <- SummarizedExperiment::SummarizedExperiment(
            assays = list(counts = localCount),
            colData = localDesign,
            rowData = txInfo
        )
    }

    ### Fitting linear models
    if(TRUE){
        if (!quiet) {
            message("Step 2 of 4: Fitting linear models...")
        }
        sumExp <- satuRn::fitDTU(
            object = sumExp,
            formula = localFormula,
            parallel = FALSE,
            BPPARAM = BiocParallel::bpparam(),
            verbose = showProgress
        )
    }

    ### Testing pairwise comparison(s)
    if(TRUE){
        if(!quiet){
            message("Step 3 of 4: Testing pairwise comparison(s)...")
        }

        # Set up contrast matrix L to perform all pairwise comparisons
        # For satuRn, this can be done easily, without looping on data splits
        L <- matrix(0, ncol = nrow(comparisonsToMake), nrow = ncol(localModel))
        rownames(L) <- colnames(localModel)
        colnames(L) <- paste0("Contrast_", seq_len(ncol(L)))
        for (i in 1:ncol(L)) {
            L[match(comparisonsToMake[i,], rownames(L)),i] <- c(-1,1)
        }

        ### Test all contrasts simultaneously
        sumExp <- satuRn::testDTU(
            object = sumExp,
            contrasts = L,
            diagplot1 = diagplots,
            diagplot2 = diagplots,
            sort = FALSE
        )
    }

    ### Preparing output
    if (TRUE) {
        if (!quiet) {
            message("\nStep 4 of 4: Preparing output...")
        }

        resultOfPairwiseTest <- rowData(sumExp)[grep("fitDTUResult",names(rowData(sumExp)))]
        for (i in 1:ncol(resultOfPairwiseTest)) {
            resultOfPairwiseTest[[i]]$condition_1 <- comparisonsToMake$condition_1[i]
            resultOfPairwiseTest[[i]]$condition_2 <- comparisonsToMake$condition_2[i]

            # ensure isoform_ids are propagated
            resultOfPairwiseTest[[i]]$isoform_id  <- rowData(sumExp)$isoform_id
        }

        ### check for which contrast(s) the empirical correction has failed
        empirical_succeed <- unlist(lapply(resultOfPairwiseTest, function(i){
            sum(!is.na(i$pval))
        }))
        if(any(empirical_succeed < 500)){ # if not all analyses allowed for computing empirical
            warning("Empirical adjustment of p-values failed, continuing analysis with raw p-values \n")
            for (i in 1:ncol(resultOfPairwiseTest)) { # work with raw p-values
                resultOfPairwiseTest[[i]]$padj <- resultOfPairwiseTest[[i]]$regular_FDR
            }
        } else { # if none failed, work with empirical p-values
            for (i in 1:ncol(resultOfPairwiseTest)) {
                resultOfPairwiseTest[[i]]$padj <- resultOfPairwiseTest[[i]]$empirical_FDR
            }
        }

        resultOfPairwiseTest <- S4Vectors::do.call(S4Vectors::rbind, resultOfPairwiseTest)
        rownames(resultOfPairwiseTest) <- NULL

        ### Replace with reference ids (which specify isoform-contrast combinations)
        resultOfPairwiseTest$iso_ref <-
            switchAnalyzeRlist$isoformFeatures$iso_ref[match(
                stringr::str_c(
                    resultOfPairwiseTest$isoform_id,
                    resultOfPairwiseTest$condition_1,
                    resultOfPairwiseTest$condition_2
                ),
                stringr::str_c(
                    switchAnalyzeRlist$isoformFeatures$isoform_id,
                    switchAnalyzeRlist$isoformFeatures$condition_1,
                    switchAnalyzeRlist$isoformFeatures$condition_2
                )
            )]

        ### Remove those without ID (below filtering treshold)
        resultOfPairwiseTest <- resultOfPairwiseTest[which(
            !is.na(resultOfPairwiseTest$iso_ref)
        ),]
        resultOfPairwiseTest$gene_ref <-
            switchAnalyzeRlist$isoformFeatures$gene_ref[match(
                resultOfPairwiseTest$iso_ref,
                switchAnalyzeRlist$isoformFeatures$iso_ref
            )]

        myDiff <- setdiff(
            colnames(resultOfPairwiseTest),
            c('iso_ref', 'gene_ref','isoform_id')
        )
        resultOfPairwiseTest <- resultOfPairwiseTest[,c('iso_ref', 'gene_ref','isoform_id', myDiff)]

        if(nrow(resultOfPairwiseTest) == 0) {
            stop('Something went wrong - please contact developers')
        }
    }

    ### Add results to switchAnalyzeRlist
    if (TRUE) {
        if (!quiet) {
            message("Result added switchAnalyzeRlist")
        }

        ### obtain gene-level results
        switchAnalyzeRlist$isoformFeatures$gene_switch_q_value <- NA
        switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <- NA
        geneQlevel <- sapply(X = split(resultOfPairwiseTest$padj,
                                       resultOfPairwiseTest$gene_ref), FUN = function(x) {
                                           min(c(1, x), na.rm = TRUE)
                                       })
        switchAnalyzeRlist$isoformFeatures$gene_switch_q_value <- geneQlevel[match(switchAnalyzeRlist$isoformFeatures$gene_ref,
                                                                                   names(geneQlevel))]
        switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <- resultOfPairwiseTest$padj[match(switchAnalyzeRlist$isoformFeatures$iso_ref,
                                                                                                     resultOfPairwiseTest$iso_ref)]

        ### add results to switchAnalyzeRlist
        switchAnalyzeRlist$isoformSwitchAnalysis <- resultOfPairwiseTest

        ### reduce to genes with at least one significant isoform if TRUE
        if (reduceToSwitchingGenes) {
            if (reduceFurtherToGenesWithConsequencePotential) {
                tmp <- extractSwitchPairs(
                    switchAnalyzeRlist,
                    alpha = alpha,
                    dIFcutoff = dIFcutoff,
                    onlySigIsoforms = onlySigIsoforms
                )
                combinedGeneIDsToKeep <- unique(tmp$gene_ref)
            } else {
                isoResTest <- any(!is.na(switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value))
                if (isoResTest) {
                    combinedGeneIDsToKeep <- switchAnalyzeRlist$isoformFeatures$gene_ref[which(switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <
                                                                                                   alpha & abs(switchAnalyzeRlist$isoformFeatures$dIF) >
                                                                                                   dIFcutoff)]
                } else {
                    combinedGeneIDsToKeep <- switchAnalyzeRlist$isoformFeatures$gene_ref[which(switchAnalyzeRlist$isoformFeatures$gene_switch_q_value <
                                                                                                   alpha & abs(switchAnalyzeRlist$isoformFeatures$dIF) >
                                                                                                   dIFcutoff)]
                }
            }
            if (length(combinedGeneIDsToKeep) == 0) {
                stop(paste("No signifcant switches were found with the supplied cutoffs",
                           "whereby we cannot reduce the switchAnalyzeRlist to only",
                           "significant genes (with consequence potential)"))
            }
            if (keepIsoformInAllConditions) {
                combinedGeneIDsToKeep <- switchAnalyzeRlist$isoformFeatures$gene_id[which(switchAnalyzeRlist$isoformFeatures$gene_ref %in%
                                                                                              combinedGeneIDsToKeep)]
                switchAnalyzeRlist <- subsetSwitchAnalyzeRlist(switchAnalyzeRlist,
                                                               switchAnalyzeRlist$isoformFeatures$gene_id %in%
                                                                   combinedGeneIDsToKeep)
            }
            if (!keepIsoformInAllConditions) {
                switchAnalyzeRlist <- subsetSwitchAnalyzeRlist(switchAnalyzeRlist,
                                                               switchAnalyzeRlist$isoformFeatures$gene_ref %in%
                                                                   combinedGeneIDsToKeep)
            }
        }
    }

    if (!quiet) {
        message("Done")
    }
    return(switchAnalyzeRlist)
}

### Test via DEXSeq
isoformSwitchTestDEXSeq <- function(
    ### Core arguments
    switchAnalyzeRlist,
    alpha = 0.05,
    dIFcutoff = 0.1,

    ### Advanced arguments
    reduceToSwitchingGenes = TRUE,
    reduceFurtherToGenesWithConsequencePotential = TRUE,
    onlySigIsoforms = FALSE,
    keepIsoformInAllConditions = TRUE,
    showProgress = TRUE,
    quiet = FALSE
) {
    tStart <- Sys.time()

    ### Test input
    if (TRUE) {
        ### Tjek arguments
        if(TRUE) {
            if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist')        {
                stop(paste(
                    'The object supplied to \'switchAnalyzeRlist\'',
                    'must be a \'switchAnalyzeRlist\''
                ))
            }
            if( switchAnalyzeRlist$sourceId == 'preDefinedSwitches' ) {
                stop(
                    paste(
                        'The switchAnalyzeRlist is made from pre-defined isoform switches',
                        'which means it is made without defining conditions (as it should be).',
                        '\nThis also means it cannot used to test for new isoform switches -',
                        'if that is your intention you need to create a new',
                        'switchAnalyzeRlist with the importRdata() function and try again.',
                        sep = ' '
                    )
                )
            }

            if (!is.logical(reduceToSwitchingGenes))  {
                stop(paste(
                    'The argument supplied to \'reduceToSwitchingGenes\'',
                    'must be an a logic'
                ))
            }
            if (dIFcutoff < 0 | dIFcutoff > 1) {
                stop('The dIFcutoff must be in the interval [0,1].')
            }
            if (reduceToSwitchingGenes) {
                if (alpha < 0 |
                    alpha > 1) {
                    stop('The alpha parameter must be between 0 and 1 ([0,1]).')
                }
                if (alpha > 0.05) {
                    warning(paste(
                        'Most journals and scientists consider an alpha',
                        'larger than 0.05 untrustworthy. We therefore recommend',
                        'using alpha values smaller than or queal to 0.05'
                    ))
                }
            }

            if( any(switchAnalyzeRlist$conditions$nrReplicates == 1)) {
                stop(
                    paste(
                        'A statistical test cannot be performed without replicates.',
                        'Please remove all conditions with only 1 replicate and try again.',
                        '\nThis can either be done before creating the switchAnalyzeRlist',
                        'in the first place or with the subsetSwitchAnalyzeRlist() function',
                        sep=' '
                    )
                )
            }

            if(any( switchAnalyzeRlist$conditions$nrReplicates > 5)) {
                warning(paste(
                    'You seem to have many replicates. We therfore recomend you use isoformSwitchTestSatuRn() instead.\n',
                    'It is just as accurate and much faster - especially for larger datasets.'
                ))
            }
        }

        ### Setup progress bar
        comaprisonsToMake <- unique(switchAnalyzeRlist$isoformFeatures[, c(
            'condition_1', 'condition_2'
        )])
        if (showProgress &  !quiet &  nrow(comaprisonsToMake) > 1) {
            progressBar <- 'text'
        } else {
            progressBar <- 'none'
        }

        ### Check Input
        countsAvailable <- !is.null( switchAnalyzeRlist$isoformCountMatrix)

        if( ! countsAvailable ) {
            stop('isoformSwitchTestDEXSeq requires count data to work. Please remake switchAnalyzeRlist and try again')
        }

        ### For messages
        nrAnalysis <- 2
        analysisDone <- 1
    }

    ### If nessesary (re-)calculate IF matrix
    if(TRUE) {
        localIF <- switchAnalyzeRlist$isoformRepIF
        rownames(localIF) <- localIF$isoform_id
        localIF$isoform_id <- NULL
    }

    ### Extract mean IFs
    if(TRUE) {
        ifMeans <- rowMeans(
            localIF[,switchAnalyzeRlist$designMatrix$sampleID,drop=FALSE],
            na.rm = TRUE
        )

        ifMeanList <- plyr::dlply(
            .data = switchAnalyzeRlist$designMatrix,
            .variables = 'condition',
            .fun = function(aDF) { # aDF <- switchAnalyzeRlist$designMatrix[1:2,]
                rowMeans(localIF[,aDF$sampleID,drop=FALSE], na.rm = TRUE)
            }
        )
    }

    ### For each pariwise comparison build and test with DEXSeq
    if(TRUE) {
        ### Estimate time
        if (!quiet) {
            message(paste(
                'Step ',
                analysisDone,
                ' of ',
                nrAnalysis,
                ': Testing each pairwise comparisons with DEXSeq (this might be a bit slow)...',
                sep=''
            ))

            ### Estimate runtime time
            if(TRUE) {
                expectedTime <- plyr::ddply(
                    .data = comaprisonsToMake,
                    .variables = c('condition_1','condition_2'),
                    .progress = 'none',
                    .fun = function(aComp) { # aComp <- comaprisonsToMake[1,]
                        sampleOverview <- switchAnalyzeRlist$conditions[which(
                            switchAnalyzeRlist$conditions$condition %in% unlist(aComp)
                        ),]

                        nIso <- length(switchAnalyzeRlist$isoformFeatures$isoform_id[which(
                            switchAnalyzeRlist$isoformFeatures$condition_1 == aComp$condition_1 &
                            switchAnalyzeRlist$isoformFeatures$condition_2 == aComp$condition_2
                        )])

                        localExp <- data.frame(
                            samples=sum(sampleOverview$nrReplicates),
                            n_iso = nIso
                        )

                        # Esitmate based on benchmarking from satuRn paper
                        localExp$min <- max(c(
                            1,  # do not repport less than 1 min
                            10^(
                                -3.5574003487257 +
                                (1.0190040124555 * log10(localExp$n_iso)) +
                                (0.0212879246701 * localExp$samples)
                            )
                        ))

                        return(localExp)
                    }
                )

                expectedTime <- sum(expectedTime$min)

                message(paste(
                    '    Estimated run time is:',
                    round(expectedTime, digits = 1),
                    'min',
                    sep=' '
                ))


            }
        }

        ### Do pairwise test
        dexseqPairwiseResults <- plyr::ddply(
            .data = comaprisonsToMake,
            .variables = c('condition_1','condition_2'),
            .progress = progressBar,
            .fun = function(aComp) { # aComp <- comaprisonsToMake[1,]
                ### Extract data
                if(TRUE) {
                    ### Extract local design
                    designSubset <- switchAnalyzeRlist$designMatrix[which(
                        switchAnalyzeRlist$designMatrix$condition %in% unlist(aComp)
                    ),]
                    designSubset$condition <- factor(
                        designSubset$condition,
                        levels = unique(c(
                            comaprisonsToMake$condition_1,
                            comaprisonsToMake$condition_2
                        )
                    ))

                    ### Extract corresponding annotation
                    localAnnot <- switchAnalyzeRlist$isoformFeatures[
                        which(
                            switchAnalyzeRlist$isoformFeatures$condition_1 == aComp$condition_1 &
                            switchAnalyzeRlist$isoformFeatures$condition_2 == aComp$condition_2
                        ),
                        c('isoform_id','iso_ref','gene_ref')
                    ]

                    ### Extract corresponding count data
                    localCount <- switchAnalyzeRlist$isoformCountMatrix[
                        which(
                            switchAnalyzeRlist$isoformCountMatrix$isoform_id %in%
                                localAnnot$isoform_id
                        ),
                        which(
                            colnames(switchAnalyzeRlist$isoformCountMatrix) %in%
                                c('isoform_id', designSubset$sampleID)
                        )
                    ]


                    ### Massage for DEXSeq
                    localCount$gene_ref <- localAnnot$gene_ref[match(
                        localCount$isoform_id, localAnnot$isoform_id
                    )]

                    localCount$iso_ref <- localAnnot$iso_ref[match(
                        localCount$isoform_id, localAnnot$isoform_id
                    )]


                    localCount$isoform_id <- NULL
                    rownames(localCount) <- localCount$iso_ref
                }

                ### Make formulas
                if(TRUE) {
                    ### Convert group of interest to factors
                    designSubset$condition <- factor(
                        designSubset$condition,
                        levels = unique(c(
                            comaprisonsToMake$condition_1,
                            comaprisonsToMake$condition_2
                        ))
                    )
                    colnames(designSubset)[1] <- 'sample'

                    ### remove non-varying features
                    if( ncol(designSubset) > 2 ) {
                        coreDesign <- designSubset[,1:2]

                        batchDesign <- designSubset[,3:ncol(designSubset), drop=FALSE]

                        batchDesignReduced <- batchDesign[,which(
                            apply(batchDesign, 2, function(x) length(unique(x))) > 1
                        ),drop=FALSE]

                        designSubset <- cbind(
                            coreDesign,
                            batchDesignReduced
                        )
                    }

                    ### remove completly co-linear features - to dangerous - could be true effects
                    if(FALSE) {
                        if( ncol(designSubset) > 2 ) {
                            toKeep <- integer()
                            for(i in 3:ncol(designSubset) ) { # i <- 3
                                if(
                                    ! identical(
                                        as.integer(as.factor(designSubset$condition)),
                                        as.integer(as.factor(designSubset[,i]))
                                    )
                                ) {
                                    toKeep <- c(toKeep, i)
                                }
                            }

                            designSubset <- designSubset[,c(1:2,toKeep)]
                        }
                    }

                    ### Check co-founders for group vs continous variables
                    if( ncol(designSubset) > 2 ) {

                        for(i in 3:ncol(designSubset) ) { # i <- 4
                            if( class(designSubset[,i]) %in% c('numeric', 'integer') ) {
                                if( uniqueLength( designSubset[,i] ) * 2 <= length(designSubset[,i]) ) {
                                    designSubset[,i] <- factor(designSubset[,i])
                                }
                            } else {
                                designSubset[,i] <- factor(designSubset[,i])
                            }
                        }
                    }

                    ### Make formula for model (exon reads as "isoform")
                    basicFormula <- '~ sample + exon'
                    if (ncol(designSubset) > 2) {
                        basicFormula <- paste(
                            basicFormula,
                            '+',
                            paste(
                                paste0(
                                    'exon:', colnames(designSubset)[3:ncol(designSubset)]
                                ),
                                collapse = ' + '
                            ),
                            sep=' '
                        )
                    }

                    ### Make full formula
                    fullFormula <- paste(basicFormula, '+ condition:exon', sep=' ')

                    fullFormula <- as.formula( fullFormula )
                    basicFormula <- as.formula( basicFormula )
                }

                ### Test with DEXSeq
                if(TRUE) {
                    ### Create data object
                    suppressWarnings(
                        suppressMessages(
                            dexList <- DEXSeqDataSet(
                                countData  = round(localCount[,which( ! colnames(localCount) %in% c('gene_ref','iso_ref'))]),      # cannot handle non-integers
                                alternativeCountData = NULL,
                                sampleData = designSubset,
                                design     = fullFormula,
                                featureID  = localCount$iso_ref,
                                groupID    = localCount$gene_ref
                            )
                        )
                    )

                    ### Estimate parmeters
                    suppressMessages(
                        dexList <- DEXSeq::estimateSizeFactors(dexList)
                    )
                    suppressMessages(
                        dexList <- DEXSeq::estimateDispersions(dexList, quiet=TRUE)
                    )

                    ### Test
                    suppressMessages(
                        dexList <- testForDEU(dexList, reducedModel = basicFormula)
                    )
                }

                ### Extract and augment test result
                if(TRUE) {
                    dexRes <- as.data.frame(
                        DEXSeqResults(dexList, independentFiltering=FALSE)
                    )
                    dexRes <- dexRes[,c('groupID','featureID','pvalue','padj')] # fdr corrected
                    #dexRes <- dexRes[which( !is.na(dexRes$pvalue)),] # outcommented 9th october 2018 by KVS
                    rownames(dexRes) <- NULL

                    colnames(dexRes)[1:2] <- c('gene_ref','iso_ref')

                    dexRes$isoform_id <- switchAnalyzeRlist$isoformFeatures$isoform_id[match(
                        dexRes$iso_ref,
                        switchAnalyzeRlist$isoformFeatures$iso_ref
                    )]

                    ### Add means IFs
                    if(TRUE) {
                        dexRes$IF1 <- ifMeanList[[ aComp$condition_1 ]][
                            match(
                                dexRes$isoform_id,
                                names(ifMeanList[[ aComp$condition_1 ]])
                            )
                            ]
                        dexRes$IF2 <- ifMeanList[[ aComp$condition_2 ]][
                            match(
                                dexRes$isoform_id,
                                names(ifMeanList[[ aComp$condition_2 ]])
                            )
                            ]
                        dexRes$dIF <- dexRes$IF2 - dexRes$IF1
                    }

                }

                return(dexRes)
            }
        )

        ### Reorder
        ofInterest <- c('iso_ref','gene_ref','isoform_id','condition_1','condition_2','dIF')
        desiredOrder <- c(
            ofInterest,
            setdiff(colnames(dexseqPairwiseResults), ofInterest)
        )

        dexseqPairwiseResults <- dexseqPairwiseResults[,match(
            desiredOrder, colnames(dexseqPairwiseResults)
        )]

        ### Update counter
        analysisDone <- analysisDone + 1

    }

    ### Integrate with switchList
    if(TRUE) {
        if (!quiet) {
            message(paste(
                'Step ',
                analysisDone,
                ' of ',
                nrAnalysis,
                ': Integrating result into switchAnalyzeRlist...',
                sep=''
            ))
        }

        ### Test result
        if(TRUE) {
            ### Overwrite previous results
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value <- NA
            switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <- NA

            ### summarize to gene level
            geneQlevel <- sapply(
                X = split(
                    dexseqPairwiseResults$padj,
                    dexseqPairwiseResults$gene_ref
                ),
                FUN = function(x) {
                    min(c(1, x), na.rm = TRUE)
                }
            )
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value <-
                geneQlevel[match(
                    switchAnalyzeRlist$isoformFeatures$gene_ref,
                    names(geneQlevel)
                )]

            ### Isoform level
            switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <-
                dexseqPairwiseResults$padj[match(
                    switchAnalyzeRlist$isoformFeatures$iso_ref,
                    dexseqPairwiseResults$iso_ref
                )]
        }

        switchAnalyzeRlist$isoformSwitchAnalysis <- dexseqPairwiseResults
    }

    ### Print status
    if (!quiet) {
        myN <-
            length(unique(switchAnalyzeRlist$isoformSwitchAnalysis$gene_ref))
        myFrac <-
            myN / length(unique(
                switchAnalyzeRlist$isoformFeatures$gene_ref
            )) * 100
        message(
            paste(
                '    Isoform switch analysis was performed for ',
                myN,
                ' gene comparisons (',
                round(myFrac, digits = 1),
                '%).',
                sep = ''
            )
        )
    }

    ### Reduce to genes with switches
    if (reduceToSwitchingGenes ) {
        if( reduceFurtherToGenesWithConsequencePotential ) {
            tmp <- extractSwitchPairs(
                switchAnalyzeRlist,
                alpha = alpha,
                dIFcutoff = dIFcutoff,
                onlySigIsoforms = onlySigIsoforms
            )
            combinedGeneIDsToKeep <- unique(tmp$gene_ref)

        } else {
            combinedGeneIDsToKeep <- unique(
                switchAnalyzeRlist$isoformFeatures$gene_ref[which(
                    switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <
                        alpha &
                        abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
                )]
            )
        }

        if (length(combinedGeneIDsToKeep) == 0) {
            stop(paste(
                'No signifcant switches were found with the supplied cutoffs',
                'whereby we cannot reduce the switchAnalyzeRlist to only',
                'significant genes (with consequence potential)'
            ))
        }

        if(   keepIsoformInAllConditions ) {
            combinedGeneIDsToKeep <- switchAnalyzeRlist$isoformFeatures$gene_id[which(
                switchAnalyzeRlist$isoformFeatures$gene_ref %in% combinedGeneIDsToKeep
            )]

            switchAnalyzeRlist <-
                subsetSwitchAnalyzeRlist(
                    switchAnalyzeRlist,
                    switchAnalyzeRlist$isoformFeatures$gene_id %in% combinedGeneIDsToKeep
                )
        }
        if( ! keepIsoformInAllConditions ) {
            switchAnalyzeRlist <-
                subsetSwitchAnalyzeRlist(
                    switchAnalyzeRlist,
                    switchAnalyzeRlist$isoformFeatures$gene_ref %in% combinedGeneIDsToKeep
                )
        }
    }

    if (!quiet) {
        message(paste('Total runtime:' ,round( difftime(Sys.time(), tStart, units = 'min'), digits = 2), 'min', sep=' '))
        message('Done')
    }
    return(switchAnalyzeRlist)
}


### Summarize switching
extractSwitchSummary <- function(
    switchAnalyzeRlist,
    filterForConsequences = FALSE,
    alpha = 0.05,
    dIFcutoff = 0.1,
    onlySigIsoforms = FALSE,
    includeCombined = nrow(unique(
        switchAnalyzeRlist$isoformFeatures[, c('condition_1', 'condition_1')]
    )) > 1
) {
    ### Test input
    if (TRUE) {
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(paste(
                'The object supplied to \'switchAnalyzeRlist\' must',
                'be a \'switchAnalyzeRlist\''
            ))
        }
        if (!any(!is.na(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value
        ))) {
            stop(paste(
                'The analsis of isoform switching must be performed before',
                'functional consequences can be analyzed. Please run \'isoformSwitchTestDEXSeq\' or \'isoformSwitchTestSatuRn\' and try again.'
            ))
        }
        if (filterForConsequences) {
            if (!'switchConsequencesGene' %in%
                colnames(switchAnalyzeRlist$isoformFeatures)) {
                stop(paste(
                    'The switchAnalyzeRlist does not contain isoform',
                    'switching analysis. Please run the',
                    '\'isoformSwitchTestDEXSeq\' or \'isoformSwitchTestSatuRn\' function first.'
                ))
            }
        }
        if (alpha < 0 |
            alpha > 1) {
            warning('The alpha parameter must be between 0 and 1 ([0,1]).')
        }
        if (alpha > 0.05) {
            warning(paste(
                'Most journals and scientists consider an alpha larger',
                'than 0.05 untrustworthy. We therefore recommend using',
                'alpha values smaller than or queal to 0.05'
            ))
        }
        if (dIFcutoff < 0 | dIFcutoff > 1) {
            stop('The dIFcutoff must be in the interval [0,1].')
        }
    }

    nrDf <-
        unique(switchAnalyzeRlist$isoformFeatures[, c(
            'condition_1', 'condition_2'
        )])

    nrDf <-
        data.frame(
            Comparison = paste(
                nrDf$condition_1,
                nrDf$condition_2, sep = ' vs '),
            nrIsoforms = 0,
            nrSwitches = 0,
            nrGenes = 0,
            stringsAsFactors = FALSE
        )

    if(includeCombined) {
        nrDf <- rbind(
            nrDf,
            data.frame(
                Comparison = 'Combined',
                nrIsoforms = 0,
                nrSwitches = 0,
                nrGenes = 0,
                stringsAsFactors = FALSE
            )
        )
    }



    ### Count genes and isoforms
    if(TRUE) {
        ### Extract data needed
        columnsToExtract <-
            c(
                'isoform_id',
                'gene_id',
                'condition_1',
                'condition_2',
                'dIF',
                'isoform_switch_q_value',
                'gene_switch_q_value',
                'switchConsequencesGene'
            )

        isoResTest <-
            any(!is.na(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
            ))
        if (isoResTest) {
            dataDF <- switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value < alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))]
        } else {
            dataDF <- switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$gene_switch_q_value    < alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))]
        }
        if (nrow(dataDF)) {
            if (filterForConsequences) {
                dataDF <- dataDF[which(dataDF$switchConsequencesGene), ]
            }
            if (nrow(dataDF) ) {

                ### Add levels
                dataDF$comparison <- paste(dataDF$condition_1, dataDF$condition_2, sep = ' vs ')
                dataDF$comparison <- factor(
                    x = dataDF$comparison,
                    levels = unique(setdiff(nrDf$Comparison, 'Combined'))
                )

                ### Summarize pr comparison
                dataList <- split(
                    dataDF,
                    f = dataDF$comparison,
                    drop = FALSE
                )

                if (includeCombined) {
                    dataList$Combined <- dataDF
                }

                isoResTest <-
                    any(!is.na(
                        switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
                    ))
                myNumbers <- plyr::ldply(dataList, function(aDF) {
                    if (isoResTest) {
                        sigGenes <-
                            unique(aDF$gene_id   [which(aDF$isoform_switch_q_value < alpha &
                                                            abs(aDF$dIF) > dIFcutoff)])
                        sigIso   <-
                            unique(aDF$isoform_id[which(aDF$isoform_switch_q_value < alpha &
                                                            abs(aDF$dIF) > dIFcutoff)])
                        return(data.frame(
                            nrIsoforms = length(sigIso),
                            nrGenes = length(sigGenes)
                        ))
                    } else {
                        sigGenes <-
                            unique(aDF$gene_id   [which(aDF$gene_switch_q_value < alpha &
                                                            abs(aDF$dIF) > dIFcutoff)])
                        sigIso   <-
                            unique(aDF$isoform_id[which(aDF$gene_switch_q_value < alpha &
                                                            abs(aDF$dIF) > dIFcutoff)])
                        return(data.frame(
                            nrIsoforms = length(sigIso),
                            nrGenes = length(sigGenes)
                        ))
                    }


                })
                colnames(myNumbers)[1] <- 'Comparison'

                ### Overwrite numbers
                nrDf$nrIsoforms <- myNumbers$nrIsoforms[match(
                    nrDf$Comparison, myNumbers$Comparison
                )]

                nrDf$nrGenes <- myNumbers$nrGenes[match(
                    nrDf$Comparison, myNumbers$Comparison
                )]

            }

        }
    }

    ### Count switches
    if(TRUE) {
        localSwitches <- extractSwitchPairs(
            switchAnalyzeRlist = switchAnalyzeRlist,
            alpha = alpha,
            dIFcutoff = dIFcutoff,
            onlySigIsoforms = onlySigIsoforms
        )
        localSwitches$Comparison = stringr::str_c(
            localSwitches$condition_1,
            ' vs ',
            localSwitches$condition_2
        )

        ### Set levels
        localSwitches$Comparison <- factor(
            localSwitches$Comparison,
            levels = unique(nrDf$Comparison)
        )

        ### Filter for consequences
        if(filterForConsequences) {
            refConseqGenes <- switchAnalyzeRlist$isoformFeatures$gene_ref[which(
                switchAnalyzeRlist$isoformFeatures$switchConsequencesGene
            )]

            localSwitches <- localSwitches[which(
                localSwitches$gene_ref %in% refConseqGenes
            ),]
        }

        ### Count
        nrSwitches <- plyr::ddply(
            .data = localSwitches,
            .variables = 'Comparison',
            .drop = FALSE,
            function(x) {
                data.frame(nrSwitches = nrow(x))
            }
        )

        if(includeCombined) {
            nrSwitches$nrSwitches[which(
                nrSwitches$Comparison == 'Combined'
            )] <- nrow( unique(
                localSwitches[,c('isoformUpregulated','isoformDownregulated')]
            ))
        }


        nrDf$nrSwitches <- nrSwitches$nrSwitches[match(
            nrDf$Comparison, nrSwitches$Comparison
        )]
    }

    return(nrDf)
}

extractSwitchOverlap <- function(
    switchAnalyzeRlist,
    filterForConsequences = FALSE,
    alpha = 0.05,
    dIFcutoff = 0.1,
    scaleVennIfPossible=TRUE,
    plotIsoforms = TRUE,
    plotSwitches = TRUE,
    plotGenes = TRUE
) {
    ### Test input
    if (TRUE) {
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(paste(
                'The object supplied to \'switchAnalyzeRlist\' must',
                'be a \'switchAnalyzeRlist\''
            ))
        }
        if (!any(!is.na(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value
        ))) {
            stop(paste(
                'The analsis of isoform switching must be performed before',
                'functional consequences can be analyzed. Please run \'isoformSwitchTestDEXSeq\' or \'isoformSwitchTestSatuRn\' and try again.'
            ))
        }
        if (filterForConsequences) {
            if (!'switchConsequencesGene' %in%
                colnames(switchAnalyzeRlist$isoformFeatures)) {
                stop(paste(
                    'The switchAnalyzeRlist does not contain isoform',
                    'switching analysis. Please run the',
                    '\'isoformSwitchTestDEXSeq\' or \'isoformSwitchTestSatuRn\' function first.'
                ))
            }
        }
        if (alpha < 0 |
            alpha > 1) {
            warning('The alpha parameter must be between 0 and 1 ([0,1]).')
        }
        if (alpha > 0.05) {
            warning(paste(
                'Most journals and scientists consider an alpha larger',
                'than 0.05 untrustworthy. We therefore recommend using',
                'alpha values smaller than or queal to 0.05'
            ))
        }
        if (dIFcutoff < 0 | dIFcutoff > 1) {
            stop('The dIFcutoff must be in the interval [0,1].')
        }

        nCon <- nrow(unique( switchAnalyzeRlist$isoformFeatures[,c('condition_1','condition_2')]))
        if( nCon > 5 ) {
            stop('Venn Diagrams unfortunatly only support up to 5 comparisons')
        }
        if( nCon < 2 ) {
            stop('One cannot make a Venn Diagram with only one comparison')
        }

        if( sum(c(plotIsoforms, plotSwitches, plotGenes)) < 1) {
            stop('You need to to plot at least one plot')
        }

    }

    ### Helper function
    mfGGplotColors <- function(n) {
        hues = seq(15, 375, length=n+1)
        hcl(h=hues, l=65, c=100)[1:n]
    }

    ### Extract and massage data
    if(TRUE) {
        columnsToExtract <-
            c(
                'isoform_id',
                'gene_ref',
                'gene_id',
                'condition_1',
                'condition_2',
                'dIF',
                'isoform_switch_q_value',
                'gene_switch_q_value',
                'switchConsequencesGene'
            )

        isoResTest <-
            any(!is.na(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
            ))
        if (isoResTest) {
            dataDF <- switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value < alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))]
        } else {
            dataDF <- switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$gene_switch_q_value    < alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))]
        }
        if (nrow(dataDF) == 0) {
            return(backUpDf)
        }

        isoPairs <- extractSwitchPairs(
            switchAnalyzeRlist,
            alpha = alpha,
            dIFcutoff = dIFcutoff,
            onlySigIsoforms = FALSE
        )

        if (filterForConsequences) {
            dataDF <- dataDF[which(dataDF$switchConsequencesGene), ]
            if (nrow(dataDF) == 0) {

                backUpDf <-
                    unique(switchAnalyzeRlist$isoformFeatures[, c(
                        'condition_1', 'condition_2'
                    )])
                backUpDf <-
                    data.frame(
                        Comparison = paste(
                            backUpDf$condition_1,
                            backUpDf$condition_2, sep = ' vs '),
                        nrIsoforms = 0,
                        nrSwitches = 0,
                        nrGenes = 0,
                        stringsAsFactors = FALSE
                    )
                return(backUpDf)
            }

            isoPairs <- isoPairs[which(
                isoPairs$gene_ref %in% dataDF$gene_ref
            ),]
        }

        isoPairs$switch <- stringr::str_c(
            isoPairs$isoformDownregulated,
            isoPairs$isoformUpregulated
        )

        isoPairs$comparison <- stringr::str_c(
            isoPairs$condition_1,
            '\nvs\n',
            isoPairs$condition_2
        )

        dataDF$comparison <- stringr::str_c(
            dataDF$condition_1,
            '\nvs\n',
            dataDF$condition_2
        )

        geneList <- split(dataDF$gene_id   , dataDF$comparison)
        isoList  <- split(
            stringr::str_c(dataDF$isoform_id, sign(dataDF$dIF)),
            dataDF$comparison
        )
        switchList <- split(isoPairs$switch   , isoPairs$comparison)

    }

    ### Assign colors and alpha
    n <- length(isoList)
    vennColors <- mfGGplotColors(n)
    localAlpha <- ifelse(n==2, 0.6, 0.4)

    ### Suppress venn log files
    futile.logger::flog.threshold(futile.logger::ERROR, name = "VennDiagramLogger")

    ### Make venn diagrams
    if(TRUE) {
        switchVenn <- VennDiagram::venn.diagram(
            x = switchList,
            euler.d=scaleVennIfPossible,
            scale=scaleVennIfPossible,
            col='transparent',
            alpha=localAlpha,
            fill=vennColors,
            filename=NULL,
            main='Overlap in Switches'
        )
        geneVenn <- VennDiagram::venn.diagram(
            x = geneList,
            euler.d=scaleVennIfPossible,
            scale=scaleVennIfPossible,
            col='transparent',
            alpha=localAlpha,
            fill=vennColors,
            filename=NULL,
            main='Overlap in Switching Genes'
        )
        isoVenn <- VennDiagram::venn.diagram(
            x = isoList,
            euler.d=scaleVennIfPossible,
            scale=scaleVennIfPossible,
            col='transparent',
            alpha=localAlpha,
            fill=vennColors,
            filename=NULL,
            main='Overlap in Isoforms'
        )

    }

    ### Figure out plot size
    nPlots <- sum(c(plotIsoforms, plotSwitches, plotGenes))
    nCols <- 3*nPlots + nPlots - 1
    positionToStart <- 1

    ### Set up viewport
    grid::grid.newpage()
    grid::pushViewport(
        grid::plotViewport(
            layout=grid::grid.layout(1, nCols)
        )
    )

    ### Plot venn diagrams
    if(plotIsoforms) {
        grid::pushViewport(
            grid::plotViewport(
                layout.pos.col = positionToStart:(positionToStart+2)
            )
        )
        grid::grid.draw(isoVenn)
        grid::popViewport()

        positionToStart <- positionToStart + 4
    }

    if(plotSwitches) {
        grid::pushViewport(
            grid::plotViewport(
                layout.pos.col = positionToStart:(positionToStart+2)
            )
        )
        grid::grid.draw(switchVenn)
        grid::popViewport()

        positionToStart <- positionToStart + 4
    }

    if(plotGenes) {
        grid::pushViewport(
            grid::plotViewport(
                layout.pos.col = positionToStart:(positionToStart+2)
            )
        )
        grid::grid.draw(geneVenn)
        grid::popViewport()
    }

}

### Extract
extractTopSwitches <- function(
    switchAnalyzeRlist,
    filterForConsequences = FALSE,
    extractGenes = TRUE,
    alpha = 0.05,
    dIFcutoff = 0.1,
    n = 10,
    inEachComparison = FALSE,
    sortByQvals = TRUE
) {
    ### Test input
    if (TRUE) {
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(paste(
                'The object supplied to \'switchAnalyzeRlist\'',
                'must be a \'switchAnalyzeRlist\''
            ))
        }
        if (!any(!is.na(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value
        ))) {
            stop(paste(
                'The analsis of isoform switching must be performed before',
                'functional consequences can be analyzed.',
                'Please run \'isoformSwitchTestDEXSeq\' or \'isoformSwitchTestSatuRn\' and try again.'
            ))
        }
        if (filterForConsequences) {
            if (!'switchConsequencesGene' %in%
                colnames(switchAnalyzeRlist$isoformFeatures)) {
                stop(paste(
                    'The switchAnalyzeRlist does not contain the consequence prediction',
                    'whereby the "filterForConsequences" argument cannot be used.',
                    '\nIf that is what you want you need to run the \'analyzeSwitchConsequences()\' function first',
                    sep = ' '
                ))
            }
        }
        if (alpha < 0 |
            alpha > 1) {
            warning('The alpha parameter should usually be between 0 and 1 ([0,1]).')
        }
        if (alpha > 0.05) {
            warning(paste(
                'Most journals and scientists consider an alpha larger',
                'than 0.05 untrustworthy. We therefore recommend using',
                'alpha values smaller than or queal to 0.05'
            ))
        }
        if (dIFcutoff < 0 | dIFcutoff > 1) {
            stop('The dIFcutoff must be in the interval [0,1].')
        }
        if (!is.logical(extractGenes)) {
            stop('The extractGenes argument supplied must be a logical')
        }
        if (!is.logical(inEachComparison)) {
            stop('The inEachComparison argument supplied must be a logical')
        }
        if (!is.logical(sortByQvals)) {
            stop('The sortByQvals argument supplied must be a logical')
        }

        if( is.na(n) ) {
            n <- Inf
        }
    }

    if (extractGenes) {
        columnsToExtract <-
            c(
                'gene_ref',
                'gene_id',
                'gene_name',
                'condition_1',
                'condition_2',
                'gene_switch_q_value',
                'switchConsequencesGene',
                'dIF'
            )
        isoResTest <-
            any(!is.na(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
            ))
        if (isoResTest) {
            dataDF <- unique(switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <
                    alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))])
        } else {
            dataDF <- unique(switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$gene_switch_q_value <
                    alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))])
        }
        if (nrow(dataDF) == 0) {
            stop('No significant switching genes were found.')
        }

        if (filterForConsequences) {
            dataDF <- dataDF[which(dataDF$switchConsequencesGene), ]
        }
        if (nrow(dataDF) == 0) {
            stop('No significant switching genes consequences were found.')
        }

        ### Sort data
        dataDF$comparison <-
            paste(dataDF$condition_1, '_vs_', dataDF$condition_2, sep = '')
        if (is.na(sortByQvals)) {
            dataDF2 <- dataDF
            rownames(dataDF2) <- NULL
        } else if (sortByQvals) {
            dataDF2 <-
                dataDF[sort.list(
                    dataDF$gene_switch_q_value,
                    decreasing = FALSE
                ), ]
            rownames(dataDF2) <- NULL
        } else {
            combinedDif <-
                split(abs(dataDF$dIF), f = dataDF$gene_ref)
            combinedDif <- sapply(combinedDif, sum)

            ### Add to df
            dataDF$combinedDIF <-
                combinedDif[match(dataDF$gene_ref , names(combinedDif))]

            dataDF2 <-
                dataDF[sort.list(dataDF$combinedDIF, decreasing = TRUE), ]
            rownames(dataDF2) <- NULL
        }

        ### reduce to collumns wanted
        columnsToExtract <-
            c(
                'gene_ref',
                'gene_id',
                'gene_name',
                'condition_1',
                'condition_2',
                'gene_switch_q_value',
                'combinedDIF',
                'switchConsequencesGene'
            )
        dataDF2 <-
            unique(dataDF2[, na.omit(
                match(columnsToExtract, colnames(dataDF))
            )])


        ### Reduce to the number wanted
        if (! is.infinite(n)) {
            if (inEachComparison) {
                dataDF2$comparison <-
                    paste(dataDF2$condition_1,
                          '_vs_',
                          dataDF2$condition_2,
                          sep = '')
            } else {
                dataDF2$comparison <- 'AllCombined'
            }

            dataDF2 <-
                plyr::ddply(
                    .data = dataDF2,
                    .variables = 'comparison',
                    .fun = function(aDF) {
                        if ( n > nrow(dataDF2) ) {
                            if (filterForConsequences) {
                                warning(paste(
                                    'Less than',n ,'genes with significant',
                                    'switches and consequences were found.',
                                    'Returning those.'
                                ))
                            } else {
                                warning(paste(
                                    'Less than',n ,'genes genes with significant',
                                    'switches were found. Returning those.'
                                ))
                            }
                            n2 <- nrow(dataDF2)
                        } else {
                            n2 <- n
                        }

                        return(aDF[1:n2, ])
                    }
                )

            dataDF2 <- dataDF2[which(
                !is.na(dataDF2$gene_ref)
            ),]

            dataDF2$comparison <- NULL

        }

        dataDF2$Rank <- 1:nrow(dataDF2)
        return(dataDF2)
    }

    ### Extract data to retun
    if (!extractGenes) {
        columnsToExtract <-
            c(
                'iso_ref',
                'gene_ref',
                'isoform_id',
                'gene_id',
                'gene_name',
                'condition_1',
                'condition_2',
                'iso_significant',
                'IF1',
                'IF2',
                'dIF',
                'isoform_switch_q_value',
                'switchConsequencesGene'
            )
        isoResTest <-
            any(!is.na(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
            ))
        if (isoResTest) {
            dataDF <- unique(switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value <
                    alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))])
        } else {
            dataDF <- unique(switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$gene_switch_q_value <
                    alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            na.omit(match(
                columnsToExtract,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))])
        }
        if (nrow(dataDF) == 0) {
            stop('No significant switching isoforms were found.')
        }

        if (filterForConsequences) {
            dataDF <- dataDF[which(dataDF$switchConsequencesGene), ]
        }
        if (nrow(dataDF) == 0) {
            stop(
                'No significant switching isoforms with consequences were found.'
            )
        }

        ### Sort data
        dataDF$comparison <-
            paste(dataDF$condition_1, '_vs_', dataDF$condition_2, sep = '')
        if (is.na(sortByQvals)) {
            dataDF2 <- dataDF
            rownames(dataDF2) <- NULL
        } else if (sortByQvals) {
            dataDF2 <-
                dataDF[sort.list(
                    dataDF$isoform_switch_q_value,
                    decreasing = FALSE
                ), ]
            rownames(dataDF2) <- NULL
        } else {
            dataDF2 <- dataDF[sort.list(abs(dataDF$dIF), decreasing = TRUE), ]
            rownames(dataDF2) <- NULL
        }

        ### Reduce to the number wanted
        if (!is.na(n)) {
            if (inEachComparison) {
                dataDF2$comparison <-
                    paste(dataDF2$condition_1,
                          '_vs_',
                          dataDF2$condition_2,
                          sep = '')
            } else {
                dataDF2$comparison <- 'AllCombined'
            }

            dataDF2 <-
                plyr::ddply(
                    .data = dataDF2,
                    .variables = 'comparison',
                    .inform = TRUE,
                    .fun = function(aDF) {
                        if ( n > nrow(dataDF2) ) {
                            if (filterForConsequences) {
                                warning(paste(
                                    'Less than',n ,'genes with significant',
                                    'switches and consequences were found.',
                                    'Returning those.'
                                ))
                            } else {
                                warning(paste(
                                    'Less than',n ,'genes genes with significant',
                                    'switches were found. Returning those.'
                                ))
                            }
                            n2 <- nrow(dataDF)
                        } else {
                            n2 <- n
                        }

                        return(aDF[1:n2, ])
                    }
                )

            dataDF2$comparison <- NULL
        }

        dataDF2$IF1 <- round(dataDF2$IF1, digits = 3)
        dataDF2$IF2 <- round(dataDF2$IF2, digits = 3)
        dataDF2$dIF <- round(dataDF2$dIF, digits = 3)

        dataDF2$Rank <- 1:nrow(dataDF2)

        return(dataDF2)
    }

}