R/analyze_switch_consequences.R

Defines functions extractGenomeWideAnalysis extractConsequenceGenomeWide extractConsequenceEnrichmentComparison extractConsequenceEnrichment extractConsequenceSummary compareAnnotationOfTwoIsoforms analyzeSwitchConsequences

Documented in analyzeSwitchConsequences extractConsequenceEnrichment extractConsequenceEnrichmentComparison extractConsequenceGenomeWide extractConsequenceSummary extractGenomeWideAnalysis

### For analyzing consequences
analyzeSwitchConsequences <- function(
    switchAnalyzeRlist,
    consequencesToAnalyze = c(
        'intron_retention',
        'coding_potential',
        'ORF_seq_similarity',
        'NMD_status',
        'domains_identified',
        'domain_isotype',
        'IDR_identified',
        'IDR_type',
        'signal_peptide_identified'
    ),
    alpha = 0.05,
    dIFcutoff = 0.1,
    onlySigIsoforms = FALSE,
    ntCutoff = 50,
    ntFracCutoff = NULL,
    ntJCsimCutoff = 0.8,
    AaCutoff = 10,
    AaFracCutoff = 0.8,
    AaJCsimCutoff = 0.9,
    removeNonConseqSwitches = TRUE,
    showProgress = TRUE,
    quiet = FALSE
) {
    ### Check input
    if (TRUE) {
        # check switchAnalyzeRlist
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(
                'The object supplied to \'switchAnalyzeRlist\' is not a \'switchAnalyzeRlist\''
            )
        }

        if (alpha < 0 |
            alpha > 1) {
            warning('The alpha parameter should usually be between 0 and 1 ([0,1]).')
        }
        if (alpha > 0.05) {
            warning(
                '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'
            )
        }

        # test wether switching have been analyzed
        if (!any(!is.na(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value
        ))) {
            stop(
                'The analsis of isoform switching must be performed before functional consequences can be analyzed. Please run any of the isoformSwitchTest*() functions and try again.'
            )
        }

        acceptedTypes <- c(
            # Transcript
            'tss',
            'tts',
            'last_exon',
            'isoform_length',
            'exon_number',
            'intron_structure',
            'intron_retention',
            'isoform_class_code',
            # cpat
            'coding_potential',
            # ORF
            'ORF_genomic',
            'ORF_length',
            '5_utr_length',
            '3_utr_length',
            # seq similarity
            'isoform_seq_similarity',
            'ORF_seq_similarity',
            '5_utr_seq_similarity',
            '3_utr_seq_similarity',
            # ORF
            'NMD_status',
            # pfam
            'domains_identified',
            'genomic_domain_position',
            'domain_length',
            'domain_isotype',

            # SignalIP
            'signal_peptide_identified',

            # IDR
            'IDR_identified',
            'IDR_length',
            'IDR_type',

            # sub cell
            'sub_cell_location',
            'sub_cell_shift_to_cell_membrane',
            'sub_cell_shift_to_cytoplasm',
            'sub_cell_shift_to_nucleus',
            'sub_cell_shift_to_Extracellular',

            # topology
            'isoform_topology',
            'extracellular_region_count',
            'intracellular_region_count',
            'extracellular_region_length',
            'intracellular_region_length'
        )

        if (!all(consequencesToAnalyze %in% c('all', acceptedTypes))) {
            stop(
                paste(
                    'The argument(s) supplied to \'typeOfconsequence\' are not accepted.',
                    'Please see ?analyzeSwitchConsequences for description of which strings are allowed.',
                    'The problem is:',
                    paste(setdiff(
                        consequencesToAnalyze , c('all', acceptedTypes)
                    ), collapse = ', '),
                    sep = ' '
                )
            )
        }

        if ('all' %in% consequencesToAnalyze) {
            consequencesToAnalyze <- acceptedTypes
        }

        ## Test whether annotation is advailable
        if ('intron_retention'  %in% consequencesToAnalyze) {
            if (is.null(switchAnalyzeRlist$intronRetentionAnalysis) & is.null( switchAnalyzeRlist$AlternativeSplicingAnalysis)) {
                stop(
                    'To test for intron retention alternative splicing must first be classified. Please run analyzeIntronRetention() and try again.'
                )
            }
        }

        if (grepl('cufflinks', switchAnalyzeRlist$sourceId)) {
            if ('isoform_class_code'  %in% consequencesToAnalyze) {
                if (!'class_code' %in%
                    colnames(switchAnalyzeRlist$isoformFeatures)
                ) {
                    stop(
                        'The switchAnalyzeRlist does not contail the calss_code information'
                    )
                }
            }
        } else {
            consequencesToAnalyze <-
                setdiff(consequencesToAnalyze, 'isoform_class_code')
        }

        if (any(
            c(
                'ORF_length',
                '5_utr_length',
                '3_utr_length',
                'ORF_seq_similarity',
                '5_utr_seq_similarity',
                '3_utr_seq_similarity',
                'domains_identified',
                'genomic_domain_position',
                'domain_length',
                'signal_peptide_identified',
                'IDR_identified',
                'IDR_length',
                'IDR_type',
                'sub_cell_location'
            ) %in% consequencesToAnalyze
        )) {
            if (!'PTC' %in% colnames(switchAnalyzeRlist$isoformFeatures)) {
                stop(
                    'To test differences in ORF or any annotation derived from these, ORF must be annotated. Please run \'addORFfromGTF()\' (and if nessesary \'analyzeNovelIsoformORF()\') and try again'
                )

            }

            if( ! is.null(switchAnalyzeRlist$orfAnalysis$orf_origin ) ) {
                if ( any( switchAnalyzeRlist$orfAnalysis$orf_origin == 'not_annotated_yet' )) {
                    stop('Some ORFs have not been annotated yet. Please return to the analyzeNovelIsoformORF() step and start again.')
                }
            }
        }
        if ('coding_potential'  %in% consequencesToAnalyze) {
            if (
                !'codingPotential' %in%
                colnames(switchAnalyzeRlist$isoformFeatures))
            {
                stop(
                    'To test differences in coding_potential, the result of the CPAT analysis must be advailable. Please run analyzeCPAT() or analyzeCPC2 and try again.'
                )
            }
        }
        if (any(
            c(
                'domains_identified',
                'domain_length',
                'genomic_domain_position'
            )  %in% consequencesToAnalyze
        )) {
            if (is.null(switchAnalyzeRlist$domainAnalysis)) {
                stop(
                    'To test differences in protein domains, the result of the Pfam analysis must be advailable. Please run analyzePFAM() and try again.'
                )
            }
        }
        if ('signal_peptide_identified'  %in% consequencesToAnalyze) {
            if (is.null(switchAnalyzeRlist$signalPeptideAnalysis)) {
                stop(
                    'To test differences in signal peptides, the result of the SignalP analysis must be advailable. Please run analyzeSignalP() and try again.'
                )
            }
        }
        if ( any(c('IDR_identified','IDR_type')  %in% consequencesToAnalyze)) {
            if (is.null(switchAnalyzeRlist$idrAnalysis)) {
                stop(
                    'To test differences in IDR, the result of the NetSurfP2 analysis must be advailable. Please run analyzeNetSurfP2() and try again,'
                )
            }
        }
        if( 'IDR_type' %in% consequencesToAnalyze ) {
            if( ! 'idr_type' %in% colnames(switchAnalyzeRlist$idrAnalysis) ) {
                stop('To analyse IDR_type the IDR analysis must have been done using IUPred2A and imported with the analyzeIUPred2A() function.')
            }
        }

        if (!is.numeric(ntCutoff)) {
            stop('The \'ntCutoff\' arugment must be an numeric')
        }
        if (ntCutoff <= 0) {
            stop('The \'ntCutoff\' arugment must be an numeric > 0')
        }

        if (!is.null(ntFracCutoff)) {
            if (ntFracCutoff <= 0 | ntFracCutoff > 1) {
                stop(
                    'The \'ntFracCutoff\' arugment must be a numeric in the interval (0,1]. Use NULL to disable.'
                )
            }
        }

        ### test sequence annotation
        if (any(
            consequencesToAnalyze %in% c(
                'isoform_seq_similarity',
                '5_utr_seq_similarity',
                '3_utr_seq_similarity'
            )
        )) {
            if (!any(names(switchAnalyzeRlist) == 'ntSequence')) {
                stop(
                    'The transcrip nucleotide sequences must be added to the switchAnalyzeRlist before overlap analysis can be performed. These can be added by using the \'extractSequence()\' function.'
                )
            }
        }
        if (any(consequencesToAnalyze %in% c('ORF_seq_similarity'))) {
            if (!any(names(switchAnalyzeRlist) == 'aaSequence')) {
                stop(
                    'The transcrip ORF amino acid sequences must be added to the switchAnalyzeRlist before ORF overlap analysis can be performed. These can be added by using the \'extractSequence()\' function.'
                )
            }
        }

        if ('sub_cell_location'  %in% consequencesToAnalyze) {
            if ( ! 'sub_cell_location' %in% colnames(switchAnalyzeRlist$isoformFeatures)) {
                stop(
                    'Cannot test for differences in sub-cellular location as such results are not annotated. Run analyzeDeepLoc2() and try again.'
                )
            }
        }

        if ('isoform_topology'  %in% consequencesToAnalyze) {
            if ( ! 'topologyAnalysis' %in% names(switchAnalyzeRlist) ) {
                stop(
                    'Cannot test for differences in topology as such results are not annotated. Run analyzeDeepTMHMM() and try again.'
                )
            }
        }

    }

    if (showProgress & !quiet) {
        progressBar <- 'text'
    } else {
        progressBar <- 'none'
    }

    ### Subset to relevant data
    if (TRUE) {
        if (!quiet) {
            message('Step 1 of 4: Extracting genes with isoform switches...')
        }

        ### Extract Iso pairs
        pairwiseIsoComparison <- extractSwitchPairs(
            switchAnalyzeRlist = switchAnalyzeRlist,
            alpha = alpha,
            dIFcutoff = dIFcutoff,
            onlySigIsoforms = onlySigIsoforms
        )

        ### Extract isoform_id pairs
        pairwiseIsoComparisonUniq <-
            unique(pairwiseIsoComparison[,c(
                'isoformUpregulated', 'isoformDownregulated'
            )])
        pairwiseIsoComparisonUniq$comparison <-
            1:nrow(pairwiseIsoComparisonUniq)

        ### Generate size reduced switchAnalyzeRList
        minimumSwitchList <- makeMinimumSwitchList(
            orgSwitchList = switchAnalyzeRlist,
            isoformsToKeep = unique(
                c(
                    pairwiseIsoComparisonUniq$isoformUpregulated,
                    pairwiseIsoComparisonUniq$isoformDownregulated
                )
            ))


    }

    ### Loop over all the the resulting genes and do a all pairwise comparison between up and down.
    if (TRUE) {
        if (!quiet) {
            message(
                paste(
                    'Step 2 of 4: Analyzing',
                    nrow(pairwiseIsoComparisonUniq),
                    'pairwise isoforms comparisons...',
                    sep = ' '
                )
            )
        }

        consequencesOfIsoformSwitching <- plyr::dlply(
            .data = pairwiseIsoComparisonUniq,
            .variables = 'comparison',
            .parallel = FALSE,
            .inform = TRUE,
            .progress = progressBar,
            .fun = function(aDF) {
                # aDF <- pairwiseIsoComparisonUniq[1,]
                compareAnnotationOfTwoIsoforms(
                    switchAnalyzeRlist    = minimumSwitchList,
                    consequencesToAnalyze = consequencesToAnalyze,
                    upIso                 = aDF$isoformUpregulated,
                    downIso               = aDF$isoformDownregulated,
                    ntCutoff              = ntCutoff,
                    ntFracCutoff          = ntFracCutoff,
                    ntJCsimCutoff         = ntJCsimCutoff,
                    AaCutoff              = AaCutoff,
                    AaFracCutoff          = AaFracCutoff,
                    AaJCsimCutoff         = AaJCsimCutoff,
                    addDescription        = TRUE,
                    testInput             = FALSE # already done by this function
                )
            }
        )

        ### Remove to those instances where there where no consequences
        if (removeNonConseqSwitches) {
            ## Remove to those instances where there where no consequences
            consequencesOfIsoformSwitching <-
                consequencesOfIsoformSwitching[which(
                    sapply(consequencesOfIsoformSwitching, function(aDF) {
                        any(aDF$isoformsDifferent)
                    })
                )]
            if (!length(consequencesOfIsoformSwitching)) {
                stop('No isoform switches with the analyzed consequences were found.')
            }
        }
    }

    ### Massage result
    if (TRUE) {
        if (!quiet) {
            message(paste(
                'Step 3 of 4: Massaging isoforms comparisons results...',
                sep = ' '
            ))
        }
        ### Convert from list to df
        consequencesOfIsoformSwitchingDf <-
            myListToDf(consequencesOfIsoformSwitching)

        ### Add the origin info
        consequencesOfIsoformSwitchingDfcomplete <- dplyr::inner_join(
            consequencesOfIsoformSwitchingDf,
            pairwiseIsoComparison,
            by = c('isoformUpregulated', 'isoformDownregulated')
        )
        #### Add additional information
        #consequencesOfIsoformSwitchingDfcomplete <- dplyr::inner_join(
        #    consequencesOfIsoformSwitchingDfcomplete,
        #    switchAnalyzeRlist$isoformFeatures[match(
        #        unique(consequencesOfIsoformSwitchingDfcomplete$gene_ref),
        #        switchAnalyzeRlist$isoformFeatures$gene_ref
        #    ),
        #    c('gene_ref',
        #      'gene_id',
        #      'gene_name',
        #      'condition_1',
        #      'condition_2')],
        #    by = 'gene_ref'
        #)

        ### reorder
        newOrder <- na.omit(match(
            c(
                'gene_ref',
                'gene_id',
                'gene_name',
                'condition_1',
                'condition_2',
                'isoformUpregulated',
                'isoformDownregulated',
                'iso_ref_up',
                'iso_ref_down',
                'featureCompared',
                'isoformsDifferent',
                'switchConsequence'
            ),
            colnames(consequencesOfIsoformSwitchingDfcomplete)
        ))
        consequencesOfIsoformSwitchingDfcomplete <-
            consequencesOfIsoformSwitchingDfcomplete[, newOrder]

        consequencesOfIsoformSwitchingDfcomplete <-
            consequencesOfIsoformSwitchingDfcomplete[order(
                consequencesOfIsoformSwitchingDfcomplete$gene_ref,
                consequencesOfIsoformSwitchingDfcomplete$isoformUpregulated,
                consequencesOfIsoformSwitchingDfcomplete$isoformDownregulated
            ), ]

    }

    ### Add result to switchAnalyzeRlist
    if (TRUE) {
        if (!quiet) {
            message('Step 4 of 4: Preparing output...')
        }
        ### Add full analysis
        switchAnalyzeRlist$switchConsequence <-
            consequencesOfIsoformSwitchingDfcomplete

        # extract indexes of those analyzed
        indexesAnalyzed <-
            which(
                switchAnalyzeRlist$isoformFeatures$gene_ref %in%
                    pairwiseIsoComparison$gene_ref
            )

        # Add indicator of switch consequence to those analyzed
        switchAnalyzeRlist$isoformFeatures$switchConsequencesGene <- NA

        switchAnalyzeRlist$isoformFeatures$switchConsequencesGene[
            indexesAnalyzed
        ] <- switchAnalyzeRlist$isoformFeatures$gene_ref[indexesAnalyzed] %in%
            consequencesOfIsoformSwitchingDfcomplete$gene_ref[which(
                consequencesOfIsoformSwitchingDfcomplete$isoformsDifferent
            )]
    }

    if (!quiet) {
        totalNrGene <-
            extractSwitchSummary(
                switchAnalyzeRlist,
                filterForConsequences = TRUE,
                includeCombined = TRUE
            )
        totalNrGene <- totalNrGene$nrGenes[which(
            totalNrGene$Comparison == 'combined'
        )]
        message(
            paste(
                'Identified',
                totalNrGene,
                'genes with containing isoforms switching with functional consequences...',
                sep = ' '
            )
        )
    }
    return(switchAnalyzeRlist)
}

compareAnnotationOfTwoIsoforms <- function(
    switchAnalyzeRlist,
    consequencesToAnalyze = 'all',
    upIso,
    downIso,
    addDescription = TRUE,
    onlyRepportDifferent = FALSE,
    ntCutoff = 50,
    ntFracCutoff = NULL,
    ntJCsimCutoff = 0.8,
    AaCutoff = 10,
    AaFracCutoff = 0.8,
    AaJCsimCutoff = 0.9,
    testInput = TRUE
) {
    ### Test and prepare input
    if (testInput) {
        # check switchAnalyzeRlist
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(
                'The object supplied to \'switchAnalyzeRlist\' is not a \'switchAnalyzeRlist\''
            )
        }

        # test wether switching have been analyzed
        if (!any(!is.na(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value
        ))) {
            stop(
                'The analsis of isoform switching must be performed before functional consequences can be analyzed. Please run ?detectIsoformSwitching and try again.'
            )
        }

        acceptedTypes <- c(
            # Transcript
            'tss',
            'tts',
            'last_exon',
            'isoform_length',
            'exon_number',
            'intron_structure',
            'intron_retention',
            'isoform_class_code',
            # cpat
            'coding_potential',
            # ORF
            'ORF_genomic',
            'ORF_length',
            '5_utr_length',
            '3_utr_length',
            # seq similarity
            'isoform_seq_similarity',
            'ORF_seq_similarity',
            '5_utr_seq_similarity',
            '3_utr_seq_similarity',
            # ORF
            'NMD_status',
            # pfam
            'domains_identified',
            'genomic_domain_position',
            'domain_length',
            'domain_isotype',
            # SignalIP
            'signal_peptide_identified',
            # IDR
            'IDR_identified',
            'IDR_length',
            'IDR_type',
            # DeepLoc2
            'sub_cell_location',
            'sub_cell_shift_to_cell_membrane',
            'sub_cell_shift_to_cytoplasm',
            'sub_cell_shift_to_nucleus',
            'sub_cell_shift_to_Extracellular',
            # topology
            'isoform_topology',
            'extracellular_region_count',
            'intracellular_region_count',
            'extracellular_region_length',
            'intracellular_region_length'
        )

        if (!all(consequencesToAnalyze %in% c('all', acceptedTypes))) {
            stop(
                paste(
                    'The argument(s) supplied to \'typeOfconsequence\' are not accepted.',
                    'Please see ?analyzeSwitchConsequences for description of which strings are allowed.',
                    'The problem is:',
                    paste(setdiff(
                        consequencesToAnalyze , c('all', acceptedTypes)
                    ), collapse = ', '),
                    sep = ' '
                )
            )
        }

        if ('all' %in% consequencesToAnalyze) {
            consequencesToAnalyze <- acceptedTypes
        }

        ## Test whether annotation is advailable
        if ('intron_retention'  %in% consequencesToAnalyze) {
            if (is.null(switchAnalyzeRlist$intronRetentionAnalysis) & is.null( switchAnalyzeRlist$AlternativeSplicingAnalysis)) {
                stop(
                    'To test for intron retention alternative splicing must first be classified. Please run analyzeIntronRetention() and try again.'
                )
            }
        }

        if (grepl('cufflinks', switchAnalyzeRlist$sourceId)) {
            if ('isoform_class_code'  %in% consequencesToAnalyze) {
                if (!'class_code' %in%
                    colnames(switchAnalyzeRlist$isoformFeatures)
                ) {
                    stop(
                        'The switchAnalyzeRlist does not contail the calss_code information'
                    )
                }
            }
        } else {
            consequencesToAnalyze <-
                setdiff(consequencesToAnalyze, 'isoform_class_code')
        }

        if (any(c('ORF_genomic', 'ORF_length', 'NMD_status') %in%
                consequencesToAnalyze
        )) {
            if (!'PTC' %in% colnames(switchAnalyzeRlist$isoformFeatures)) {
                stop(
                    'To test differences in ORFs or PCT, ORF must be annotated. Please run \'addORFfromGTF()\' (and if nessesary \'analyzeNovelIsoformORF()\') and try again'
                )
            }
        }
        if ('coding_potential'  %in% consequencesToAnalyze) {
            if (!'codingPotential' %in%
                colnames(switchAnalyzeRlist$isoformFeatures)
            ) {
                stop(
                    'To test differences in coding_potential, the result of the CPAT analysis must be advailable. Please run addCPATanalysis() and try again'
                )
            }
        }
        if (any(
            c(
                'domains_identified',
                'domain_length',
                'genomic_domain_position'
            )  %in% consequencesToAnalyze
        )) {
            if (is.null(switchAnalyzeRlist$domainAnalysis)) {
                stop(
                    'To test differences in protein domains, the result of the Pfam analysis must be advailable. Please run addPFAManalysis() and try again'
                )
            }
        }
        if ('signal_peptide_identified'  %in% consequencesToAnalyze) {
            if (is.null(switchAnalyzeRlist$signalPeptideAnalysis)) {
                stop(
                    'To test differences in signal peptides, the result of the SignalP analysis must be advailable. Please run addSignalIPanalysis() and try again'
                )
            }
        }

        if ( any(c(
            'sub_cell_location',
            'sub_cell_shift_to_cell_membrane',
            'sub_cell_shift_to_cytoplasm',
            'sub_cell_shift_to_nucleus',
            'sub_cell_shift_to_Extracellular'
            ) %in% consequencesToAnalyze)) {
            if ( ! 'sub_cell_location' %in% colnames(switchAnalyzeRlist$isoformFeatures)) {
                stop(
                    'Cannot test for differences in sub-cellular location as such results are not annotated'
                )
            }
        }

        if ( any(
            c(
                'isoform_topology',
                'intracellular_region_length',
                'extracellular_region_length',
                'extracellular_region_count',
                'intracellular_region_count'
            )  %in% consequencesToAnalyze
        ) %in% consequencesToAnalyze) {
            if ( ! 'topologyAnalysis' %in% names(switchAnalyzeRlist) ) {
                stop(
                    'Cannot test for differences in topology as such results are not annotated. Run analyzeDeepTMHMM() and try again.'
                )
            }
        }

        if (!is.numeric(ntCutoff)) {
            stop('The ntCutoff arugment must be an numeric')
        }
        if (ntCutoff <= 0) {
            stop('The ntCutoff arugment must be an numeric > 0')
        }

        if (!is.null(ntFracCutoff)) {
            if (ntFracCutoff <= 0 | ntFracCutoff > 1) {
                stop(
                    'The ntFracCutoff arugment must be a numeric in the interval (0,1]. Use NULL to disable.'
                )
            }
        }

        ### test sequence annotation
        if (any(
            consequencesToAnalyze %in% c(
                'isoform_seq_similarity',
                '5_utr_seq_similarity',
                '3_utr_seq_similarity'
            )
        )) {
            if (!any(names(switchAnalyzeRlist) == 'ntSequence')) {
                stop(
                    'The transcrip nucleotide sequences must be added to the switchAnalyzeRlist before overlap analysis can be performed. Please run \'extractSequence\' and try again.'
                )
            }
        }
        if (any(consequencesToAnalyze %in% c('ORF_seq_similarity'))) {
            if (!any(names(switchAnalyzeRlist) == 'aaSequence')) {
                stop(
                    'The transcrip ORF amino acid sequences must be added to the switchAnalyzeRlist before ORF overlap analysis can be performed. Please run \'extractSequence\' and try again.'
                )
            }
        }

    }

    ### Extract and massage data
    if (TRUE) {
        if (!is.null(ntFracCutoff)) {
            fractionFilter <- TRUE
        } else {
            fractionFilter <- FALSE
        }

        isoformsToAnalyze <-  c(upIso, downIso)
        names(isoformsToAnalyze) <- c('up', 'down')

        ### transcript structure
        exonData <-
            switchAnalyzeRlist$exons[which(
                switchAnalyzeRlist$exons$isoform_id %in% isoformsToAnalyze
                ),
                'isoform_id'
            ] # tss, tts, isoform_length, exon_number
        exonDataList <- split(exonData, f = exonData$isoform_id)

        ### transcript annotation
        columnsToExtract <- 'isoform_id'
        if ('isoform_class_code' %in% consequencesToAnalyze) {
            columnsToExtract <- c(columnsToExtract, 'class_code')
        }
        if ('intron_retention'  %in% consequencesToAnalyze) {
            columnsToExtract <- c(columnsToExtract, 'IR')
        }
        if ('NMD_status'         %in% consequencesToAnalyze) {
            columnsToExtract <- c(columnsToExtract, 'PTC')
        }
        if ('coding_potential'   %in% consequencesToAnalyze) {
            columnsToExtract <- c(columnsToExtract, 'codingPotential')
        }
        if (any(c(
            'sub_cell_location',
            'sub_cell_shift_to_cell_membrane',
            'sub_cell_shift_to_cytoplasm',
            'sub_cell_shift_to_nucleus',
            'sub_cell_shift_to_Extracellular' ) %in% consequencesToAnalyze) ) {
            columnsToExtract <- c(columnsToExtract, 'sub_cell_location')
        }

        transcriptData <-
            unique(switchAnalyzeRlist$isoformFeatures[
                which(
                    switchAnalyzeRlist$isoformFeatures$isoform_id %in%
                        isoformsToAnalyze
                ), columnsToExtract,
                drop =FALSE
            ])

        ### ORF data
        columnsToExtract2 <- 'isoform_id'
        if ('ORF_genomic'  %in% consequencesToAnalyze) {
            columnsToExtract2 <-
                c(columnsToExtract2,
                  c('orfStartGenomic', 'orfEndGenomic'))
        }
        if (any(
            c(
                'ORF_length',
                '5_utr_length',
                '3_utr_length',
                'ORF_seq_similarity',
                '5_utr_seq_similarity',
                '3_utr_seq_similarity',
                'domains_identified',
                'genomic_domain_position',
                'domain_length',
                'signal_peptide_identified',
                'IDR_identified',
                'IDR_length',
                'IDR_type',
                'sub_cell_location',
                'sub_cell_shift_to_cell_membrane',
                'sub_cell_shift_to_cytoplasm',
                'sub_cell_shift_to_nucleus',
                'sub_cell_shift_to_Extracellular',
                'soform_topology',
                'intracellular_region_length',
                'extracellular_region_length',
                'extracellular_region_count',
                'intracellular_region_count'
            ) %in% consequencesToAnalyze
        )) {
            columnsToExtract2 <-
                c(
                    columnsToExtract2,
                    c(
                        'orfTransciptStart',
                        'orfTransciptEnd',
                        'orfTransciptLength'
                    )
                )
        }
        if (length(columnsToExtract2) > 1) {
            orfData <- unique(switchAnalyzeRlist$orfAnalysis[
                which(
                    switchAnalyzeRlist$orfAnalysis$isoform_id %in%
                        isoformsToAnalyze
                    ),
                columnsToExtract2
            ])

            transcriptData <-
                dplyr::left_join(
                    transcriptData,
                    orfData,
                    by = 'isoform_id',
                )
        }

        ### intron retention data
        if ('intron_retention' %in% consequencesToAnalyze) {
            if( is.null( switchAnalyzeRlist$AlternativeSplicingAnalysis) ) {
                localIRdata <-
                    switchAnalyzeRlist$intronRetentionAnalysis[which(
                        switchAnalyzeRlist$intronRetentionAnalysis$isoform_id %in%
                            isoformsToAnalyze
                    ), ]
            } else {
                localIRdata <-
                    switchAnalyzeRlist$AlternativeSplicingAnalysis[which(
                        switchAnalyzeRlist$AlternativeSplicingAnalysis$isoform_id %in%
                            isoformsToAnalyze
                    ), ]
            }

            if (nrow(localIRdata) != 2) {
                warning(
                    paste(
                        'There was a problem with the extraction if intron retentions -',
                        'please contact the developers with this example so they can fix it.',
                        'For now they are ignored. The isoforms affected are:',
                        paste(isoformsToAnalyze, collapse = ', '),
                        sep = ' '
                    )
                )
                consequencesToAnalyze <-
                    consequencesToAnalyze[which(
                        !consequencesToAnalyze %in% c('intron_retention')
                    )]
            } else {
                localIRdata$irCoordinats <-
                    paste(localIRdata$IR_genomic_start,
                          localIRdata$IR_genomic_end,
                          sep = ':')
            }
        }

        ### If nessesary extract data to remove
        onPlusStrand <- as.character(strand(exonData)[1]) == '+'
        if( 'wasTrimmed' %in% colnames(switchAnalyzeRlist$orfAnalysis) ) {

            if(onPlusStrand) {
                localTrimmed <- switchAnalyzeRlist$orfAnalysis[
                    which(switchAnalyzeRlist$orfAnalysis$isoform_id %in% isoformsToAnalyze),
                    c('isoform_id','wasTrimmed','trimmedStartGenomic','orfEndGenomic')
                ]
            } else {
                localTrimmed <- switchAnalyzeRlist$orfAnalysis[
                    which(switchAnalyzeRlist$orfAnalysis$isoform_id %in% isoformsToAnalyze),
                    c('isoform_id','wasTrimmed','trimmedStartGenomic','orfStartGenomic')
                    ]
                colnames(localTrimmed) <- c('isoform_id','wasTrimmed','trimmedStartGenomic','orfEndGenomic')
            }

            if(any(localTrimmed$wasTrimmed, na.rm = TRUE)) {
                localTrimmed <- localTrimmed[which(
                    localTrimmed$wasTrimmed
                ),]

                regionToOmmit <- GenomicRanges::reduce(IRanges(localTrimmed$trimmedStartGenomic, localTrimmed$orfEndGenomic))
            }
        }

        ### domain data
        if (any(
            c(
                'domains_identified',
                'genomic_domain_position',
                'domain_length',
                'domain_isotype'
            )  %in% consequencesToAnalyze
        )) {
            domanData <-
                switchAnalyzeRlist$domainAnalysis[which(
                    switchAnalyzeRlist$domainAnalysis$isoform_id %in%
                        isoformsToAnalyze
                ), ]
            domanData$isoform_id <-
                factor(domanData$isoform_id, levels = isoformsToAnalyze)


            colIndex <- na.omit(match(
                c(
                    'hmm_name',
                    'pfamStartGenomic',
                    'pfamEndGenomic',
                    'domain_isotype_simple',
                    'orf_aa_start',
                    'orf_aa_end'
                ),
                colnames(domanData)
            ))

            domanDataSplit <-
                split(
                    domanData[,colIndex],
                    f = domanData$isoform_id
                )

            ### Remove those overlapping trimmed regions
            if( exists('regionToOmmit') ) {
                domanDataSplit <- lapply(domanDataSplit, function(aSet) { # aSet <- domanDataSplit[[2]]
                    if( onPlusStrand ) {
                        aSet[which(
                            ! overlapsAny(
                                IRanges(aSet$pfamStartGenomic, aSet$pfamEndGenomic),
                                regionToOmmit
                            )
                        ),]
                    } else {
                        aSet[which(
                            ! overlapsAny(
                                IRanges(aSet$pfamEndGenomic, aSet$pfamStartGenomic),
                                regionToOmmit
                            )
                        ),]
                    }

                })
            }

            # overwrite if no ORF is detected
            isNAnames <-
                transcriptData$isoform_id[which(
                    is.na(transcriptData$orfTransciptLength)
                )]
            if (length(isNAnames)) {
                isNAindex <- which(names(domanDataSplit) %in% isNAnames)
                domanDataSplit[isNAindex] <-
                    lapply(domanDataSplit[isNAindex], function(aDF) {
                        aDF[0, ]
                    })
            }
        }

        ### IDR data
        if ( any( c('IDR_identified','IDR_type','IDR_length') %in% consequencesToAnalyze ) ) {
            idrData <-
                switchAnalyzeRlist$idrAnalysis[which(
                    switchAnalyzeRlist$idrAnalysis$isoform_id %in%
                        isoformsToAnalyze
                ), ]
            idrData$isoform_id <-
                factor(idrData$isoform_id, levels = isoformsToAnalyze)
            idrDataSplit <-
                split(idrData[, c(
                    'idrStartGenomic',
                    'idrEndGenomic',
                    'orf_aa_start',
                    'orf_aa_end',
                    'idr_type'
                )], f = idrData$isoform_id)

            ### Remove those overlapping trimmed regions
            if( exists('regionToOmmit') ) {
                idrDataSplit <- lapply(idrDataSplit, function(aSet) { # aSet <- idrDataSplit[[2]]
                    if( onPlusStrand ) {
                        aSet[which(
                            ! overlapsAny(
                                IRanges(aSet$idrStartGenomic, aSet$idrEndGenomic),
                                regionToOmmit
                            )
                        ),]
                    } else {
                        aSet[which(
                            ! overlapsAny(
                                IRanges(aSet$idrEndGenomic, aSet$idrStartGenomic),
                                regionToOmmit
                            )
                        ),]
                    }
                })
            }

            ### overwrite if no ORF is detected
            isNAnames <-
                transcriptData$isoform_id[which(
                    is.na(transcriptData$orfTransciptLength)
                )]
            if (length(isNAnames)) {
                isNAindex <- which(names(idrDataSplit) %in% isNAnames)
                idrDataSplit[isNAindex] <-
                    lapply(idrDataSplit[isNAindex], function(aDF) {
                        aDF[0, ]
                    })
            }
        }

        ### peptide data
        if ('signal_peptide_identified'  %in% consequencesToAnalyze) {
            peptideData <-
                switchAnalyzeRlist$signalPeptideAnalysis[which(
                    switchAnalyzeRlist$signalPeptideAnalysis$isoform_id %in%
                        isoformsToAnalyze
                ), ]
            peptideData$isoform_id <-
                factor(peptideData$isoform_id, levels = isoformsToAnalyze)
            peptideDataSplit <-
                split(peptideData, f = peptideData$isoform_id)


            ### Remove those overlapping trimmed regions
            if( exists('regionToOmmit') ) {
                peptideDataSplit <- lapply(peptideDataSplit, function(aSet) { # aSet <- peptideDataSplit[[2]]
                    aSet[which(
                        ! overlapsAny(
                            IRanges(aSet$genomicClevageAfter, aSet$genomicClevageAfter),
                            regionToOmmit
                        )
                    ),]
                })
            }


            # overwrite if no ORF is detected
            isNAnames <-
                transcriptData$isoform_id[which(
                    is.na(transcriptData$orfTransciptLength)
                )]
            if (length(isNAnames)) {
                isNAindex <- which(names(peptideDataSplit) %in% isNAnames)
                peptideDataSplit[isNAindex] <-
                    lapply(peptideDataSplit[isNAindex], function(aDF) {
                        aDF[0, ]
                    })
            }
        }

        ### Topology
        if (any(
            c(
                'isoform_topology',
                'intracellular_region_length',
                'extracellular_region_length',
                'extracellular_region_count',
                'intracellular_region_count'
            )  %in% consequencesToAnalyze
        )) {
            topData <-
                switchAnalyzeRlist$topologyAnalysis[which(
                    switchAnalyzeRlist$topologyAnalysis$isoform_id %in%
                        isoformsToAnalyze
                ), ]
            topData$isoform_id <-
                factor(topData$isoform_id, levels = isoformsToAnalyze)


            colIndex <- na.omit(match(
                c(
                    'region_type',
                    'regionStartGenomic',
                    'regionEndGenomic'
                ),
                colnames(topData)
            ))

            topDataSplit <-
                split(
                    topData[,colIndex],
                    f = topData$isoform_id
                )

            ### Remove those overlapping trimmed regions
            if( exists('regionToOmmit') ) {
                topDataSplit <- lapply(topDataSplit, function(aSet) { # aSet <- topDataSplit[[2]]
                    if( onPlusStrand ) {
                        aSet[which(
                            ! overlapsAny(
                                IRanges(aSet$regionStartGenomic, aSet$regionEndGenomic),
                                regionToOmmit
                            )
                        ),]
                    } else {
                        aSet[which(
                            ! overlapsAny(
                                IRanges(aSet$regionEndGenomic, aSet$regionStartGenomic),
                                regionToOmmit
                            )
                        ),]
                    }

                })
            }

            # overwrite if no ORF is detected
            isNAnames <-
                transcriptData$isoform_id[which(
                    is.na(transcriptData$orfTransciptLength)
                )]
            if (length(isNAnames)) {
                isNAindex <- which(names(topDataSplit) %in% isNAnames)
                topDataSplit[isNAindex] <-
                    lapply(topDataSplit[isNAindex], function(aDF) {
                        aDF[0, ]
                    })
            }

            topDataSplit <- lapply(topDataSplit, function(x) {
                x$regionLength <- abs(x$regionEndGenomic - x$regionStartGenomic)
                return(x)
            })

        }
        # topDataSplit


        ### Extract isoform length
        if (any(
            c(
                'isoform_length',
                '5_utr_length',
                '3_utr_length',
                'ORF_length',
                'isoform_seq_similarity',
                '5_utr_seq_similarity',
                '3_utr_seq_similarity',
                'ORF_seq_similarity'
            ) %in% consequencesToAnalyze
        )) {
            isoform_length <- sapply(exonDataList, function(x)
                sum(width(x)))
            transcriptData$length <-
                isoform_length[match(transcriptData$isoform_id,
                                     names(isoform_length))]
        }

        ### Sequence overlap
        if (any(
            c(
                'isoform_seq_similarity',
                '5_utr_seq_similarity',
                '3_utr_seq_similarity'
            ) %in% consequencesToAnalyze
        )) {
            ntSeq <-
                switchAnalyzeRlist$ntSequence[which(
                    names(switchAnalyzeRlist$ntSequence) %in% c(upIso, downIso)
                )]
            if (length(ntSeq) == 2) {
                upNtSeq   <- ntSeq[upIso]
                downNtSeq <- ntSeq[downIso]
            } else {
                #warning('The amino acid sequence of some of the transcripts to compare is not pressent - \'isoform_seq_similarity\',\'5_utr_seq_similarity\' and \'3_utr_seq_similarity\' will therefor be ignored.')
                consequencesToAnalyze <-
                    consequencesToAnalyze[which(
                        !consequencesToAnalyze %in% c(
                            'isoform_seq_similarity',
                            '5_utr_seq_similarity',
                            '3_utr_seq_similarity'
                        )
                    )]
            }
        }
        if ('ORF_seq_similarity' %in% consequencesToAnalyze) {
            aaSeq <-
                switchAnalyzeRlist$aaSequence[which(
                    names(switchAnalyzeRlist$aaSequence) %in%
                        transcriptData$isoform_id[which(
                            !is.na(transcriptData$orfTransciptLength)
                        )]
                )]
            if (length(aaSeq) == 2) {
                upAAseq   <- aaSeq[upIso]
                downAAseq <- aaSeq[downIso]
            }
        }

        if (length(consequencesToAnalyze) == 0) {
            return(NULL)
        }
        if (nrow(transcriptData) != 2) {
            return(NULL)
        }

    }

    ### Analyze differences
    if (TRUE) {
        ### Make result data.frame
        isoComparison <-
            data.frame(
                isoformUpregulated = upIso,
                isoformDownregulated = downIso,
                featureCompared = consequencesToAnalyze,
                isoformsDifferent = NA
            )
        if (addDescription) {
            isoComparison$switchConsequence <- NA
        }

        ### analyze  isoform differnces
        if ('tss'                       %in% consequencesToAnalyze) {
            localExonData <- unlist(range(exonDataList))

            if (as.character(localExonData@strand[1]) == '+') {
                tssCoordinats <- start(localExonData)

                ### test
                tssDifferent <-
                    abs(tssCoordinats[1] - tssCoordinats[2]) > ntCutoff
                mostUpstream <-
                    names(localExonData)[which.min(tssCoordinats)]
            } else {
                tssCoordinats <- end(localExonData)

                # test
                tssDifferent <-
                    abs(tssCoordinats[1] - tssCoordinats[2]) > ntCutoff
                mostUpstream <-
                    names(localExonData)[which.max(tssCoordinats)]
            }

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'tss')
            isoComparison$isoformsDifferent[localIndex] <-
                tssDifferent

            # add description
            if (tssDifferent & addDescription) {
                switchMoreUpstram <- mostUpstream == upIso

                if (switchMoreUpstram) {
                    isoComparison$switchConsequence[localIndex] <-
                        'Tss more upstream'
                } else {
                    isoComparison$switchConsequence[localIndex] <-
                        'Tss more downstream'
                }
            }
        }

        if ('tts'                       %in% consequencesToAnalyze) {
            localExonData <- unlist(range(exonDataList))

            if (as.character(localExonData@strand[1]) == '+') {
                ttsCoordinats <- end(localExonData)

                # test
                ttsDifferent <-
                    abs(ttsCoordinats[1] - ttsCoordinats[2]) > ntCutoff
                mostDownstream <-
                    names(localExonData)[which.max(ttsDifferent)]
            } else {
                ttsCoordinats <- start(localExonData)

                # test
                ttsDifferent <-
                    abs(ttsCoordinats[1] - ttsCoordinats[2]) > ntCutoff
                mostDownstream <-
                    names(localExonData)[which.min(ttsDifferent)]
            }

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'tts')
            isoComparison$isoformsDifferent[localIndex] <-
                ttsDifferent

            # add description
            if (ttsDifferent & addDescription) {
                switchMoreDownstream <- mostDownstream == upIso

                if (switchMoreDownstream) {
                    isoComparison$switchConsequence[localIndex] <-
                        'Tts more downstream'
                } else {
                    isoComparison$switchConsequence[localIndex] <-
                        'Tts more upstream'
                }
            }
        }

        if ('last_exon'                 %in% consequencesToAnalyze) {
            isPlusStrand <- as.logical(exonData@strand[1] == '+')
            if (isPlusStrand) {
                lastExons <- lapply(exonDataList, function(x)
                    x[length(x), 0])
            } else {
                lastExons <- lapply(exonDataList, function(x)
                    x[1        , 0])
            }

            lastExonDifferent <-
                !overlapsAny(lastExons[[1]], lastExons[[2]])

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'last_exon')
            isoComparison$isoformsDifferent[localIndex] <-
                lastExonDifferent

            # add description
            if (lastExonDifferent & addDescription) {
                if (isPlusStrand) {
                    localEndCoordinats <-
                        sapply(lastExons, function(aGRange)
                            end(aGRange))
                    mostDownstream <-
                        names(localEndCoordinats)[which.max(localEndCoordinats)]
                } else {
                    localEndCoordinats <-
                        sapply(lastExons, function(aGRange)
                            start(aGRange))
                    mostDownstream <-
                        names(localEndCoordinats)[which.min(localEndCoordinats)]
                }

                switchMoreDownstream <- mostDownstream == upIso

                if (switchMoreDownstream) {
                    isoComparison$switchConsequence[localIndex] <-
                        'Last exon more downstream'
                } else {
                    isoComparison$switchConsequence[localIndex] <-
                        'Last exon more upstream'
                }
            }
        }

        if ('isoform_length'            %in% consequencesToAnalyze) {
            differentLength <- abs(diff(isoform_length)) > ntCutoff

            if (fractionFilter) {
                downLength <-
                    isoform_length[which(names(isoform_length) == downIso)]
                fractionDifference <-
                    abs(diff(isoform_length)) / downLength > ntFracCutoff

                differentLength <-
                    differentLength & fractionDifference
            }

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'isoform_length')
            isoComparison$isoformsDifferent[localIndex] <-
                differentLength

            # make description
            if (differentLength & addDescription) {
                lengthGain <-
                    names(isoform_length)[which.max(isoform_length)] == upIso

                if (lengthGain) {
                    isoComparison$switchConsequence[localIndex] <- 'Length gain'
                } else {
                    isoComparison$switchConsequence[localIndex] <- 'Length loss'
                }

            }
        }

        if ('isoform_seq_similarity'    %in% consequencesToAnalyze) {
            localAlignment <-
                Biostrings::pairwiseAlignment(pattern = upNtSeq,
                                  subject = downNtSeq,
                                  type = 'global')

            overlapSize <- min(c(
                nchar(gsub(
                    '-', '', as.character(Biostrings::alignedSubject(localAlignment))
                )),
                nchar(gsub(
                    '-', '', as.character(Biostrings::alignedPattern(localAlignment))
                ))
            ))
            totalWidth <- width(localAlignment@subject@unaligned) +
                width(localAlignment@pattern@unaligned) -
                overlapSize + 1

            jcDist <- overlapSize / totalWidth

            differentIsoformOverlap <-
                jcDist < ntJCsimCutoff & totalWidth - overlapSize > ntCutoff

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'isoform_seq_similarity')
            isoComparison$isoformsDifferent[localIndex] <-
                differentIsoformOverlap

            # make description
            if (differentIsoformOverlap & addDescription) {
                lengthGain <-
                    names(isoform_length)[which.max(isoform_length)] == upIso

                if (lengthGain) {
                    isoComparison$switchConsequence[localIndex] <- 'Length gain'
                } else {
                    isoComparison$switchConsequence[localIndex] <- 'Length loss'
                }

            }
        }

        if ('exon_number'               %in% consequencesToAnalyze) {
            localNrExons <- sapply(exonDataList, length)
            exonDifferent <- localNrExons[1] != localNrExons[2]

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'exon_number')
            isoComparison$isoformsDifferent[localIndex] <-
                exonDifferent

            # make description
            if (exonDifferent & addDescription) {
                switchGainsExons <-
                    names(localNrExons)[which.max(localNrExons)] == upIso

                if (switchGainsExons) {
                    isoComparison$switchConsequence[localIndex] <- 'Exon gain'
                } else {
                    isoComparison$switchConsequence[localIndex] <- 'Exon loss'
                }
            }

        }

        if ('intron_structure'          %in% consequencesToAnalyze) {
            localIntrons <-
                lapply(exonDataList, function(aGR) {
                    gaps(ranges(aGR))
                })

            differentintron_structure <- !any(
                all(localIntrons[[1]] %in% localIntrons[[2]]),
                all(localIntrons[[2]] %in% localIntrons[[1]])
            )

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'intron_structure')
            isoComparison$isoformsDifferent[localIndex] <-
                differentintron_structure
        }

        if ('intron_retention'          %in% consequencesToAnalyze) {
            if (all(!is.na(transcriptData$IR))) {
                differentNrIR <-
                    localIRdata$irCoordinats[1] != localIRdata$irCoordinats[2]

                # Make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'intron_retention')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentNrIR

                # make description
                if (differentNrIR & addDescription) {
                    if (transcriptData$IR[1] == transcriptData$IR[2]) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Intron retention switch'
                    } else if (transcriptData$isoform_id[which.max(
                        transcriptData$IR
                    )] == upIso) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Intron retention gain'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Intron retention loss'
                    }
                }
            }
        }

        if ('isoform_class_code'        %in% consequencesToAnalyze) {
            differentAnnotation <-
                transcriptData$class_code[1] != transcriptData$class_code[2]

            # make repport
            localIndex <-
                which(isoComparison$featureCompared == 'isoform_class_code')
            isoComparison$isoformsDifferent[localIndex] <-
                differentAnnotation
        }

        if ('ORF_genomic'               %in% consequencesToAnalyze) {
            localIndex <- which(isoComparison$featureCompared == 'ORF_genomic')

            # if both have ORFs
            if (all(!is.na(transcriptData$orfStartGenomic))) {
                genomicOrfDifferent <-
                    transcriptData$orfStartGenomic[1] !=
                    transcriptData$orfStartGenomic[2] |
                    transcriptData$orfEndGenomic[1] !=
                    transcriptData$orfEndGenomic[2]

                # make repport
                isoComparison$isoformsDifferent[localIndex] <-
                    genomicOrfDifferent
                # If only 1 ORF
            } else if (sum(!is.na(transcriptData$orfStartGenomic)) == 1) {
                # make repport
                isoComparison$isoformsDifferent[localIndex] <- TRUE

            } else {
                isoComparison$isoformsDifferent[localIndex] <- FALSE
            }
        }

        if ('ORF_length'                %in% consequencesToAnalyze) {
            localIndex <- which(isoComparison$featureCompared == 'ORF_length')
            # if both ORFs are annotated
            if (all(!is.na(transcriptData$orfTransciptLength))) {
                orfDifferent <-
                    abs(diff(transcriptData$orfTransciptLength)) > ntCutoff

                if (fractionFilter) {
                    downLength <-
                        isoform_length[which(names(isoform_length) == downIso)]
                    fractionDifference <-
                        abs(diff(transcriptData$orfTransciptLength)) /
                        downLength > ntFracCutoff

                    orfDifferent <-
                        orfDifferent & fractionDifference
                }

                # make repport
                isoComparison$isoformsDifferent[localIndex] <-
                    orfDifferent

                # make description
                if (orfDifferent & addDescription) {
                    orfGain <-
                        transcriptData$isoform_id[which.max(
                            transcriptData$orfTransciptLength
                        )] == upIso

                    if (orfGain) {
                        isoComparison$switchConsequence[localIndex] <-
                            'ORF is longer'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'ORF is shorter'
                    }
                }
                #If only one ORF is annotated
            } else if (sum(!is.na(transcriptData$orfTransciptLength)) == 1) {
                # make repport
                isoComparison$isoformsDifferent[localIndex] <- TRUE

                if (addDescription) {
                    orfLoss <-
                        is.na(transcriptData$orfTransciptLength[which(
                            transcriptData$isoform_id == upIso
                        )])

                    if (orfLoss) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Complete ORF loss'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Complete ORF gain'
                    }
                }
            } else {
                isoComparison$isoformsDifferent[localIndex] <- FALSE
            }
        }

        if ('ORF_seq_similarity'        %in% consequencesToAnalyze) {
            localIndex <-
                which(isoComparison$featureCompared == 'ORF_seq_similarity')
            # if both have ORFs
            if (all(!is.na(transcriptData$orfTransciptLength))) {
                if (length(aaSeq) == 2) {
                    localAlignment <-
                        Biostrings::pairwiseAlignment(
                            pattern = upAAseq,
                            subject = downAAseq,
                            type = 'global'
                        )

                    overlapSize <- min(c(nchar(
                        gsub(
                            '-',
                            '',
                            as.character(Biostrings::alignedSubject(localAlignment))
                        )
                    ),
                        nchar(
                            gsub(
                                '-',
                                '',
                                as.character(Biostrings::alignedPattern(localAlignment))
                            )
                        ))
                    )
                    totalWidth <-
                        width(localAlignment@subject@unaligned) +
                        width(localAlignment@pattern@unaligned) -
                        overlapSize + 1

                    jcDist <- overlapSize / totalWidth

                    differentORFoverlap <-
                        jcDist < AaJCsimCutoff &
                        (totalWidth - overlapSize) > AaCutoff

                    # make repport
                    isoComparison$isoformsDifferent[localIndex] <-
                        differentORFoverlap

                    # make description
                    if (differentORFoverlap & addDescription) {
                        lengthGain <-
                            transcriptData$isoform_id[which.max(
                                transcriptData$orfTransciptLength
                            )] == upIso

                        if (lengthGain) {
                            isoComparison$switchConsequence[localIndex] <-
                                'ORF is longer'
                        } else {
                            isoComparison$switchConsequence[localIndex] <-
                                'ORF is shorter'
                        }

                    }
                }
                # If only 1 ORF is annotated
            } else if (sum(!is.na(transcriptData$orfTransciptLength)) == 1) {
                # make repport
                isoComparison$isoformsDifferent[localIndex] <- TRUE

                if (addDescription) {
                    orfLoss <-
                        is.na(transcriptData$orfTransciptLength[which(
                            transcriptData$isoform_id == upIso
                        )])

                    if (orfLoss) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Complete ORF loss'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Complete ORF gain'
                    }
                }
            } else {
                isoComparison$isoformsDifferent[localIndex] <- FALSE
            }
        }

        if ('5_utr_length'              %in% consequencesToAnalyze) {
            if (all(!is.na(transcriptData$orfTransciptStart))) {
                utr5Different <-
                    abs(diff(transcriptData$orfTransciptStart)) > ntCutoff

                if (fractionFilter) {
                    downLength <-
                        isoform_length[which(names(isoform_length) == downIso)]
                    fractionDifference <-
                        abs(diff(transcriptData$orfTransciptStart)) /
                        downLength > ntFracCutoff

                    utr5Different <-
                        utr5Different & fractionDifference
                }


                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == '5_utr_length')
                isoComparison$isoformsDifferent[localIndex] <-
                    utr5Different

                # make description
                if (utr5Different & addDescription) {
                    utr5Gain <-
                        transcriptData$isoform_id[which.max(
                            transcriptData$orfTransciptStart
                        )] == upIso

                    if (utr5Gain) {
                        isoComparison$switchConsequence[localIndex] <-
                            '5UTR is longer'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            '5UTR is shorter'
                    }
                }
            }
        }

        if ('5_utr_seq_similarity'      %in% consequencesToAnalyze) {
            if (all(!is.na(transcriptData$orfTransciptStart))) {
                localUpNt   <-
                    Biostrings::subseq(upNtSeq,
                           1,
                           transcriptData$orfTransciptStart[which(
                               transcriptData$isoform_id == upIso
                           )] - 1)
                localDownNt <-
                    Biostrings::subseq(downNtSeq,
                           1,
                           transcriptData$orfTransciptStart[which(
                               transcriptData$isoform_id == downIso
                           )] -1)

                # if one of them have no UTR
                if (width(localUpNt) > 0 &
                    width(localDownNt) > 0) {
                    localAlignment <-
                        Biostrings::pairwiseAlignment(
                            pattern = localUpNt,
                            subject = localDownNt,
                            type = 'overlap'
                        )

                    overlapSize <- min(c(nchar(
                        gsub(
                            '-',
                            '',
                            as.character(Biostrings::alignedSubject(localAlignment))
                        )
                    ),
                    nchar(
                        gsub(
                            '-',
                            '',
                            as.character(Biostrings::alignedPattern(localAlignment))
                        )
                    )))
                    totalWidth <-
                        width(localAlignment@subject@unaligned) +
                        width(localAlignment@pattern@unaligned) -
                        overlapSize + 1

                    jcDist <- overlapSize / totalWidth
                } else {
                    overlapSize <- 0
                    totalWidth <-
                        abs(width(localUpNt) - width(localDownNt))
                    jcDist <- 0
                }

                differenttUTRoverlap <-
                    jcDist < ntJCsimCutoff & totalWidth - overlapSize > ntCutoff

                # make repport
                localIndex <-
                    which(
                        isoComparison$featureCompared == '5_utr_seq_similarity'
                    )
                isoComparison$isoformsDifferent[localIndex] <-
                    differenttUTRoverlap

                # make description
                if (differenttUTRoverlap & addDescription) {
                    utr5Gain <-
                        transcriptData$isoform_id[which.max(
                            transcriptData$orfTransciptStart
                        )] == upIso

                    if (utr5Gain) {
                        isoComparison$switchConsequence[localIndex] <-
                            '5UTR is longer'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            '5UTR is shorter'
                    }

                }
            }
        }

        if ('3_utr_length'              %in% consequencesToAnalyze) {
            if (all(!is.na(transcriptData$orfTransciptEnd))) {
                transcriptData$utr3length <-
                    transcriptData$length - (transcriptData$orfTransciptEnd + 1)

                utr3Different <-
                    abs(diff(transcriptData$utr3length)) > ntCutoff

                if (fractionFilter) {
                    downLength <-
                        isoform_length[which(names(isoform_length) == downIso)]
                    fractionDifference <-
                        abs(diff(transcriptData$utr3length)) /
                        downLength > ntFracCutoff

                    utr3Different <-
                        utr3Different & fractionDifference
                }

                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == '3_utr_length')
                isoComparison$isoformsDifferent[localIndex] <-
                    utr3Different

                # make description
                if (utr3Different & addDescription) {
                    utr3Gain <-
                        transcriptData$isoform_id[which.max(
                            transcriptData$utr3length
                        )] == upIso

                    if (utr3Gain) {
                        isoComparison$switchConsequence[localIndex] <-
                            '3UTR is longer'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            '3UTR is shorter'
                    }
                }
            }
        }

        if ('3_utr_seq_similarity'      %in% consequencesToAnalyze) {
            if (all(!is.na(transcriptData$orfTransciptEnd))) {
                transcriptData$utr3length <-
                    transcriptData$length - (transcriptData$orfTransciptEnd + 1)

                up3UTRstart   <-
                    transcriptData$orfTransciptEnd[which(
                        transcriptData$isoform_id == upIso
                    )] + 4 # +4 because the stop codon is not included in the ORF
                down3UTRstart <-
                    transcriptData$orfTransciptEnd[which(
                        transcriptData$isoform_id == downIso
                    )] + 4 # +4 because the stop codon is not included in the ORF

                # if 3UTR somehow is annotated as longer set it to the length +1 (then the sequence will be length 0)
                if (up3UTRstart   > width(upNtSeq)   + 1) {
                    up3UTRstart   <- width(upNtSeq)  + 1
                }
                if (down3UTRstart > width(downNtSeq) + 1) {
                    down3UTRstart <- width(downNtSeq) + 1
                }

                localUpNt   <-
                    subseq(upNtSeq,   up3UTRstart,   width(upNtSeq))
                localDownNt <-
                    subseq(downNtSeq, down3UTRstart, width(downNtSeq))

                ### if one of them have no UTR
                if (width(localUpNt) > 0 &
                    width(localDownNt) > 0) {
                    localAlignment <-
                        Biostrings::pairwiseAlignment(
                            pattern = localUpNt,
                            subject = localDownNt,
                            type = 'overlap'
                        )

                    overlapSize <- min(c(nchar(
                        gsub(
                            '-',
                            '',
                            as.character(Biostrings::alignedSubject(localAlignment))
                        )
                    ),
                    nchar(
                        gsub(
                            '-',
                            '',
                            as.character(Biostrings::alignedPattern(localAlignment))
                        )
                    )))
                    totalWidth <-
                        width(localAlignment@subject@unaligned) +
                        width(localAlignment@pattern@unaligned) -
                        overlapSize + 1

                    jcDist <- overlapSize / totalWidth
                } else {
                    overlapSize <- 0
                    totalWidth <-
                        abs(width(localUpNt) - width(localDownNt))
                    jcDist <- 0
                }

                different3UTRoverlap <-
                    jcDist < ntJCsimCutoff &
                    totalWidth - overlapSize > ntCutoff

                # make repport
                localIndex <-
                    which(
                        isoComparison$featureCompared == '3_utr_seq_similarity'
                    )
                isoComparison$isoformsDifferent[localIndex] <-
                    different3UTRoverlap

                # make description
                if (different3UTRoverlap & addDescription) {
                    utr3Gain <-
                        transcriptData$isoform_id[which.max(
                            transcriptData$utr3length
                        )] == upIso

                    if (utr3Gain) {
                        isoComparison$switchConsequence[localIndex] <-
                            '3UTR is longer'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            '3UTR is shorter'
                    }

                }
            }
        }

        if ('NMD_status'                      %in% consequencesToAnalyze) {
            if (all(!is.na(transcriptData$PTC))) {
                ptcDifferent <- transcriptData$PTC[1] != transcriptData$PTC[2]

                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'NMD_status')
                isoComparison$isoformsDifferent[localIndex] <-
                    ptcDifferent

                # make description
                if (ptcDifferent & addDescription) {
                    upIsSensitive <-
                        transcriptData$PTC[which(
                            transcriptData$isoform_id == upIso
                        )]

                    if (upIsSensitive) {
                        isoComparison$switchConsequence[localIndex] <-
                            'NMD sensitive'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'NMD insensitive'
                    }
                }
            }
        }

        if ('coding_potential'                %in% consequencesToAnalyze) {
            if (all(!is.na(transcriptData$codingPotential))) {
                cpDifferent <-
                    transcriptData$codingPotential[1] !=
                    transcriptData$codingPotential[2]

                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'coding_potential')
                isoComparison$isoformsDifferent[localIndex] <-
                    cpDifferent

                # make description
                if (cpDifferent & addDescription) {
                    upIsCoding <-
                        transcriptData$codingPotential[which(
                            transcriptData$isoform_id == upIso
                        )]

                    if (upIsCoding) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Transcript is coding'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Transcript is Noncoding'
                    }
                }
            }
        }

        if ('domains_identified'              %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {
                domianNames <- lapply(domanDataSplit, function(x)
                    x$hmm_name)

                #t1 <- table(domianNames[[1]])
                #t2 <- table(domianNames[[2]])

                t1 <- rle(domianNames[[1]])
                t2 <- rle(domianNames[[2]])

                differentDomains <- !identical(t1, t2)
                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'domains_identified')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentDomains

                if (differentDomains & addDescription) {
                    if (sum(t1$lengths) != sum(t2$lengths)) {
                        nrDomains <- sapply(domanDataSplit, nrow)
                        upHasMoreDomains <-
                            names(nrDomains)[which.max(nrDomains)] == upIso

                        if (upHasMoreDomains) {
                            isoComparison$switchConsequence[localIndex] <-
                                'Domain gain'
                        } else {
                            isoComparison$switchConsequence[localIndex] <-
                                'Domain loss'
                        }
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Domain switch'
                    }
                }
            }
        }

        if ('genomic_domain_position'         %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {
                localDomanDataSplit <-
                    lapply(domanDataSplit, function(aDF) {
                        aDF[
                            sort.list(aDF$pfamStartGenomic),
                            c('hmm_name','pfamStartGenomic','pfamEndGenomic')
                        ]
                    })
                genomicDomainDifferent <-
                    !identical(
                        localDomanDataSplit[[1]],
                        localDomanDataSplit[[2]]
                    )

                # make repport
                localIndex <-
                    which(
                        isoComparison$featureCompared ==
                            'genomic_domain_position'
                    )
                isoComparison$isoformsDifferent[localIndex] <-
                    genomicDomainDifferent
            }
        }

        if ('domain_length'                   %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {
                domanDataSplit <- lapply(
                    domanDataSplit,
                    function(x) {
                        x$length <- x$orf_aa_end - x$orf_aa_start + 1
                        return(x)
                    }
                )

                ### Analyze overlap
                nDom <- sapply(domanDataSplit, nrow)
                if( all(nDom) ) {
                    domainRanges <- lapply(
                        domanDataSplit,
                        function(x) {
                            if( x$pfamStartGenomic[1] < x$pfamEndGenomic[1]) {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$pfamStartGenomic,
                                        end   = x$pfamEndGenomic
                                    ),
                                    hmm_name  = x$hmm_name,
                                    orfLength = x$length
                                )
                            } else {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$pfamEndGenomic,
                                        end   = x$pfamStartGenomic
                                    ),
                                    hmm_name  = x$hmm_name,
                                    orfLength = x$length
                                )
                            }

                        }
                    )

                    ### Calculate overlap
                    localOverlap <- findOverlaps(
                        query   = domainRanges[[ upIso   ]],
                        subject = domainRanges[[ downIso ]]
                    )

                    if(length(localOverlap)) {

                        localOverlapDf <- as.data.frame(localOverlap)
                        localOverlapDf$upName <- domainRanges[[ upIso   ]]$hmm_name[queryHits   (localOverlap)]
                        localOverlapDf$dnName <- domainRanges[[ downIso ]]$hmm_name[subjectHits (localOverlap)]
                        localOverlapDf$upLength <- domainRanges[[ upIso   ]]$orfLength[queryHits   (localOverlap)]
                        localOverlapDf$dnLength <- domainRanges[[ downIso ]]$orfLength[subjectHits (localOverlap)]

                        ### Subset to same domain
                        localOverlapDf <- localOverlapDf[which(
                            localOverlapDf$upName == localOverlapDf$dnName
                        ),]

                        localOverlapDf$lengthDiff <-
                            abs(localOverlapDf$upLength - localOverlapDf$dnLength) > AaCutoff

                        if (!is.null(AaFracCutoff)) {
                            localOverlapDf$minLength <- apply(localOverlapDf[,c('upLength','dnLength')], 1, min)
                            localOverlapDf$maxLength <- apply(localOverlapDf[,c('upLength','dnLength')], 1, max)

                            localOverlapDf$lengthDiff <-
                                localOverlapDf$minLength / localOverlapDf$maxLength < AaFracCutoff & localOverlapDf$lengthDiff

                        }

                        localOverlapDfDiff <- localOverlapDf[which(
                            localOverlapDf$lengthDiff
                        ),]

                        differentDomainLength <- nrow(localOverlapDfDiff) > 0

                        localIndex <-
                            which(isoComparison$featureCompared == 'domain_length')
                        isoComparison$isoformsDifferent[localIndex] <-
                            differentDomainLength

                    } else {
                        differentDomainLength <- FALSE
                    }

                } else {
                    differentDomainLength <- FALSE
                }

                ### Repport if any difference
                if (differentDomainLength & addDescription) {

                    localOverlapDfDiff$maxIsUp <- localOverlapDfDiff$maxLength == localOverlapDfDiff$upLength

                    ### Deside consequence
                    if(
                        all( c(TRUE, FALSE) %in% localOverlapDfDiff$maxIsUp )
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            #'IDR length gain and loss'
                            'Mixed Domain length differences'
                    } else if(
                        all(localOverlapDfDiff$maxIsUp)
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Domain length gain'
                    } else if(
                        all( ! localOverlapDfDiff$maxIsUp )
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Domain length loss'
                    } else {
                        stop('Something with Domain length analysis went wrong')
                    }

                }

            }
        }

        if ('domain_isotype'                  %in% consequencesToAnalyze) {
            if (
                sum(!is.na(transcriptData$orfTransciptLength)) > 0
            ) {
                domanDataSplit <- lapply(
                    domanDataSplit,
                    function(x) {
                        x$length <- x$orf_aa_end - x$orf_aa_start + 1
                        return(x)
                    }
                )

                ### Analyze structure
                nDom <- sapply(domanDataSplit, nrow)
                if( all(nDom) ) {
                    domainRanges <- lapply(
                        domanDataSplit,
                        function(x) {
                            if( x$pfamStartGenomic[1] < x$pfamEndGenomic[1]) {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$pfamStartGenomic,
                                        end   = x$pfamEndGenomic
                                    ),
                                    hmm_name  = x$hmm_name,
                                    domain_isotype  = x$domain_isotype_simple
                                )
                            } else {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$pfamEndGenomic,
                                        end   = x$pfamStartGenomic
                                    ),
                                    hmm_name  = x$hmm_name,
                                    domain_isotype  = x$domain_isotype_simple
                                )
                            }

                        }
                    )

                    ### Calculate overlap
                    localOverlap <- findOverlaps(
                        query   = domainRanges[[ upIso   ]],
                        subject = domainRanges[[ downIso ]]
                    )

                    if(length(localOverlap)) {

                        localOverlapDf <- as.data.frame(localOverlap)
                        localOverlapDf$upName <- domainRanges[[ upIso   ]]$hmm_name[queryHits   (localOverlap)]
                        localOverlapDf$dnName <- domainRanges[[ downIso ]]$hmm_name[subjectHits (localOverlap)]
                        localOverlapDf$upStructure <- domainRanges[[ upIso   ]]$domain_isotype[queryHits   (localOverlap)]
                        localOverlapDf$dnStructure <- domainRanges[[ downIso ]]$domain_isotype[subjectHits (localOverlap)]

                        ### Subset to same domain
                        localOverlapDf <- localOverlapDf[which(
                            localOverlapDf$upName == localOverlapDf$dnName
                        ),]

                        if(nrow(localOverlapDf)) {
                            localOverlapDf$upCombined <- stringr::str_c(localOverlapDf$upName, '_', localOverlapDf$upStructure)
                            localOverlapDf$dnCombined <- stringr::str_c(localOverlapDf$dnName, '_', localOverlapDf$dnStructure)

                            differentDomainStructure <- any( localOverlapDf$upCombined != localOverlapDf$dnCombined)
                        } else {
                            differentDomainStructure <- FALSE
                        }
                    } else {
                        differentDomainStructure <- FALSE
                    }

                } else {
                    differentDomainStructure <- FALSE
                }

                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'domain_isotype')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentDomainStructure

                ### Report if any difference
                if (differentDomainStructure & addDescription) {

                    upCount <- sum(localOverlapDf$upStructure != 'Reference')
                    dnCount <- sum(localOverlapDf$dnStructure != 'Reference')

                    ### Deside consequence
                    if(
                        all( c(upCount, dnCount) > 0 ) & upCount == dnCount
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Mixed domain isotype changes'
                    } else if(
                        upCount > dnCount
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Domain non-reference isotype gain'
                    } else if(
                        dnCount > upCount
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Domain non-reference isotype loss'
                    }

                }
            }
        }

        if ('IDR_identified'                  %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {

                nIdr <- sapply(idrDataSplit, nrow)
                if( all(nIdr) ) {
                    idrRanges <- lapply(
                        idrDataSplit,
                        function(x) {
                            if( x$idrStartGenomic[1] < x$idrEndGenomic[1]) {
                                IRanges(
                                    start = x$idrStartGenomic,
                                    end   = x$idrEndGenomic
                                )
                            } else {
                                IRanges(
                                    start = x$idrEndGenomic,
                                    end   = x$idrStartGenomic
                                )
                            }

                        }
                    )

                    ### Calculate overlap
                    overlap1 <- overlapsAny(idrRanges[[1]], idrRanges[[2]])
                    overlap2 <- overlapsAny(idrRanges[[2]], idrRanges[[1]])
                    #overlap1 <- grangesFracOverlap(idrRanges[[1]], idrRanges[[2]])$fracOverlap >= maxIdrFracOverlap
                    #overlap2 <- grangesFracOverlap(idrRanges[[2]], idrRanges[[1]])$fracOverlap >= maxIdrFracOverlap

                } else if( nIdr[1] == 0 & nIdr[2] == 0 ) {
                    overlap1 <- NA
                    overlap2 <- NA
                } else if( nIdr[1] == 0 ) {
                    overlap1 <- TRUE
                    overlap2 <- FALSE
                } else {
                    overlap1 <- FALSE
                    overlap2 <- TRUE
                }

                ### Test overlap
                differentIdr <- any( c(
                    ! overlap1,
                    ! overlap2
                ), na.rm = TRUE)

                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'IDR_identified')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentIdr

                if (differentIdr & addDescription) {

                    if( any(!overlap1) & any(!overlap2) ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'IDR switch'
                    } else {
                        if( any(!overlap1) ) {
                            theDiff <- 1
                        } else {
                            theDiff <- 2
                        }

                        upHasMoreIDR <- names(idrDataSplit)[theDiff] == upIso

                        if (upHasMoreIDR) {
                            isoComparison$switchConsequence[localIndex] <-
                                'IDR gain'
                        } else {
                            isoComparison$switchConsequence[localIndex] <-
                                'IDR loss'
                        }
                    }
                }
            }
        }

        if ('IDR_type'                        %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {

                nIdr <- sapply(idrDataSplit, nrow)
                if( all(nIdr) ) {
                    idrRanges <- lapply(
                        idrDataSplit,
                        function(x) {
                            if( x$idrStartGenomic[1] < x$idrEndGenomic[1]) {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$idrStartGenomic,
                                        end   = x$idrEndGenomic
                                    ),
                                    type  = x$idr_type
                                )
                            } else {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$idrEndGenomic,
                                        end   = x$idrStartGenomic
                                    ),
                                    type  = x$idr_type
                                )
                            }

                        }
                    )

                    ### Calculate overlap
                    localOverlap <- findOverlaps(
                        query   = idrRanges[[ upIso   ]],
                        subject = idrRanges[[ downIso ]]
                    )

                    if(length(localOverlap)) {
                        localOverlapDf <- as.data.frame(localOverlap)
                        localOverlapDf$upType <- idrRanges[[ upIso   ]]$type[queryHits   (localOverlap)]
                        localOverlapDf$dnType <- idrRanges[[ downIso ]]$type[subjectHits (localOverlap)]

                        differentIdrType <- any( na.omit(
                            localOverlapDf$upType != localOverlapDf$dnType
                        ))

                        localIndex <-
                            which(isoComparison$featureCompared == 'IDR_type')
                        isoComparison$isoformsDifferent[localIndex] <-
                            differentIdrType
                    } else {
                        differentIdrType <- FALSE
                    }
                } else {
                    differentIdrType <- FALSE
                }

                if (differentIdrType & addDescription) {

                    localOverlapDf <- localOverlapDf[which(
                        localOverlapDf$upType != localOverlapDf$dnType
                    ),]

                    localOverlapDf$bindingGain <-
                        localOverlapDf$dnType == 'IDR' & localOverlapDf$upType == 'IDR_w_binding_region'

                    localOverlapDf$bindingLoss <-
                        localOverlapDf$dnType == 'IDR_w_binding_region' & localOverlapDf$upType == 'IDR'

                    ### Deside consequence
                    if(
                        any( localOverlapDf$bindingGain) &
                        any(localOverlapDf$bindingLoss)
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'IDR w binding region switch'
                    } else if(
                          any( localOverlapDf$bindingGain ) &
                        ! any( localOverlapDf$bindingLoss )
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'IDR w binding region gain'
                    } else if(
                        ! any( localOverlapDf$bindingGain ) &
                        any( localOverlapDf$bindingLoss )
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'IDR w binding region loss'
                    } else {
                        stop('Something with idr binding regions went wrong')
                    }

                }
            }
        }

        if ('IDR_length'                      %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {

                idrDataSplit <- lapply(
                    idrDataSplit,
                    function(x) {
                        x$length <- x$orf_aa_end - x$orf_aa_start + 1
                        return(x)
                    }
                )

                ### Analyze overlap
                nIdr <- sapply(idrDataSplit, nrow)
                if( all(nIdr) ) {
                    idrRanges <- lapply(
                        idrDataSplit,
                        function(x) {
                            if( x$idrStartGenomic[1] < x$idrEndGenomic[1]) {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$idrStartGenomic,
                                        end   = x$idrEndGenomic
                                    ),
                                    type  = x$idr_type,
                                    orfLength = x$length
                                )
                            } else {
                                GRanges(
                                    seqnames = 'artificial',
                                    IRanges(
                                        start = x$idrEndGenomic,
                                        end   = x$idrStartGenomic
                                    ),
                                    type  = x$idr_type,
                                    orfLength = x$length
                                )
                            }

                        }
                    )

                    ### Calculate overlap
                    localOverlap <- findOverlaps(
                        query   = idrRanges[[ upIso   ]],
                        subject = idrRanges[[ downIso ]]
                    )

                    if(length(localOverlap)) {

                        localOverlapDf <- as.data.frame(localOverlap)
                        localOverlapDf$upLength <- idrRanges[[ upIso   ]]$orfLength[queryHits   (localOverlap)]
                        localOverlapDf$dnLength <- idrRanges[[ downIso ]]$orfLength[subjectHits (localOverlap)]


                        localOverlapDf$lengthDiff <-
                            abs(localOverlapDf$upLength - localOverlapDf$dnLength) > AaCutoff

                        if (!is.null(AaFracCutoff)) {
                            localOverlapDf$minLength <- apply(localOverlapDf[,c('upLength','dnLength')], 1, min)
                            localOverlapDf$maxLength <- apply(localOverlapDf[,c('upLength','dnLength')], 1, max)

                            localOverlapDf$lengthDiff <-
                                localOverlapDf$minLength / localOverlapDf$maxLength < AaFracCutoff & localOverlapDf$lengthDiff

                        }


                        localOverlapDfDiff <- localOverlapDf[which(
                            localOverlapDf$lengthDiff
                        ),]

                        differentIdrLength <- nrow(localOverlapDfDiff) > 0

                        localIndex <-
                            which(isoComparison$featureCompared == 'IDR_length')
                        isoComparison$isoformsDifferent[localIndex] <-
                            differentIdrLength

                    } else {
                        differentIdrType <- FALSE
                    }

                } else {
                    differentIdrType <- FALSE
                }

                ### Repport if any difference
                if (differentIdrType & addDescription) {

                    localOverlapDfDiff$maxIsUp <- localOverlapDfDiff$maxLength == localOverlapDfDiff$upLength

                    ### Deside consequence
                    if(
                        all( c(TRUE, FALSE) %in% localOverlapDfDiff$maxIsUp )
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            #'IDR length gain and loss'
                            'Mixed IDR length differences'
                    } else if(
                        all(localOverlapDfDiff$maxIsUp)
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'IDR length gain'
                    } else if(
                        all( ! localOverlapDfDiff$maxIsUp )
                    ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'IDR length loss'
                    } else {
                        stop('Something with idr length went wrong')
                    }

                }
            }
        }

        if ('signal_peptide_identified'       %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {
                nrSignalPeptides <- sapply(peptideDataSplit, nrow)

                differentSignaPeptied <-
                    nrSignalPeptides[1] != nrSignalPeptides[2]

                # make repport
                localIndex <-
                    which(
                        isoComparison$featureCompared ==
                            'signal_peptide_identified'
                    )
                isoComparison$isoformsDifferent[localIndex] <-
                    differentSignaPeptied

                if (differentSignaPeptied & addDescription) {
                    upHasSignal <-
                        names(nrSignalPeptides)[which.max(
                            nrSignalPeptides
                        )] == upIso

                    if (upHasSignal) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Signal peptide gain'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Signal peptide loss'
                    }
                }
            }
        }

        if ('sub_cell_location'               %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) == 2) {
                if( sum(!is.na(transcriptData$sub_cell_location)) == 2 ) {

                    upIsoIndex <- which(
                        transcriptData$isoform_id == upIso
                    )
                    dnIsoIndex <- which(
                        transcriptData$isoform_id != upIso
                    )

                    locUp <- sort( unlist( strsplit(transcriptData$sub_cell_location[upIsoIndex], ',') ))
                    locDn <- sort( unlist( strsplit(transcriptData$sub_cell_location[dnIsoIndex], ',') ))

                    differentLoc <- ! identical(locUp, locDn)

                    # make repport
                    localIndex <-
                        which(
                            isoComparison$featureCompared ==
                                'sub_cell_location'
                        )
                    isoComparison$isoformsDifferent[localIndex] <-
                        differentLoc

                    if (differentLoc & addDescription) {
                        ### Analyze overlap
                        upUniqueLength <- length(setdiff(locUp, locDn))
                        #shared   <- intersect(locUp, locDn)
                        dnUniqueLength <- length(setdiff(locDn, locUp))

                        if( upUniqueLength  > 0 & dnUniqueLength > 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location switch'
                        }
                        if( upUniqueLength  > 0 & dnUniqueLength == 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location gain'
                        }
                        if( upUniqueLength == 0 & dnUniqueLength > 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location loss'
                        }

                    }
                }
            }
        }

        if ('sub_cell_shift_to_cell_membrane' %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) == 2) {
                if( sum(!is.na(transcriptData$sub_cell_location)) == 2 ) {

                    upIsoIndex <- which(
                        transcriptData$isoform_id == upIso
                    )
                    dnIsoIndex <- which(
                        transcriptData$isoform_id != upIso
                    )

                    upLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[upIsoIndex], ',') ),
                        'Cell_membrane'
                    ))
                    dnLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[dnIsoIndex], ',') ),
                        'Cell_membrane'
                    ))

                    differentLoc <- upLength != dnLength

                    # make repport
                    localIndex <-
                        which(
                            isoComparison$featureCompared ==
                                'sub_cell_shift_to_cell_membrane'
                        )
                    isoComparison$isoformsDifferent[localIndex] <-
                        differentLoc

                    if (differentLoc & addDescription) {

                        if( upLength  > 0 & dnLength == 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location memb gain'
                        }
                        if( upLength == 0 & dnLength > 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location memb loss'
                        }

                    }
                }
            }
        }

        if ('sub_cell_shift_to_cytoplasm'     %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) == 2) {
                if( sum(!is.na(transcriptData$sub_cell_location)) == 2 ) {

                    upIsoIndex <- which(
                        transcriptData$isoform_id == upIso
                    )
                    dnIsoIndex <- which(
                        transcriptData$isoform_id != upIso
                    )

                    upLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[upIsoIndex], ',') ),
                        'Cytoplasm'
                    ))
                    dnLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[dnIsoIndex], ',') ),
                        'Cytoplasm'
                    ))

                    differentLoc <- upLength != dnLength

                    # make repport
                    localIndex <-
                        which(
                            isoComparison$featureCompared ==
                                'sub_cell_shift_to_cytoplasm'
                        )
                    isoComparison$isoformsDifferent[localIndex] <-
                        differentLoc

                    if (differentLoc & addDescription) {

                        if( upLength  > 0 & dnLength == 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location cyto gain'
                        }
                        if( upLength == 0 & dnLength > 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location cyto loss'
                        }

                    }
                }
            }
        }

        if ('sub_cell_shift_to_nucleus'       %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) == 2) {
                if( sum(!is.na(transcriptData$sub_cell_location)) == 2 ) {

                    upIsoIndex <- which(
                        transcriptData$isoform_id == upIso
                    )
                    dnIsoIndex <- which(
                        transcriptData$isoform_id != upIso
                    )

                    upLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[upIsoIndex], ',') ),
                        'Nucleus'
                    ))
                    dnLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[dnIsoIndex], ',') ),
                        'Nucleus'
                    ))

                    differentLoc <- upLength != dnLength

                    # make report
                    localIndex <-
                        which(
                            isoComparison$featureCompared ==
                                'sub_cell_shift_to_nucleus'
                        )
                    isoComparison$isoformsDifferent[localIndex] <-
                        differentLoc

                    if (differentLoc & addDescription) {

                        if( upLength  > 0 & dnLength == 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location nucl gain'
                        }
                        if( upLength == 0 & dnLength > 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location nucl loss'
                        }

                    }
                }
            }
        }

        if ('sub_cell_shift_to_Extracellular' %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) == 2) {
                if( sum(!is.na(transcriptData$sub_cell_location)) == 2 ) {

                    upIsoIndex <- which(
                        transcriptData$isoform_id == upIso
                    )
                    dnIsoIndex <- which(
                        transcriptData$isoform_id != upIso
                    )

                    upLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[upIsoIndex], ',') ),
                        'Extracellular'
                    ))
                    dnLength <- length(intersect(
                        unlist( strsplit(transcriptData$sub_cell_location[dnIsoIndex], ',') ),
                        'Extracellular'
                    ))

                    differentLoc <- upLength != dnLength

                    # make report
                    localIndex <-
                        which(
                            isoComparison$featureCompared ==
                                'sub_cell_shift_to_Extracellular'
                        )
                    isoComparison$isoformsDifferent[localIndex] <-
                        differentLoc

                    if (differentLoc & addDescription) {

                        if( upLength  > 0 & dnLength == 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location ext cell gain'
                        }
                        if( upLength == 0 & dnLength > 0 ) {
                            isoComparison$switchConsequence[localIndex] <- 'SubCell location ext cell loss'
                        }

                    }
                }
            }
        }

        if ('isoform_topology'                %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {
                toplogy <- lapply(topDataSplit, function(x)
                    x$region_type)

                t1 <- rle(toplogy[[1]])
                t2 <- rle(toplogy[[2]])

                differentTopology <- !identical(t1, t2)
                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'isoform_topology')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentTopology

                if (differentTopology & addDescription) {
                    if (sum(t1$lengths) != sum(t2$lengths)) {
                        nrTop <- sapply(topDataSplit, nrow)
                        upHasMoreTop <-
                            names(nrTop)[which.max(nrTop)] == upIso

                        if (upHasMoreTop) {
                            isoComparison$switchConsequence[localIndex] <-
                                'Topology complexity gain'
                        } else {
                            isoComparison$switchConsequence[localIndex] <-
                                'Topology complexity loss'
                        }
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            NA
                    }
                }
            }
        }

        if ('extracellular_region_count'      %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {

                upCount <- sum( topDataSplit[[upIso  ]]$region_type == 'outside')
                dnCount <- sum( topDataSplit[[downIso]]$region_type == 'outside')

                differentTopology <- upCount != dnCount
                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'extracellular_region_count')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentTopology

                if (differentTopology & addDescription) {
                    if ( upCount > dnCount ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Extracellular region gain'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Extracellular region loss'
                    }
                }
            }
        }

        if ('intracellular_region_count'      %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {

                upCount <- sum( topDataSplit[[upIso  ]]$region_type == 'inside')
                dnCount <- sum( topDataSplit[[downIso]]$region_type == 'inside')

                differentTopology <- upCount != dnCount
                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'intracellular_region_count')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentTopology

                if (differentTopology & addDescription) {
                    if ( upCount > dnCount ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Intracellular region gain'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Intracellular region loss'
                    }
                }
            }
        }

        if ('extracellular_region_length'     %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {

                upLength <- sum( topDataSplit[[upIso  ]]$regionLength[which(
                    topDataSplit[[upIso  ]]$region_type == 'outside'
                )] )
                dnLength <- sum( topDataSplit[[downIso ]]$regionLength[which(
                    topDataSplit[[downIso]]$region_type == 'outside'
                )] )

                lengthGain <- upLength - dnLength
                minLength <- min(c(upLength, dnLength))
                maxLength <- max(c(upLength, dnLength))

                differentLength <- abs(lengthGain) > AaCutoff & (minLength / maxLength) < AaFracCutoff

                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'extracellular_region_length')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentLength

                if (differentLength & addDescription) {
                    if ( lengthGain > 0 ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Extracellular length gain'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Extracellular length loss'
                    }
                }
            }
        }

        if ('intracellular_region_length'     %in% consequencesToAnalyze) {
            if (sum(!is.na(transcriptData$orfTransciptLength)) > 0) {

                upLength <- sum( topDataSplit[[upIso  ]]$regionLength[which(
                    topDataSplit[[upIso  ]]$region_type == 'inside'
                )] )
                dnLength <- sum( topDataSplit[[downIso ]]$regionLength[which(
                    topDataSplit[[downIso]]$region_type == 'inside'
                )] )

                lengthGain <- upLength - dnLength
                minLength <- min(c(upLength, dnLength))
                maxLength <- max(c(upLength, dnLength))

                differentLength <- abs(lengthGain) > AaCutoff & (minLength / maxLength) < AaFracCutoff

                # make repport
                localIndex <-
                    which(isoComparison$featureCompared == 'intracellular_region_length')
                isoComparison$isoformsDifferent[localIndex] <-
                    differentLength

                if (differentLength & addDescription) {
                    if ( lengthGain > 0 ) {
                        isoComparison$switchConsequence[localIndex] <-
                            'Intracellular length gain'
                    } else {
                        isoComparison$switchConsequence[localIndex] <-
                            'Intracellular length loss'
                    }
                }
            }
        }


    }

    if (onlyRepportDifferent) {
        isoComparison <-
            isoComparison[which(
                isoComparison$isoformsDifferent
        ),]
    }

    return(isoComparison)
}


### Summarizing consequences
extractConsequenceSummary <- function(
    switchAnalyzeRlist,
    consequencesToAnalyze = 'all',
    includeCombined = FALSE,
    asFractionTotal = FALSE,
    alpha = 0.05,
    dIFcutoff = 0.1,
    plot = TRUE,
    plotGenes = FALSE,
    simplifyLocation = TRUE,
    removeEmptyConsequences = FALSE,
    localTheme = theme_bw(),
    returnResult = FALSE
) {
    ### Test input
    if (TRUE) {
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(
                'The object supplied to \'switchAnalyzeRlist\' must be a \'switchAnalyzeRlist\''
            )
        }

        # test wether switching have been analyzed
        if (!any(!is.na(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value
        ))) {
            stop(
                'The analsis of isoform switching must be performed before functional consequences can be analyzed. Please run detectIsoformSwitching() and try again.'
            )
        }
        # test whether switches have been predicted
        if (is.null(switchAnalyzeRlist$switchConsequence)) {
            stop(
                'The analsis of isoform switch consequences must be performed before it can be summarized. Please use analyzeSwitchConsequences() and try again.'
            )
        }

        # input format
        if (dIFcutoff < 0 | dIFcutoff > 1) {
            stop('The dIFcutoff must be in the interval [0,1].')
        }
        if (alpha < 0 |
            alpha > 1) {
            warning('The alpha parameter should usually be between 0 and 1 ([0,1]).')
        }
        if (alpha > 0.05) {
            warning(
                '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'
            )
        }

        ### Consequences to analyze
        acceptedTypes <- c(
            # Transcript
            'tss',
            'tts',
            'last_exon',
            'isoform_length',
            'exon_number',
            'intron_structure',
            'intron_retention',
            'isoform_class_code',
            # cpat
            'coding_potential',
            # ORF
            'ORF_genomic',
            'ORF_length',
            '5_utr_length',
            '3_utr_length',
            # seq similarity
            'isoform_seq_similarity',
            'ORF_seq_similarity',
            '5_utr_seq_similarity',
            '3_utr_seq_similarity',
            # ORF
            'NMD_status',
            # pfam
            'domains_identified',
            'genomic_domain_position',
            'domain_length',
            'domain_isotype',

            # SignalIP
            'signal_peptide_identified',

            # IDR
            'IDR_identified',
            'IDR_length',
            'IDR_type',

            # sub cell
            'sub_cell_location',
            'sub_cell_shift_to_cell_membrane',
            'sub_cell_shift_to_cytoplasm',
            'sub_cell_shift_to_nucleus',
            'sub_cell_shift_to_Extracellular',

            # topology
            'isoform_topology',
            'extracellular_region_count',
            'intracellular_region_count',
            'extracellular_region_length',
            'intracellular_region_length'
        )

        if (!all(consequencesToAnalyze %in% c('all', acceptedTypes))) {
            stop(
                'The argument(s) supplied to \'typeOfconsequence\' are not accepted. Please see ?analyzeSwitchConsequences under details for description of which strings are allowed.'
            )
        }

        consequencesAnalyzed <-
            unique(switchAnalyzeRlist$switchConsequence$featureCompared)
        if ('all' %in% consequencesToAnalyze) {
            consequencesToAnalyze <- consequencesAnalyzed
        }

        consequencesNotAnalyzed <-
            setdiff(consequencesToAnalyze, consequencesAnalyzed)
        if (length(consequencesNotAnalyzed)) {
            warning(
                paste(
                    'The following consequences appear not to have been analyzed and will therefor not be summarized:',
                    paste(consequencesNotAnalyzed, collapse = ', '),
                    sep = ' '
                )
            )
        }


    }

    ### Massage for plotting
    localSwitchConsequences <-
        switchAnalyzeRlist$switchConsequence[which(
            switchAnalyzeRlist$switchConsequence$isoformsDifferent &
                switchAnalyzeRlist$switchConsequence$featureCompared %in%
                consequencesToAnalyze
        ),]
    if (!nrow(localSwitchConsequences)) {
        stop('No swithces with consequences were found')
    }

    ### Subset to those with switches
    localSwitchConsequences <-
        localSwitchConsequences[which(!is.na(
            localSwitchConsequences$switchConsequence
        )), ]
    if (!nrow(localSwitchConsequences)) {
        stop('No swithces with describable consequences were found')
    }
    localSwitchConsequences$Comparison <-
        paste(
            localSwitchConsequences$condition_1,
            'vs',
            localSwitchConsequences$condition_2,
            sep = ' '
        )

    if (includeCombined) {
        tmp <- localSwitchConsequences

        tmp$featureCompared <- 'combined'
        tmp$switchConsequence <- 'any consequence'
        tmp <- unique(tmp)

        localSwitchConsequences <-
            rbind(localSwitchConsequences, tmp)
    }

    ### Extract Sig iso
    isoResTest <-
        any(!is.na(
            switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
        ))
    if (isoResTest) {
        sigIso <- switchAnalyzeRlist$isoformFeatures[which(
            switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value < alpha &
                abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
        ),
        c('iso_ref', 'gene_ref')]
    } else {
        sigIso <- switchAnalyzeRlist$isoformFeatures[which(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value < alpha &
                abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
        ),
        c('iso_ref', 'gene_ref')]
    }

    ### summarize count
    myNumbers <-
        plyr::ddply(
            localSwitchConsequences,
            .variables = c(
                'switchConsequence',
                'featureCompared',
                'Comparison'
            ),
            .drop = TRUE,
            .fun = function(aDF) {
                data.frame(
                    featureCompared = aDF$featureCompared[1],
                    nrGenesWithConsequences = length(
                        intersect(aDF$gene_ref, sigIso$gene_ref)
                    ),
                    nrIsoWithConsequences = length(intersect(
                        c(aDF$iso_ref_up, aDF$iso_ref_down), sigIso$iso_ref
                    )),
                    stringsAsFactors = FALSE
                )
            }
        )

    ### Add those analyzed but with zero consequences
    if( ! removeEmptyConsequences ) {
        notInData <- setdiff(
            consequencesAnalyzed,
            unique(myNumbers$featureCompared)
        )
        if(length(notInData)) {

            missingData <- unique(
                switchAnalyzeRlist$switchConsequence[
                    which(switchAnalyzeRlist$switchConsequence$featureCompared %in% notInData),
                    c('condition_1','condition_2','featureCompared','switchConsequence')
                ]
            )

            missingData$Comparison <-
                paste(
                    missingData$condition_1,
                    'vs',
                    missingData$condition_2,
                    sep = ' '
                )

            missingData$switchConsequence <- 'Any consequence'

            missingData <- missingData[,c(
                'switchConsequence',
                'Comparison',
                'featureCompared'
            )]

            missingData$nrGenesWithConsequences <- 0
            missingData$nrIsoWithConsequences <- 0

            myNumbers <- rbind(
                myNumbers,
                missingData
            )
        }
    }


    myNumbers$featureCompared <-
        startCapitalLetter(gsub('_', ' ', myNumbers$featureCompared))
    myNumbers <-
        myNumbers[, c(
            'Comparison',
            'featureCompared',
            'switchConsequence',
            'nrGenesWithConsequences',
            'nrIsoWithConsequences'
        )]

    if (asFractionTotal) {
        totalNrSwitches <-
            extractSwitchSummary(switchAnalyzeRlist,
                                 alpha = alpha,
                                 dIFcutoff = dIFcutoff)

        ### add to numbers
        matchVec <-
            match(myNumbers$Comparison, totalNrSwitches$Comparison)
        myNumbers$nrIsoforms <- totalNrSwitches$nrIsoforms[matchVec]
        myNumbers$nrGenes <- totalNrSwitches$nrGenes[matchVec]

        myNumbers$geneFraction <-
            round(myNumbers$nrGenesWithConsequences / myNumbers$nrGenes,
                  digits = 4)
        myNumbers$isoFraction  <-
            round(myNumbers$nrIsoWithConsequences   / myNumbers$nrIsoforms,
                  digits = 4)
    }

    if (plot) {
        myNumbers$plotComparison <-
            gsub(' vs ', '\nvs\n', myNumbers$Comparison)

        myNumbers$featureCompared <- newLineAtMiddel(
            myNumbers$featureCompared
        )

        ### Make plot
        if (plotGenes) {
            if (asFractionTotal) {
                g1 <- ggplot(myNumbers, aes(
                    x = switchConsequence,
                    y = geneFraction
                )) +
                    labs(
                        x = 'Consequence of switch\n(features of the upregulated isoform)',
                        y = 'Fraction of genes'
                    )
            } else {
                g1 <-
                    ggplot(myNumbers, aes(
                        x = switchConsequence,
                        y = nrGenesWithConsequences
                    )) +
                    labs(
                        x = 'Consequence of switch\n(features of the upregulated isoform)',
                        y = 'Number of genes'
                        )
            }
        } else {
            if (asFractionTotal) {
                g1 <- ggplot(
                    myNumbers,
                    aes(x = switchConsequence, y = isoFraction)
                ) + labs(
                    x = 'Consequence of switch\n(features of the upregulated isoform)',
                    y = 'Fraction of isoforms')
            } else {
                g1 <- ggplot(
                    myNumbers,
                    aes(x = switchConsequence, y = nrIsoWithConsequences)
                ) + labs(
                    x = 'Consequence of switch\n(features of the upregulated isoform)',
                    y = 'Number of isoforms'
                )
            }
        }

        g1 <- g1 + geom_bar(stat = "identity") +
            facet_grid(plotComparison ~ featureCompared,
                       scales = 'free',
                       space = 'free_x') +
            localTheme +
            theme(strip.text.y = element_text(angle = 0)) +
            theme(axis.text.x = element_text(
                angle = -45,
                hjust = 0,
                vjust = 1
            ))

        myNumbers$plotComparison <- NULL
    }

    if (returnResult) {
        if(plot) {
            print(g1)
        }

        return(myNumbers)
    } else {
        if(plot) {
            return(g1)
        }
    }

}

extractConsequenceEnrichment <- function(
    switchAnalyzeRlist,
    consequencesToAnalyze = 'all',
    alpha=0.05,
    dIFcutoff = 0.1,
    countGenes = TRUE,
    analysisOppositeConsequence=FALSE,
    plot=TRUE,
    localTheme = theme_bw(base_size = 12),
    minEventsForPlotting = 10,
    returnResult=TRUE,
    returnSummary=TRUE
) {
    ### Test input
    if (TRUE) {
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(
                'The object supplied to \'switchAnalyzeRlist\' must be a \'switchAnalyzeRlist\''
            )
        }

        # test wether switching have been analyzed
        if (!any(!is.na(
            switchAnalyzeRlist$isoformFeatures$gene_switch_q_value
        ))) {
            stop(
                'The analsis of isoform switching must be performed before functional consequences can be analyzed. Please run detectIsoformSwitching() and try again.'
            )
        }
        # test whether switches have been predicted
        if (is.null(switchAnalyzeRlist$switchConsequence)) {
            stop(
                'The analsis of isoform switch consequences must be performed before it can be summarized. Please use analyzeSwitchConsequences() and try again.'
            )
        }

        # input format
        if (dIFcutoff < 0 | dIFcutoff > 1) {
            stop('The dIFcutoff must be in the interval [0,1].')
        }
        if (alpha < 0 |
            alpha > 1) {
            warning('The alpha parameter should usually be between 0 and 1 ([0,1]).')
        }
        if (alpha > 0.05) {
            warning(
                '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'
            )
        }

        ### Consequences to analyze
        acceptedTypes <- c(
            # Transcript
            'tss',
            'tts',
            'last_exon',
            'isoform_length',
            'exon_number',
            'intron_structure',
            'intron_retention',
            'isoform_class_code',
            # cpat
            'coding_potential',
            # ORF
            'ORF_genomic',
            'ORF_length',
            '5_utr_length',
            '3_utr_length',
            # seq similarity
            'isoform_seq_similarity',
            'ORF_seq_similarity',
            '5_utr_seq_similarity',
            '3_utr_seq_similarity',
            # ORF
            'NMD_status',
            # pfam
            'domains_identified',
            'genomic_domain_position',
            'domain_length',
            'domain_isotype',

            # SignalIP
            'signal_peptide_identified',

            # IDR
            'IDR_identified',
            'IDR_length',
            'IDR_type',

            # sub cell
            'sub_cell_location',
            'sub_cell_shift_to_cell_membrane',
            'sub_cell_shift_to_cytoplasm',
            'sub_cell_shift_to_nucleus',
            'sub_cell_shift_to_Extracellular',

            # topology
            'isoform_topology'
        )

        if (!all(consequencesToAnalyze %in% c('all', acceptedTypes))) {
            stop(
                'The argument(s) supplied to \'typeOfconsequence\' are not accepted. Please see ?summarizeSwitchConsequences under details for description of which strings are allowed.'
            )
        }

        consequencesAnalyzed <-
            unique(switchAnalyzeRlist$switchConsequence$featureCompared)
        if ('all' %in% consequencesToAnalyze) {
            consequencesToAnalyze <- consequencesAnalyzed
        }

        consequencesNotAnalyzed <-
            setdiff(consequencesToAnalyze, consequencesAnalyzed)
        if (length(consequencesNotAnalyzed)) {
            warning(
                paste(
                    'The following consequences appear not to have been analyzed and will therefor not be summarized:',
                    paste(consequencesNotAnalyzed, collapse = ', '),
                    sep = ' '
                )
            )
        }

    }

    ### Extract non-location consequences
    if(TRUE) {
        localConseq <- switchAnalyzeRlist$switchConsequence[
            which( !is.na(
                switchAnalyzeRlist$switchConsequence$switchConsequence
            ))
            ,]
        localConseq <- localConseq[which(
            ! grepl('switch', localConseq$switchConsequence)
        ),]

        ### make list with levels
        levelList <- list(
            tss=c('Tss more upstream','Tss more downstream'),
            tts=c('Tts more downstream','Tts more upstream'),
            last_exon=c('Last exon more downstream','Last exon more upstream'),
            isoform_length=c('Length gain','Length loss'),
            isoform_seq_similarity=c('Length gain','Length loss'),
            exon_number=c('Exon gain','Exon loss'),
            intron_retention=c('Intron retention gain','Intron retention loss'),
            ORF_length=c('ORF is longer','ORF is shorter'),
            ORF=c('Complete ORF loss','Complete ORF gain'),
            x5_utr_length=c('5UTR is longer','5UTR is shorter'),
            x3_utr_length=c('3UTR is longer','3UTR is shorter'),
            NMD_status=c('NMD sensitive','NMD insensitive'),
            coding_potential=c('Transcript is coding','Transcript is Noncoding'),
            domains_identified=c('Domain gain','Domain loss'),
            domain_length=c('Domain length gain','Domain length loss'),
            domain_isotype=c('Domain non-reference isotype gain','Domain non-reference isotype loss'),
            IDR_identified = c('IDR gain','IDR loss'),
            IDR_length = c('IDR length gain', 'IDR length loss'),
            IDR_type = c('IDR w binding region gain', 'IDR w binding region loss'),
            signal_peptide_identified=c('Signal peptide gain','Signal peptide loss'),
            sub_cell_location = c('SubCell location gain','SubCell location loss'),
            sub_cell_shift_to_cell_membrane = c('SubCell location memb gain','SubCell location memb loss'),
            sub_cell_shift_to_cytoplasm     = c('SubCell location cyto gain','SubCell location cyto loss'),
            sub_cell_shift_to_nucleus       = c('SubCell location nucl gain','SubCell location nucl loss'),
            sub_cell_shift_to_Extracellular = c('SubCell location ext cell gain','SubCell location ext cell loss'),
            isoform_topology = c('Topology complexity gain','Topology complexity loss'),
            extracellular_region_count = c('Extracellular region gain', 'Extracellular region loss'),
            intracellular_region_count = c('Intracellular region gain', 'Intracellular region loss'),
            extracellular_region_length = c('Extracellular length gain','Extracellular length loss'),
            intracellular_region_length = c('Intracellular length gain','Intracellular length loss')
        )
        levelListDf <- plyr::ldply(levelList, function(x) data.frame(feature=x, stringsAsFactors = FALSE))

        ### Add consequence paris
        localConseq$conseqPair <- levelListDf$.id[match(localConseq$switchConsequence, levelListDf$feature)]
        localConseq <- localConseq[which( !is.na(localConseq$conseqPair)),]

        ### Subset to consequences analyzed
        localConseq <- localConseq[which(
            localConseq$featureCompared %in% consequencesToAnalyze
        ),]
    }

    ### Subset to significant features
    if(TRUE) {
        ### Extract Sig iso
        isoResTest <-
            any(!is.na(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
            ))
        if (isoResTest) {
            sigIso <- switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value < alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            c('iso_ref', 'gene_ref')]
        } else {
            sigIso <- switchAnalyzeRlist$isoformFeatures[which(
                switchAnalyzeRlist$isoformFeatures$gene_switch_q_value < alpha &
                    abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
            ),
            c('iso_ref', 'gene_ref')]
        }

        if(isoResTest) {
            localConseq <- localConseq[which(
                localConseq$iso_ref_down %in% sigIso$iso_ref |
                localConseq$iso_ref_up   %in% sigIso$iso_ref
            ),]
        } else {
            localConseq <- localConseq[which(
                localConseq$gene_ref %in% sigIso$gene_ref
            ),]
        }
    }


    ### Summarize gain vs loss for each consequence in each condition
    consequenceBalance <- plyr::ddply(
        .data = localConseq,
        .variables = c('condition_1','condition_2','conseqPair'),
        #.inform = TRUE,
        .fun = function(aDF) { # aDF <- localConseq[1:20,]
            ### Add levels
            if(analysisOppositeConsequence) {
                localLvl <- rev(sort(
                    levelList[[ aDF$conseqPair[1] ]]
                ))
            } else {
                localLvl <- sort(
                    levelList[[ aDF$conseqPair[1] ]]
                )
            }
            aDF$switchConsequence <- factor(
                aDF$switchConsequence,
                levels=localLvl
            )

            ### Summarize category
            if( countGenes ) {
                df2 <- aDF[
                    which(!is.na(aDF$switchConsequence)),
                    c('gene_id','switchConsequence')
                ]
                localNumber <- plyr::ddply(df2, .drop = FALSE, .variables = 'switchConsequence', function(x) {
                    data.frame(
                        Freq = length(unique(x$gene_id))
                    )
                })
                colnames(localNumber)[1] <- 'Var1'
            } else {
                localNumber <- as.data.frame(table(aDF$switchConsequence))
            }

            if(nrow(localNumber) == 2) {
                localTest <- suppressWarnings(
                    stats::binom.test(localNumber$Freq[1], sum(localNumber$Freq))
                )

                localRes <- data.frame(
                    feature=stringr::str_c(
                        localNumber$Var1[1],
                        ' (paired with ',
                        localNumber$Var1[2],
                        ')'
                    ),
                    propOfRelevantEvents=localTest$estimate,
                    stringsAsFactors = FALSE
                )

                localRes$propCiLo <- min(localTest$conf.int)
                localRes$propCiHi <- max(localTest$conf.int)
                localRes$propPval <- localTest$p.value
            } else {
                warning('Somthing strange happend - contact developer with reproducible example')
            }

            localRes$nUp   <- localNumber$Freq[which( localNumber$Var1 == levels(localNumber$Var1)[1] )]
            localRes$nDown <- localNumber$Freq[which( localNumber$Var1 == levels(localNumber$Var1)[2] )]

            return(localRes)
        }
    )

    consequenceBalance$propQval <- p.adjust(consequenceBalance$propPval, method = 'fdr')
    consequenceBalance$Significant <- consequenceBalance$propQval < alpha
    consequenceBalance$Significant <- factor(
        consequenceBalance$Significant,
        levels=c(FALSE,TRUE)
    )

    ### Plot result
    if(plot) {
        if(countGenes) {
            xText <- 'Fraction of Genes Having the Consequence Indicated\n(of Genes Affected by Either of Opposing Consequences)\n(With 95% Confidence Interval)'
        } else {
            xText <- 'Fraction of Switches Having the Consequence Indicated\n(of Switches Affected by Either of Opposing Consequences)\n(With 95% Confidence Interval)'
        }

        consequenceBalance2 <- consequenceBalance[which(
            (consequenceBalance$nUp + consequenceBalance$nDown) >= minEventsForPlotting
        ),]
        if(nrow(consequenceBalance2) == 0) {
            stop('No features left for ploting after filtering with via "minEventsForPlotting" argument.')
        }

        consequenceBalance2$nTot <- consequenceBalance2$nDown + consequenceBalance2$nUp

        ### Add comparison
        consequenceBalance2$Comparison <- paste(
            consequenceBalance2$condition_1,
            'vs',
            consequenceBalance2$condition_2,
            sep='\n'
        )

        ### Massage
        consequenceBalance2$feature2 <- gsub(' \\(', '\n(', consequenceBalance2$feature)

        consequenceBalance2$feature2 <- factor(
            consequenceBalance2$feature2,
            levels = rev(sort(unique(as.character(consequenceBalance2$feature2))))
        )

        g1 <- ggplot(data=consequenceBalance2, aes(y=feature2, x=propOfRelevantEvents, color=Significant)) +
            #geom_point(size=4) +
            geom_errorbarh(aes(xmax = propCiLo, xmin=propCiHi), height = .3) +
            geom_point(aes(size=nTot)) +
            facet_wrap(~Comparison) +
            geom_vline(xintercept=0.5, linetype='dashed') +
            labs(
                x=xText,
                y='Consequence of Isoform Switch\n(and the opposing consequence)'
            ) +
            localTheme +
            theme(axis.text.x=element_text(angle=-45, hjust = 0, vjust=1)) +
            scale_color_manual(name = paste0('FDR < ', alpha), values=c('black','red'), drop=FALSE) +
            scale_radius(limits=c(
                0,
                # Ensure largest number is always incuded
                max(
                    roundUpToNearsTenOrHundred( consequenceBalance2$nTot )
                )
            )) +
            guides(
                color = guide_legend(order=1),
                size = guide_legend(order=2)
            ) +
            coord_cartesian(xlim=c(0,1))

        if( countGenes ) {
            g1 <- g1 + labs(size = 'Genes')
        } else {
            g1 <- g1 + labs(size = 'Switches')
        }
    }

    if(returnResult) {
        if(plot) {
            print(g1)
        }

        if(returnSummary) {
            consequenceBalance$Comparison <- NULL

            return(consequenceBalance)

        } else {
            return(localConseq)
        }


    } else {
        if(plot) {
            return(g1)
        }
    }
}

extractConsequenceEnrichmentComparison <- function(
    switchAnalyzeRlist,
    consequencesToAnalyze = 'all',
    alpha=0.05,
    dIFcutoff = 0.1,
    countGenes = TRUE,
    analysisOppositeConsequence=FALSE,
    plot=TRUE,
    localTheme = theme_bw(base_size = 14),
    minEventsForPlotting = 10,
    returnResult=TRUE
) {
    ### Test
    if(TRUE) {
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(
                'The object supplied to \'switchAnalyzeRlist\' must be a \'switchAnalyzeRlist\''
            )
        }
        nComp <- nrow( unique(
            switchAnalyzeRlist$isoformFeatures[,c('condition_1','condition_2')]
        ))

        if( nComp == 1) {
            stop('Cannot do a contrast of different comparisons since only one comparison is analyzed in the switchAnalyzeRlist.')
        }
    }


    ### Extract splicing enrichment
    conseqCount <- extractConsequenceEnrichment(
        switchAnalyzeRlist = switchAnalyzeRlist,
        consequencesToAnalyze = consequencesToAnalyze,
        alpha = alpha,
        dIFcutoff = dIFcutoff,
        countGenes = countGenes,
        minEventsForPlotting = 0,
        analysisOppositeConsequence=analysisOppositeConsequence,
        plot = FALSE,
        returnResult = TRUE
    )
    conseqCount$Comparison <- stringr::str_c(
        conseqCount$condition_1,
        '\nvs\n',
        conseqCount$condition_2
    )
    conseqCount$nTot <- conseqCount$nDown + conseqCount$nUp

    ### Extract pairs
    conseqPairs <- split(as.character(unique(conseqCount$feature)), unique(conseqCount$feature))

    ### Make each pairwise comparison
    myComparisons <- allPairwiseFeatures( unique(conseqCount$Comparison) )

    ### Loop over each pariwise comparison
    fisherRes <- plyr::ddply(myComparisons, c('var1','var2'), function(localComparison) { # localComparison <- myComparisons[1,]
        ### Loop over each
        localAsRes <- plyr::ldply(conseqPairs, .inform = TRUE, function(localConseq) { # localConseq <- 'Membrane tethering gain (paired with Membrane tethering loss)'
            ### Extract local data
            localSpliceCount1 <- conseqCount[which(
                conseqCount$Comparison %in% c(localComparison$var1, localComparison$var2) &
                    conseqCount$feature == localConseq
            ),]

            if(
                nrow(localSpliceCount1) != 2 |
                any(is.na(localSpliceCount1[,c('nUp','nDown')]))
            ) {
                return(NULL)
            }

            ### Test difference
            fisherResult <- fisher.test(localSpliceCount1[,c('nUp','nDown')])

            ###
            fisherTestResult <- data.frame(odds_ratio=fisherResult$estimate, p_value=fisherResult$p.value, lowCI=fisherResult$conf.int[1], highCI=fisherResult$conf.int[2])
            rownames(fisherTestResult) <- NULL

            localSpliceCount1$fisherOddsRatio <- fisherTestResult$odds_ratio
            localSpliceCount1$fisherPvalue    <- fisherTestResult$p_value

            localSpliceCount1$pair <- 1:2
            localSpliceCount1$forPlotting <- any(
                (localSpliceCount1$nUp + localSpliceCount1$nDown) >= minEventsForPlotting
            )

            return(
                localSpliceCount1[,c('Comparison','propOfRelevantEvents','propCiLo','propCiHi','fisherOddsRatio','fisherPvalue','pair','forPlotting','nTot')]
            )
        })

        colnames(localAsRes)[1] <- 'consequence'

        return(localAsRes)
    })

    ### Add comparison
    fisherRes$comp <- paste(
        gsub('\\n',' ',fisherRes$var1),
        gsub('\\n',' ',fisherRes$var2),
        sep='\ncompared to\n'
    )

    ### Multiple test correction
    # perform correction for pair = 1 (to avoid correcting twice for the same pair)
    tmp <- fisherRes[which(fisherRes$pair == 1),]
    tmp$fisherQvalue <- p.adjust(tmp$fisherPvalue, method = 'fdr')

    fisherRes$fisherQvalue <- tmp$fisherQvalue[match(
        stringr::str_c(fisherRes$var1, fisherRes$var2, fisherRes$consequence),
        stringr::str_c(tmp$var1, tmp$var2, tmp$consequence)
    )]
    fisherRes$Significant <- fisherRes$fisherQvalue < alpha
    fisherRes$Significant <- factor(
        fisherRes$Significant,
        levels=c(FALSE,TRUE)
    )

    ### Plot
    if(plot) {
        fisherRes2 <- fisherRes[which(fisherRes$forPlotting),]
        if(nrow(fisherRes2) == 0) {
            stop('No features left to plot after subsetting with \'minEventsForPlotting\'.')
        }

        fisherRes2$consequence <- gsub(' \\(paired','\n\\(paired', fisherRes2$consequence)
        fisherRes2$consequence <- gsub('with ','with\n', fisherRes2$consequence)

        if(countGenes) {
            xText <- 'Fraction of genes having the consequence indicated\n(of the genes affected by either of opposing consequences)\n(with 95% confidence interval)'
        } else {
            xText <- 'Fraction of switches having the consequence indicated\n(of the switches affected by either of opposing consequences)\n(with 95% confidence interval)'
        }

        g1 <- ggplot(data=fisherRes2, aes(y=Comparison, x=propOfRelevantEvents, color=Significant)) +
            #geom_point(aes(size=4)) +
            geom_point(aes(size=nTot)) +
            geom_errorbarh(aes(xmax = propCiHi, xmin=propCiLo), height = .3) +
            facet_grid(comp~consequence, scales = 'free_y') +
            geom_vline(xintercept=0.5, linetype='dashed') +
            labs(x=xText, y='Comparison') +
            #scale_color_manual('Fraction in\nComparisons\nSignifcantly different', values=c('black','red'), drop=FALSE) +
            scale_color_manual(name = paste0('FDR across\ncomparisons\n< ', alpha), values=c('black','red'), drop=FALSE) +
            guides(
                color = guide_legend(order=1),
                size = guide_legend(order=2)
            ) +
            localTheme +
            theme(axis.text.x=element_text(angle=-45, hjust = 0, vjust=1), strip.text.y = element_text(angle = 0)) +
            coord_cartesian(xlim=c(0,1))

        if( countGenes ) {
            g1 <- g1 + scale_size_continuous(name = 'Genes')
        } else {
            g1 <- g1 + scale_size_continuous(name = 'Switches')
        }
    }

    if(returnResult) {
        if(plot) {
            print(g1)
        }

        fisherRes$nTot <- NULL

        fisherRes$pair <- stringr::str_c('propUp_comparison_', fisherRes$pair)
        colnames(fisherRes)[which(colnames(fisherRes) == 'propOfRelevantEvents')] <- 'propUp'

        fisherRes2 <- reshape2::dcast(data = fisherRes, comp + consequence ~ pair, value.var=c('propUp'))

        matchIndex <- match(
            stringr::str_c(fisherRes2$comp, fisherRes2$consequence),
            stringr::str_c(tmp$comp, tmp$consequence)
        )
        #fisherRes2$fisherOddsRatio <- tmp$fisherOddsRatio[matchIndex]
        fisherRes2$fisherQvalue    <- tmp$fisherQvalue   [matchIndex]

        fisherRes2$Significant <- fisherRes2$fisherQvalue < alpha
        colnames(fisherRes2)[1] <- 'comparisonsCompared'
        fisherRes2$comparisonsCompared <- gsub('\\n', ' ', fisherRes2$comparisonsCompared)
        return(fisherRes2)
    } else {
        if(plot) {
            return(g1)
        }
    }

}

extractConsequenceGenomeWide <- function(
    switchAnalyzeRlist,
    featureToExtract = 'isoformUsage',
    annotationToAnalyze = 'all',
    alpha=0.05,
    dIFcutoff = 0.1,
    log2FCcutoff = 1,
    violinPlot=TRUE,
    alphas=c(0.05, 0.001),
    localTheme=theme_bw(),
    plot=TRUE,
    returnResult=TRUE
) {
    ### Test input
    if (TRUE) {
        if (class(switchAnalyzeRlist) != 'switchAnalyzeRlist') {
            stop(
                '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 plot conditional expression -',
                    'if that is your intention you need to create a new',
                    'switchAnalyzeRlist with the importRdata() function and start over.',
                    sep = ' '
                )
            )
        }

        if (alpha < 0 |
            alpha > 1) {
            warning('The alpha parameter should usually be between 0 and 1 ([0,1]).')
        }
        if (alpha > 0.05) {
            warning(
                '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 (log2FCcutoff < 0) {
            stop(
                'The log2FCcutoff cannot be negative (as the cutoff is applied to absolute values)'
            )
        }
        if (length(alphas) != 2) {
            stop('A vector of length 2 must be given to the argument \'alphas\'')
        }
        if (any(alphas < 0) |
            any(alphas > 1)) {
            stop('The \'alphas\' parameter must be numeric between 0 and 1 ([0,1]).')
        }
        if (any(alphas > 0.05)) {
            warning(
                '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 (!featureToExtract %in%
            c('isoformUsage', 'isoformExp', 'geneExp', 'all')
        ) {
            stop(
                'The \'featureToExtract\' argument must be  \'isoformUsage\', \'isoformExp\',  \'geneExp\' or \'all\''
            )
        }

        ### Check annotation
        okAnnot <-
            c(
                'ORF',
                'NMD_status',
                'coding_potential',
                'signal_peptide_identified',
                'domains_identified',
                'intron_retention',
                'switch_consequences',
                'isoform_class_code',
                'domain_isotype'
            )
        if ('all' %in% annotationToAnalyze) {
            annotationToAnalyze <- okAnnot
        }
        if (!all(annotationToAnalyze %in% okAnnot)) {
            stop(
                paste(
                    'The \'annotationToAnalyze\' argument must be a one (or multiple) of: \'',
                    paste(okAnnot, collapse = '\', \''),
                    '\'',
                    sep = ''
                )
            )
        }

        ### Replace with annotaion names with collum names
        annotToExtract <-
            unique(gsub('NMD_status|ORF', 'PTC', annotationToAnalyze))
        annotToExtract <- gsub(
            'coding_potential',
            'codingPotential',
            annotToExtract
        )
        annotToExtract <- gsub(
            'domains_identified'  ,
            'domain_identified',
            annotToExtract
        )
        annotToExtract <-gsub(
            'intron_retention',
            'IR',
            annotToExtract
        )
        annotToExtract <- gsub(
            'switch_consequences' ,
            'switchConsequencesGene',
            annotToExtract
        )
        annotToExtract <- gsub(
            'isoform_class_code',
            'class_code',
            annotToExtract
        )

        if (length(annotToExtract) == 0) {
            stop(
                'Somthing in the annoation decoding went wrong - please send a small example dataset reconstructing the mistake to the developers.'
            )
        }
    }

    ### Extract annotation
    if (TRUE) {
        switchAnalyzeRlist$isoformFeatures$comparison <- paste(
            switchAnalyzeRlist$isoformFeatures$condition_1,
            switchAnalyzeRlist$isoformFeatures$condition_2,
            sep = ' vs '
        )

        isoformsToAnalyze <- extractSigData(
            switchAnalyzeRlist = switchAnalyzeRlist,
            alpha = alpha,
            dIFcutoff = dIFcutoff,
            log2FCcutoff = log2FCcutoff,
            featureToExtract = featureToExtract
        )

        colToExtract <-
            c('iso_ref',
              'isoform_id',
              'comparison',
              'IF1',
              'IF2',
              annotToExtract)
        colToExtract <-
            intersect(colToExtract,
                      colnames(switchAnalyzeRlist$isoformFeatures))

        isoData <- switchAnalyzeRlist$isoformFeatures[
            which(
                switchAnalyzeRlist$isoformFeatures$iso_ref %in%
                    isoformsToAnalyze
            ),
            na.omit(match(
                colToExtract ,
                colnames(switchAnalyzeRlist$isoformFeatures)
            ))
            ]

        if('domain_isotype' %in% annotationToAnalyze) {
            structureVariantIso <- unique(
                switchAnalyzeRlist$domainAnalysis$isoform_id[which(
                    switchAnalyzeRlist$domainAnalysis$domain_isotype_simple == "Non-reference"
                )]
            )

            isoData$domain_isotype_identified <- isoData$isoform_id %in% structureVariantIso
        }


        if (nrow(isoData) == 0) {
            stop('No data left after filtering')
        }
    }

    ### Overwrite annotation with propper categories
    if (TRUE) {
        # ORF
        if (!is.null(isoData$PTC) &
            'ORF' %in% annotationToAnalyze) {
            # ORF
            isoData$ORF <- 'With ORF'
            isoData$ORF[which(is.na(isoData$PTC))] <- 'Without ORF'
        }
        # PTC
        if (!is.null(isoData$PTC) &
            'NMD_status' %in% annotationToAnalyze) {
            # PTC
            isoData$PTC <-
                ifelse(test = isoData$PTC,
                       yes = 'NMD sensitive',
                       no = 'NMD insensitive')
        } else {
            isoData$PTC <- NULL
        }
        # coding potential
        if (!is.null(isoData$codingPotential)) {
            isoData$codingPotential <-
                ifelse(test = isoData$codingPotential,
                       yes = 'Isoform is coding',
                       no = 'Isoform is non-coding')
        }
        # signal peptide
        if (!is.null(isoData$signal_peptide_identified)) {
            isoData$signal_peptide_identified <-
                ifelse(
                    test = isoData$signal_peptide_identified == 'yes',
                    yes = 'With signal peptide',
                    no = 'Without signal peptide'
                )
        }
        # protein domains
        if (!is.null(isoData$domain_identified)) {
            isoData$domain_identified <-
                ifelse(
                    test = isoData$domain_identified == 'yes',
                    yes = 'With domain',
                    no = 'Without domain'
                )
        }
        # domain structure
        if (!is.null(isoData$domain_isotype_identified)) {
            isoData$domain_isotype_identified <-
                ifelse(
                    test = isoData$domain_isotype_identified,
                    yes = 'With non-reference domain isotype',
                    no  = 'Without non-reference domain isotype'
                )
        }

        # intron retention
        if (!is.null(isoData$IR)) {
            isoData$IR <-
                ifelse(test = isoData$IR > 0,
                       yes = 'With intron retention',
                       no = 'Without intron retention')
        }
        # switch consequences
        if (!is.null(isoData$switchConsequences)) {
            isoData$switchConsequences <-
                ifelse(
                    test = isoData$switchConsequences,
                    yes = 'With switch consequence',
                    no = 'Without switch consequence'
                )
        }
        # class code
        if (!is.null(isoData$class_code)) {
            isoData$class_code <-
                paste('class code: \"', isoData$class_code, '\"', sep = '')
        }
    }

    ### Prepare for plotting
    if (TRUE) {
        # melt categories
        isoDataMelt <-
            reshape2::melt(isoData,
                 id.vars = c(
                     'iso_ref', 'isoform_id', 'comparison', 'IF1', 'IF2'
                 ))
        isoDataMelt <-
            isoDataMelt[which(!is.na(isoDataMelt$value)), ]

        # melt IF
        colnames(isoDataMelt)[match(
            c('variable', 'value'), colnames(isoDataMelt)
        )] <- c('category', 'isoform_feature')
        isoDataMelt <-
            reshape2::melt(
                isoDataMelt,
                id.vars = c(
                    'iso_ref',
                    'isoform_id',
                    'comparison',
                    'category',
                    'isoform_feature'
                )
            )

        # massage category
        isoDataMelt$category <- as.character(isoDataMelt$category)

        isoDataMelt$category <- gsub(
            'PTC',
            'NMD Status',
            isoDataMelt$category
        )

        isoDataMelt$category <-
            gsub('codingPotential',
                 'Coding Potential',
                 isoDataMelt$category)

        isoDataMelt$category <-
            gsub('signal_peptide_identified',
                 'Signal Peptide',
                 isoDataMelt$category)

        isoDataMelt$category <-
            gsub('domains_identified',
                 'Protein Domains',
                 isoDataMelt$category)


        isoDataMelt$category <-
            gsub('IR','Intron Retention',isoDataMelt$category)
        isoDataMelt$category <-
            gsub('switchConsequences',
                 'Switch Consequence',
                 isoDataMelt$category)
        isoDataMelt$category <-
            gsub('class_code','Isoform Class',isoDataMelt$category)

        isoDataMelt$comparison2 <-
            paste(isoDataMelt$comparison, '\n(IF1 vs IF2)', sep = '')
    }

    ### Calculate statistics
    if (TRUE) {
        mySigTest <-
            plyr::ddply(
                isoDataMelt,
                .variables = c('comparison', 'category', 'isoform_feature'),
                .drop = TRUE,
                .fun = function(aDF) {
                    data1 <- aDF$value[which(aDF$variable == 'IF1')]
                    data2 <- aDF$value[which(aDF$variable == 'IF2')]

                    myTest <- suppressWarnings(wilcox.test(data1,
                                                           data2))

                    myResult <- data.frame(
                        n = length(data1),
                        medianIF1 = median(data1, na.rm = TRUE),
                        medianIF2 = median(data2, na.rm = TRUE)
                    )
                    myResult$medianDIF <-
                        myResult$medianIF2 - myResult$medianIF1

                    myResult$wilcoxPval <- myTest$p.value
                    myResult$ymax <- max(c(data1, data2))

                    return(myResult)
                }
            )
        mySigTest$ymax <- mySigTest$ymax * 1.05
        mySigTest$wilcoxQval <-
            p.adjust(mySigTest$wilcoxPval, method = 'fdr')
        mySigTest$significance <-
            sapply(mySigTest$wilcoxQval, function(x)
                evalSig(x, alphas))


        mySigTest <-
            plyr::ddply(
                mySigTest,
                .variables = 'category',
                .fun = function(aDF) {
                    aDF$isoform_feature <- factor(aDF$isoform_feature)
                    aDF$idNr <- as.numeric(aDF$isoform_feature)
                    return(aDF)
                }
            )

        mySigTest$comparison2 <-
            paste(mySigTest$comparison, '\n(IF1 vs IF2)', sep = '')

    }

    ### Plot result
    if (plot) {
        # start plot
        if (violinPlot) {
            p1 <- ggplot() +
                geom_violin(
                    data = isoDataMelt,
                    aes(
                        x = isoform_feature,
                        y = value,
                        fill = variable
                    ),
                    scale = 'area'
                ) +
                stat_summary(
                    data = isoDataMelt,
                    aes(
                        x = isoform_feature,
                        y = value,
                        fill = variable
                    ),
                    fun = medianQuartile,
                    geom = 'point',
                    position = position_dodge(width = 0.9),
                    size = 2
                )
        } else {
            p1 <- ggplot() +
                geom_boxplot(data = isoDataMelt, aes(
                    x = isoform_feature,
                    y = value,
                    fill = variable
                ))
        }

        # add significance
        p1 <- p1 +
            geom_text(
                data = mySigTest,
                aes(x = isoform_feature, y = ymax, label = significance),
                vjust = -0.2,
                size = localTheme$text$size * 0.3
            ) +
            geom_segment(data = mySigTest, aes(
                x = idNr - 0.25,
                xend = idNr + 0.25,
                y = ymax,
                yend = ymax
            ))

        # build rest of plot
        p1 <- p1 +
            facet_grid(comparison2 ~ category,
                       scales = 'free_x',
                       space = 'free_x') +
            localTheme +
            theme(strip.text.y = element_text(angle = 0)) +
            theme(axis.text.x = element_text(
                angle = -45,
                hjust = 0,
                vjust = 1
            )) +
            scale_fill_discrete(name = NULL) + theme(legend.position = "top") +
            labs(x = 'Isoform feature', y = 'Isoform Usage (IF)') +
            #coord_cartesian(ylim=c(0,1.25))
            coord_cartesian(ylim = c(0, 1.1 + max(c(
                0, 0.02 * (length(unique(
                    mySigTest$comparison2
                )) - 2)
            ))))
    }

    ### Return result
    if (returnResult) {
        if(plot) {
            print(p1)
        }

        mySigTest2 <-
            mySigTest[, c(
                'comparison',
                'category',
                'isoform_feature',
                'n',
                'medianIF1',
                'medianIF2',
                'medianDIF',
                'wilcoxPval',
                'wilcoxQval',
                'significance'
            )]
        mySigTest2 <-
            mySigTest2[order(mySigTest2$comparison,
                             mySigTest2$category,
                             mySigTest2$isoform_feature), ]
        return(mySigTest2)
    } else {
        if(plot) {
            return(p1)
        }
    }
}

# for backward compatability
extractGenomeWideAnalysis <- function(
    switchAnalyzeRlist,
    featureToExtract = 'isoformUsage',
    annotationToAnalyze = 'all',
    alpha=0.05,
    dIFcutoff = 0.1,
    log2FCcutoff = 1,
    violinPlot=TRUE,
    alphas=c(0.05, 0.001),
    localTheme=theme_bw(),
    plot=TRUE,
    returnResult=TRUE
) {
    extractConsequenceGenomeWide(
        switchAnalyzeRlist=switchAnalyzeRlist,
        featureToExtract=featureToExtract,
        annotationToAnalyze=annotationToAnalyze,
        alpha=alpha,
        dIFcutoff=dIFcutoff,
        log2FCcutoff=log2FCcutoff,
        violinPlot=violinPlot,
        alphas=alphas,
        localTheme=localTheme,
        plot=plot,
        returnResult=returnResult
    )
}
kvittingseerup/IsoformSwitchAnalyzeR documentation built on June 28, 2024, 5:41 p.m.