R/methods-VariantFilteringResults.R

In VariantFiltering: Filtering of coding and non-coding genetic variants

##
## show methodss
##

## provide a concise show method for CompressedVRangesList objects
## that result from doing split(vr, sampleNames(vr)) on a VRanges object 'vr'
setMethod("show", signature(object="CompressedVRangesList"),
          function(object) {
            lo <- length(object)
            cat(classNameForDisplay(object), " of length ", lo, "\n",
                sep = "")
            if (!is.null(names(object)))
              cat(S4Vectors:::labeledLine("names", names(object)))
          })

setMethod("show", signature(object="VariantFilteringResults"),
          function(object) {
            cat("\nVariantFiltering results object\n")
            cat(sprintf("\n  Genome version(s):"))
            for (i in seq_along(param(object)$seqInfos))
              if (length(param(object)$seqInfos[[i]]) > 0)
                cat(sprintf(" %s(%s)", paste(unique(genome(param(object)$seqInfos[[i]])), collapse=","),
                            seqlevelsStyle(param(object)$seqInfos[[i]])[1]))

            sampleNames <- ifelse(length(samples(object)) <= 4, paste(samples(object), collapse=", "),
                                  paste(paste(head(samples(object), n=4), collapse=", "), "...", sep=", "))
            cat(sprintf("\n  Number of individuals: %d (%s)\n", length(samples(object)), sampleNames))

            if (inheritanceModel(object) != "any")
              cat(sprintf("  Variants segregate according to a(n) %s inheritance model\n", inheritanceModel(object)))
            else
              cat("  Variants are not filtered by inheritance model\n")
            if (length(param(object)@qualityFilterNames) > 0) {
              cat("  Quality filters\n")
              print(active(filters(object)[param(object)@qualityFilterNames]))
            }
            cat("  Functional annotation filters\n")
            print(active(filters(object))[setdiff(names(filters(object)), param(object)@qualityFilterNames)])
          })

setMethod("summary", signature(object="VariantFilteringResults"),
          function(object, method=c("SO", "SOfull", "bioc")) {

            method <- match.arg(method)

            fv <- filteredVariants(object)
            fv1 <- fv[[1]]
            nvars <- length(unique(fv1$VCFIDX))

            res <- data.frame()
            if (method == "bioc") {
              total <- length(fv1)
              vxloc <- table(fv1$LOCATION)
              vxcon <- table(fv1$CONSEQUENCE)
              vxloc <- table(unique(mcols(fv1)[, c("VCFIDX", "LOCATION")])$LOCATION)
              vxcon <- table(unique(mcols(fv1)[, c("VCFIDX", "CONSEQUENCE")])$CONSEQUENCE)
              res <- data.frame("BIOCID"=c(names(vxloc), names(vxcon)),
                                "Nr. Variants"=c(as.integer(vxloc), as.integer(vxcon)),
                                check.names=FALSE)
              res <- res[res[[2]] > 0, ]
              f <- res[[2]] / nvars
              nd <- max(-floor(log10(f[f > 0]))-1)
              res$"% Variants" <- round(100*f, digits=nd)
            } else if (method == "SO" || method == "SOfull") {
              if (!"SOterms" %in% names(cutoffs(object)))
                stop("There is no 'SOterms' entry in the list of cutoff values.")

              soterms <- cutoffs(object)$SOterms
              gSO <- sog(object)

              ## no SO terms specified or no active SOterms filter, imply no restriction
              if (is.null(soterms) || all(is.na(soterms)) || !active(filters(object))["SOterms"])
                soterms <- nodes(gSO)[sapply(nodeData(gSO, nodes(gSO), "vcfIdx"), length) > 0]
              else
                soterms <- .findSOIDs(soterms, gSO)

              description <- character(0)
              numvariants <- integer(0)
              if (length(soterms) > 0) {
                if (method == "SO") { ## show only given SO terms
                  soterms <- soterms[order(match(soterms, tsort(gSO)))] ## show SO terms in topological order
                  numvariants <- sapply(nodeData(gSO, soterms, "vcfIdx"),
                                        function(vcfidx, allvcfidx) length(unique(vcfidx[vcfidx %in% allvcfidx])),
                                        fv1$VCFIDX)
                  description <- unlist(nodeData(gSO, soterms, "label"))
                  res <- data.frame("SOID"=soterms[numvariants > 0],
                                    "Description"=description[numvariants > 0],
                                    "Nr. Variants"=numvariants[numvariants > 0],
                                    check.names=FALSE, row.names=NULL, stringsAsFactors=FALSE)
                } else { ## show hierarchy involving the given SO terms

                  asoterms <- unique(c(soterms, .ancestorsSO(param(object), soterms)))
                  avcfidx <- nodeData(gSO, asoterms, "vcfIdx")
                  asoterms <- names(avcfidx)[sapply(avcfidx, length) > 0]
                  soterms <- unique(c(soterms, asoterms))
                  subgSO <- sog(param(object))
                  nodeDataDefaults(subgSO, "vcfIdx") <- integer(0)
                  nodeData(subgSO, soterms, "vcfIdx") <- nodeData(gSO, soterms, "vcfIdx")

                  dsoterms <- unique(c(soterms, .descendantsSO(param(object), soterms)))
                  subgSO <- subGraph(dsoterms, subgSO)
                  to <- tsort(subgSO)
                  soterms <- character(0)
                  for (v in to) {
                    soterms <- c(soterms, v)
                    parentterms <- inEdges(subgSO)[[v]]
                    parentvcfidx <- unlist(lapply(parentterms, function(x) nodeData(subgSO, x, "vcfIdx")), use.names=FALSE)
                    nodeData(subgSO, v, "vcfIdx")[[v]] <- unique(c(nodeData(subgSO, v, "vcfIdx")[[1]], parentvcfidx))
                    vcfidx <- nodeData(subgSO, v, "vcfIdx")[[v]]
                    numvariants <- c(numvariants, length(unique(vcfidx[vcfidx %in% fv1$VCFIDX])))
                    description <- c(description, unlist(nodeData(subgSO, v, "label"), use.names=FALSE))
                  }
                  alld <- dijkstra.sp(g=as(t(as(subgSO, "matrix")), "graphNEL"), start=to[length(to)])$distances
                  res <- data.frame("SOID"=soterms,
                                    "Level"=alld[soterms],
                                    "Description"=description,
                                    "Nr. Variants"=numvariants,
                                    check.names=FALSE, row.names=NULL, stringsAsFactors=FALSE)
                }
              }

              if (nrow(res) > 0) {
                f <- res[["Nr. Variants"]] / nvars
                nd <- max(-floor(log10(f[f > 0]))-1)
                res$"% Variants" <- round(100*f, digits=nd)
              } else
                res$"% Variants" <- numeric(0)
            }

            res
          })


##
## getter and setter methods
##

setMethod("length", "VariantFilteringResults", function(x) length(x@variants))

setMethod("filters", signature(x="VariantFilteringResults"),
          function(x) {
            x@filters
          })

setReplaceMethod("filters", signature(x="VariantFilteringResults"),
                 function(x, value) {
                   x@filters <- value
                   x
                 })

setMethod("filtersMetadata", signature(x="VariantFilteringResults"),
          function (x) {
            x@filtersMetadata
          })

setMethod("cutoffs", signature(x="VariantFilteringResults"),
          function(x) {
            x@cutoffs
          })

setReplaceMethod("cutoffs", signature(x="VariantFilteringResults", value="logical"),
                 function(x, value) {
                   x@cutoffs <- value
                   x
                 })

setMethod("sortings", signature(x="VariantFilteringResults"),
          function(x) {
            x@sortings
          })

setMethod("softFilterMatrix", signature(x="VariantFilteringResults"),
          function(x) {
            softFilterMatrix(allVariants(x, groupBy="nothing"))
          })

setReplaceMethod("softFilterMatrix", signature(x="VariantFilteringResults"),
                 function(x, value) {
                   softFilterMatrix(x@variants) <- value
                   x
                 })

setMethod("sog", signature(x="VariantFilteringResults"),
          function(x, reverse=FALSE) {
            g <- x@gSO
            if (reverse) {
              grev <- as(as(t(as(as(g, "graphAM"), "matrix")), "graphAM"), "graphNEL")
              nodeDataDefaults(grev, "label") <- NA_character_
              allterms <- unlist(nodeData(g, nodes(g), "label"))
              nodeData(grev, nodes(grev), "label") <- allterms[nodes(grev)]
              g <- grev
            }

            g
          })

setMethod("samples", signature(object="VariantFilteringResults"),
          function(object) {
            object@activeSamples
          })

setReplaceMethod("samples", signature(object="VariantFilteringResults"),
                 function(object, value) {
                   if (inheritanceModel(object) != "unrelated individuals")
                     stop("Active samples can only be changed when no inheritance model is used (unrelated individuals).")

                   mask <- !value %in% param(object)$sampleNames
                   if (any(mask)) {
                     if (sum(mask) > 1)
                       stop(sprintf("%s are not valid sample names.", value[mask]))
                     else
                       stop(sprintf("%s is not a valid sample name.", value[mask]))
                   }

                   object@activeSamples <- value

                   object
                 })

setMethod("resetSamples", signature(object="VariantFilteringResults"),
          function(object) {
            object@activeSamples <- param(object)$sampleNames

            object
          })

setMethod("bamFiles", signature(object="VariantFilteringResults"),
          function(object) {
            object@bamViews
          })

setReplaceMethod("bamFiles", signature(object="VariantFilteringResults", value="BamViews"),
                 function(object, value) {
                   mask <- rownames(bamSamples(value)) %in% param(object)$sampleNames
                   if (!all(mask))
                     stop("Sample names do not match.")

                   object@bamViews <- value

                   object
                 })

setMethod("param", signature(x="VariantFilteringResults"),
          function(x) {
            x@inputParameters
          })

setMethod("inheritanceModel", signature(x="VariantFilteringResults"),
          function(x) {
            x@inheritanceModel
          })

setMethod("annoGroups", signature(x="VariantFilteringResults"),
          function(x) {
            x@annoGroups
          })

setMethod("filters", signature(x="VariantFilteringResults"),
          function(x) {
            x@filters
          })

setMethod("referenceGenome", signature=(x="VariantFilteringResults"),
          function(x) x@genomeDescription)

##
## methods for retrieving variants
##

## get all variants without applying any filter
setMethod("allVariants", signature(x="VariantFilteringResults"), 
          function(x, groupBy="sample") {
            vars <- x@variants
            if (groupBy[1] %in% "sample" && length(sampleNames(x@variants)) > 0) {
              f <- sampleNames(x@variants)
              if (all(is.na(f)))
                f <- rep("nosample", length(x))
              vars <- split(x@variants, f)
            } else if (groupBy[1] %in% colnames(mcols(x@variants)))
              vars <- split(x@variants, mcols(x@variants)[, groupBy])

            vars
          })

## get variants after applying all filters
setMethod("filteredVariants", signature(x="VariantFilteringResults"), 
          function(x, groupBy=c("sample", "nothing"),
                   unusedColumns.rm=FALSE) {
            groupBy <- match.arg(groupBy)

            if (length(filters(x)) > 0)
              x <- softFilter(x, filters(x))
            varsxsam <- allVariants(x)

            if (length(varsxsam) == 0)
              return(varsxsam)

            vars <- varsxsam[[1]]
            selcols <- names(active(filters(x)))[active(filters(x))] ## default VCF filters may or may not show up
            rowsMask <- apply(softFilterMatrix(vars)[, selcols, drop=FALSE], 1, all, na.rm=TRUE)
            if (!all(param(x)$sampleNames %in% samples(x)) && all(samples(x) %in% names(varsxsam))) { ## not all samples are active
              varsxsam <- varsxsam[samples(x)]

              ## discard variants that are not present in active samples
              gt <- sapply(varsxsam, function(x) x$GT)
              gt <- apply(gt, 1, paste, collapse="")
              gt <- gsub("/|\\|", "", gt)
              gt <- gsub(".", "0", gt, fixed=TRUE)
              gt <- gsub("0+", "0", gt)
              rowsMask <- rowsMask & !gt %in% "0"
            }

            ## if any of the 5' cryptic ss or 3' cryptic ss meet the cutoff, select the row
            ## mtNoSCORE5ss <- mtNoSCORE3ss <- NULL
            ## if (all(!is.na(param(x)$spliceSiteMatricesFilenames))) {
            ##   crypssMask <- rep(FALSE, length(vars))
            ##   if (is.na(minScore5ss(x)))
            ##     mtNoSCORE5ss <- grep("SCORE5ss", colnames(mcols(vars)))
            ##   else {
            ##     cryp5ssMask <- vars$SCORE5ssALT >= minScore5ss(x)
            ##     cryp5ssMask[is.na(cryp5ssMask)] <- FALSE
            ##     crypssMask <- crypssMask | cryp5ssMask
            ##   }
            ##   if (is.na(minScore3ss(x)))
            ##     mtNoSCORE3ss <- grep("SCORE3ss", colnames(mcols(vars)))
            ##   else {
            ##     cryp3ssMask <- vars$SCORE3ssALT >= minScore3ss(x)
            ##     cryp3ssMask[is.na(cryp3ssMask)] <- FALSE
            ##     crypssMask <- crypssMask | cryp3ssMask
            ##   }
            ##   ## if no filter on 5' and 3' cryptic ss is set, then select all rows
            ##   if (is.na(minScore5ss(x)) && is.na(minScore3ss(x)))
            ##     crypssMask <- rep(TRUE, length(vars))
            ##
            ##   rowsMask <- rowsMask & crypssMask
            ## }

            ## select variant annotations using the logical mask of the filters
            colsIdx <- 1:ncol(mcols(vars))
            ## if (unusedColumns.rm) ## remove data columns that are not used for filtering
            ##   colsIdx <- setdiff(colsIdx, c(mtNoMAF, mtNoMinPhastCons, mtNoMinPhylostratum, mtNoSCORE5ss, mtNoSCORE3ss))

            colsMask <- rep(FALSE, ncol(mcols(vars)))
            colsMask[colsIdx] <- TRUE

            vars <- unlist(varsxsam, use.names=FALSE)
            sampleNames(vars) <- Rle(names(varsxsam), elementNROWS(varsxsam))
            rowsMask <- rep(rowsMask, length(varsxsam))

            vars <- vars[rowsMask, colsMask]

            if (sortings(x)$criterion != "position")
              vars <- vars[order(mcols(vars)[[as.character(sortings(x)$criterion)]], decreasing=sortings(x)$decreasing), ]
            else if (sortings(x)$decreasing)
              vars <- rev(vars)

            if (groupBy[1] %in% "sample" && length(sampleNames(vars)) > 0)
              vars <- split(vars, sampleNames(vars))
            else if (groupBy[1] %in% colnames(mcols(vars)))
              vars <- split(vars, mcols(vars)[, groupBy])

            vars
          })

##
## shiny app to filter and visualize variants
##

setMethod("reportVariants", signature(vfResultsObj="VariantFilteringResults"),
          function(vfResultsObj, type=c("shiny", "csv", "tsv"), file=NULL, UCSCorg=NA_character_) {
              
  type <- match.arg(type)

  if (class(vfResultsObj) != "VariantFilteringResults")
    stop("Input argument 'vfResultsObj' should be an object of class 'VariantFilteringResults'.")

  if ((type == "csv" || type == "tsv") && is.null(file))
    stop("If type=\"csv\" or type=\"tsv\" then the input argument 'file' cannot be NULL.")

  if (type == "csv" || type == "tsv") {
    varsdf <- filteredVariants(vfResultsObj, groupBy="nothing")
    varsdf <- as.data.frame(DataFrame(CHR=seqnames(varsdf),
                                      POS=start(varsdf),
                                      SAMPLEID=sampleNames(varsdf),
                                      mcols(varsdf)))
    firstcols <- c("VARID", "dbSNP", "CHR", "POS", "SAMPLEID", "GT")
    varsdf <- varsdf[, c(firstcols, setdiff(colnames(varsdf), firstcols))]
    if (type == "csv")
      write.csv(varsdf, file=file, quote=FALSE, row.names=FALSE, col.names=TRUE)
    else
      write.table(varsdf, file, row.names=FALSE, col.names=TRUE, quote=FALSE, sep="\t")
   

    bn <- basename(file)
    message(sprintf("Written %s", bn))
    fullpath <- file
    if (dirname(fullpath) == ".")
      fullpath <- file.path(getwd(), bn)
    return(fullpath)
  }
  
  ## parameters for the UCSC genome browser. the organism name is derived from
  ## the 'organism()' getter for BSgenome objects when not specified with the
  ## function argument 'UCSCorg'

  if (is.na(UCSCorg)) {
    UCSCorg <- organism(referenceGenome(vfResultsObj))
    UCSCorg <- paste0(substr(UCSCorg, 1, 1), strsplit(UCSCorg, " ")[[1]][2])
  }
  genomeBuild <- providerVersion(referenceGenome(vfResultsObj))

  annotationObjClasses <- param(vfResultsObj)$otherAnnotationsClass
  mtMafDb <- match("MafDb", annotationObjClasses)
  mtPhastConsDb <- match("GScores", annotationObjClasses)
  mtGenePhylostrataDb <- match("GenePhylostrataDb", annotationObjClasses)
  phylostrata <- "Unavailable"
  if (!is.na(mtGenePhylostrataDb))
    phylostrata <- rev(genePhylostrata(get(param(vfResultsObj)$otherAnnotations[mtGenePhylostrataDb]))$Description)
  ## crypsplice <- param(vfResultsObj)$spliceSiteMatricesFilenames
  
  app <- list(ui=NULL, server=NULL)

  app$ui <- pageWithSidebar(
              headerPanel('R/BioC VariantFiltering Web App'),
              sidebarPanel(
                ## genome tab
                conditionalPanel(condition="input.tsp == 'genome'",
                                 checkboxInput('dbSNPpresentFlag',
                                               'Filter by presence in dbSNP', FALSE)),
                conditionalPanel(condition="input.tsp == 'genome' && input.dbSNPpresentFlag == true",
                                 selectInput("dbSNPpresent", "Present in dbSNP:",
                                             choices=c("Yes", "No"))),
                conditionalPanel(condition="input.tsp == 'genome'", helpText(strong("Restrict variants to the following types:"))),
                conditionalPanel(condition="input.tsp == 'genome'", checkboxInput("SNV", "SNV", TRUE)),
                conditionalPanel(condition="input.tsp == 'genome'", checkboxInput("Insertion", "Insertion", TRUE)),
                conditionalPanel(condition="input.tsp == 'genome'", checkboxInput("Deletion", "Deletion", TRUE)),
                conditionalPanel(condition="input.tsp == 'genome'", checkboxInput("MNV", "MNV", TRUE)),
                conditionalPanel(condition="input.tsp == 'genome'", checkboxInput("Delins", "Delins", TRUE)),
                ## transcript tab
                conditionalPanel(condition="input.tsp == 'transcript'", numericInput('minCUFC', 'Minimum Codon Usage Absolute log2-Fold Change:', 0.00)),
                conditionalPanel(condition="input.tsp == 'transcript'", helpText(strong("Restrict variants to the following locations:"))),
                conditionalPanel(condition="input.tsp == 'transcript'", checkboxInput("coding", "Coding", TRUE)),
                conditionalPanel(condition="input.tsp == 'transcript'", checkboxInput("fiveUTR", "5' UTR", TRUE)),
                conditionalPanel(condition="input.tsp == 'transcript'", checkboxInput("threeUTR", "3' UTR", TRUE)),
                conditionalPanel(condition="input.tsp == 'transcript'", checkboxInput("intron", "Intronic", TRUE)),
                conditionalPanel(condition="input.tsp == 'transcript'", checkboxInput("spliceSite", "Known splice site", TRUE)),
                conditionalPanel(condition="input.tsp == 'transcript'", checkboxInput("promoter", "Promoter", TRUE)),
                conditionalPanel(condition="input.tsp == 'transcript'", checkboxInput("intergenic", "Intergenic", TRUE)),
                
                ## gene tab
                conditionalPanel(condition="input.tsp == 'gene'",
                                 checkboxInput('OMIMpresentFlag',
                                               'Filter by presence in OMIM', FALSE)),
                conditionalPanel(condition="input.tsp == 'gene' && input.OMIMpresentFlag == true",
                                 selectInput("OMIMpresent", "Present in OMIM:",
                                             choices=c("Yes", "No"))),
                ## protein tab
                conditionalPanel(condition="input.tsp == 'protein'", helpText(strong("Restrict variants to the following consequences:"))),
                conditionalPanel(condition="input.tsp == 'protein'", checkboxInput("synonymous", "Synonymous", TRUE)),
                conditionalPanel(condition="input.tsp == 'protein'", checkboxInput("nonsynonymous", "Nonsynonymous", TRUE)),
                conditionalPanel(condition="input.tsp == 'protein'", checkboxInput("frameshift", "Frameshift", TRUE)),
                conditionalPanel(condition="input.tsp == 'protein'", checkboxInput("nonsense", "Nonsense", TRUE)),
                conditionalPanel(condition="input.tsp == 'protein'", checkboxInput("not tranlated", "Not translated", TRUE)),
                conditionalPanel(condition="input.tsp == 'protein'", selectInput("aaChangeType", "Amino acid change type:",
                                                                     choices=c("Any", "Radical", "Conservative"))),
                ## MAF tab
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('naMAF', 'Keep variants without MAF', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", numericInput('maxMAF', 'Maximum MAF:', 1.00)),
                conditionalPanel(condition="input.tsp == 'maf'", helpText("Note: the maximum MAF cutoff is applied on",
                                                                          "the following selected human populations:")),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AFKG', 'All MAF KG', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AMR_AFKG', 'AMR MAF KG', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('ASN_AFKG', 'ASN MAF KG', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AFR_AFKG', 'AFR MAF KG', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('EUR_AFKG', 'EUR MAF KG', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AFESP', 'All MAF ESP', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('EA_AFESP', 'EA MAF ESP', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AA_AFESP', 'AA MAF ESP', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AFExAC', 'All MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AFR_AFExAC', 'AFR MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('AMR_AFExAC', 'AMR MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('Adj_AFExAC', 'Adj MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('EAS_AFExAC', 'EAS MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('FIN_AFExAC', 'FIN MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('NFE_AFExAC', 'NFE MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('OTH_AFExAC', 'OTH MAF ExAC', TRUE)),
                conditionalPanel(condition="input.tsp == 'maf'", checkboxInput('SAS_AFExAC', 'SAS MAF ExAC', TRUE)),
                ## conservation tab
                conditionalPanel(condition="input.tsp == 'conservation'", tags$hr()),
                ## if (!is.na(mtPhastConsDb))
                conditionalPanel(condition="input.tsp == 'conservation'",
                                 checkboxInput('minPhastConsFlag',
                                               'Filter by phastCons score', FALSE)),
                conditionalPanel(condition="input.tsp == 'conservation' && input.minPhastConsFlag == true",
                                 numericInput('minPhastCons', 'Minimum phastCons score:', 0.00)),
                ## if (!is.na(mtGenePhylostrataDb))
                conditionalPanel(condition="input.tsp == 'conservation'",
                                 checkboxInput('minPhylostratumFlag', 'Filter by conserved gene phylostratum', FALSE)),
                conditionalPanel(condition="input.tsp == 'conservation' && input.minPhylostratumFlag == true",
                                 selectInput("minPhylostratum", "Minimum conserved gene phylostratum:",
                                             choices=phylostrata)),
                ## cryptic splice sites tab
                conditionalPanel(condition="input.tsp == 'cryp'", tags$hr()),
                conditionalPanel(condition="input.tsp == 'cryp'", checkboxInput('minScore5ssFlag', 'Filter by cryptic 5\'ss', FALSE)),
                conditionalPanel(condition="input.tsp == 'cryp' && input.minScore5ssFlag == true",
                                 numericInput('minScore5ss', 'Minimum cryptic 5\'ss score:', -99)),
                conditionalPanel(condition="input.tsp == 'cryp'", checkboxInput('minScore3ssFlag', 'Filter by cryptic 3\'ss', FALSE)),
                conditionalPanel(condition="input.tsp == 'cryp' && input.minScore3ssFlag == true",
                                 numericInput('minScore3ss', 'Minimum cryptic 3\'ss score:', -99)),
                tags$hr(),
                downloadButton('downloadData', 'Download Variants'),
                downloadButton('generateReport', 'Generate Report'),
                actionButton('closesavebutton', 'Save & Close')
              ),
              mainPanel(
                uiOutput('mytabs')
              )
            )

  app$server <- function(input, output) {

    observe({
      if (input$closesavebutton == 0)
        return()
      isolate({
        ## presence in dbSNP
        ## if (input$dbSNPpresentFlag)
        ##   dbSNPpresent(vfResultsObj) <- input$dbSNPpresent
        ## else
        ##   dbSNPpresent(vfResultsObj) <- NA_character_

        ## presence in OMIM
        ## if (input$OMIMpresentFlag)
        ##   OMIMpresent(vfResultsObj) <- input$OMIMpresent
        ## else
        ##   OMIMpresent(vfResultsObj) <- NA_character_

        ## type of variant
        ## varTypMask <- variantType(vfResultsObj)
        ## for (i in names(varTypMask))
        ##   varTypMask[i] <- input[[i]]
        ## variantType(vfResultsObj) <- varTypMask

        ## variant location
        ## locMask <- variantLocation(vfResultsObj)
        ## for (i in names(locMask))
        ##   locMask[i] <- input[[i]]
        ## variantLocation(vfResultsObj) <- locMask

        ## variant consequence
        conMask <- variantConsequence(vfResultsObj)
        for (i in names(conMask))
          conMask[i] <- input[[i]]
        variantConsequence(vfResultsObj) <- conMask

        ## type of amino acid change
        ## aaChangeType(vfResultsObj) <- input$aaChangeType

        ## minimum allele frequency
        if (!is.na(mtMafDb)) {
          mafMask <- MAFpop(vfResultsObj)
          for (i in names(mafMask)) ## unlisting a reactivevalues object does not work :(
            mafMask[i] <- input[[i]]
          MAFpop(vfResultsObj) <- mafMask
          maxMAF(vfResultsObj) <- as.numeric(input$maxMAF)
          naMAF(vfResultsObj) <- input$naMAF
        }

        ## nucleotide conservation
        if (!is.na(mtPhastConsDb)) {
          if (input$minPhastConsFlag)
            minPhastCons(vfResultsObj) <- input$minPhastCons
          else
            minPhastCons(vfResultsObj) <- NA_real_
        }

        ## gene conservation
        if (!is.na(mtGenePhylostrataDb)) {
          if (input$minPhylostratumFlag)
            minPhylostratum(vfResultsObj) <- input$minPhylostratum
          else
            minPhylostratum(vfResultsObj) <- NA_integer_
        }

        ## cryptic splice sites
        ## if (all(!is.na(crypsplice))) {
        ##   if (input$minScore5ssFlag)
        ##     minScore5ss(vfResultsObj) <- input$minScore5ss
        ##   else
        ##     minScore5ss(vfResultsObj) <- NA_real_
        ##
        ##   if (input$minScore3ssFlag)
        ##     minScore3ss(vfResultsObj) <- input$minScore3ss
        ##   else
        ##     minScore3ss(vfResultsObj) <- NA_real_
        ## }

        ## codon-usage fold-change
        minCUFC(vfResultsObj) <- as.numeric(input$minCUFC)

        stopApp(returnValue=vfResultsObj)
      })
    })

    filteredVariantsReact <- reactive({
      vdf <- data.frame()
      withProgress(message="Filtering variants", value=0, {
      incProgress(1/3, detail="setting filters ...")
      ## presence in dbSNP
      ## if (input$dbSNPpresentFlag)
      ##   dbSNPpresent(vfResultsObj) <- input$dbSNPpresent
      ## else
      ##   dbSNPpresent(vfResultsObj) <- NA_character_

      ## presence in OMIM
      ## if (input$OMIMpresentFlag)
      ##   OMIMpresent(vfResultsObj) <- input$OMIMpresent
      ## else
      ##   OMIMpresent(vfResultsObj) <- NA_character_

      ## type of variant
      ## varTypMask <- variantType(vfResultsObj)
      ## for (i in names(varTypMask))
      ##   varTypMask[i] <- input[[i]]
      ## variantType(vfResultsObj) <- varTypMask

      ## variant location
      ## locMask <- variantLocation(vfResultsObj)
      ## for (i in names(locMask))
      ##   locMask[i] <- input[[i]]
      ## variantLocation(vfResultsObj) <- locMask

      ## variant consequence
      ## conMask <- variantConsequence(vfResultsObj)
      ## for (i in names(conMask))
      ##   conMask[i] <- input[[i]]
      ## variantConsequence(vfResultsObj) <- conMask

      ## type of amino acid change
      ## aaChangeType(vfResultsObj) <- input$aaChangeType

      ## minimum allele frequency
      ## if (!is.na(mtMafDb)) {
      ##   mafMask <- MAFpop(vfResultsObj)
      ##   for (i in names(mafMask)) ## unlisting a reactivevalues object does not work :(
      ##     mafMask[i] <- input[[i]]
      ##   MAFpop(vfResultsObj) <- mafMask
      ##   maxMAF(vfResultsObj) <- as.numeric(input$maxMAF)
      ##   naMAF(vfResultsObj) <- input$naMAF
      ## }

      ## nucleotide conservation
      ## if (!is.na(mtPhastConsDb)) {
      ##   if (input$minPhastConsFlag)
      ##     minPhastCons(vfResultsObj) <- input$minPhastCons
      ##   else
      ##     minPhastCons(vfResultsObj) <- NA_real_
      ## }

      ## gene conservation
      ## if (!is.na(mtGenePhylostrataDb)) {
      ##   if (input$minPhylostratumFlag)
      ##     minPhylostratum(vfResultsObj) <- input$minPhylostratum
      ##   else
      ##     minPhylostratum(vfResultsObj) <- NA_integer_
      ## }

      ## cryptic splice sites
      ## if (all(!is.na(crypsplice))) {
      ##   if (input$minScore5ssFlag)
      ##     minScore5ss(vfResultsObj) <- as.numeric(input$minScore5ss)
      ##   else
      ##     minScore5ss(vfResultsObj) <- NA_real_
      ##
      ##   if (input$minScore3ssFlag)
      ##     minScore3ss(vfResultsObj) <- as.numeric(input$minScore3ss)
      ##   else
      ##     minScore3ss(vfResultsObj) <- NA_real_
      ## }

      ## codon-usage fold-change
      ## minCUFC(vfResultsObj) <- as.numeric(input$minCUFC)

      incProgress(2/3, detail="applying filters ...")
      fvxsam <- filteredVariants(vfResultsObj)
      incProgress(3/3, detail="retrieving variants ...")
      fv <- fvxsam[[1]]
      vdf <- DataFrame(VarID=fv$VARID, CHR=seqnames(fv),
                       POS=start(fv), POSITION=start(fv),
                       mcols(fv))
      vdf$REF <- ref(fv)
      vdf$ALT <- alt(fv)
      if (length(fv) > 1) {
        vdf$DP <- as.integer(round(rowMeans(sapply(fvxsam, totalDepth)), digits=0))
        vdf$REFDP <- as.integer(round(rowMeans(sapply(fvxsam, refDepth)), digits=0))
        vdf$ALTDP <- as.integer(round(rowMeans(sapply(fvxsam, altDepth)), digits=0))
      } else if (length(fv) == 1) {
        vdf$DP <- as.integer(round(mean(sapply(fvxsam, totalDepth)), digits=0))
        vdf$REFDP <- as.integer(round(mean(sapply(fvxsam, refDepth)), digits=0))
        vdf$ALTDP <- as.integer(round(mean(sapply(fvxsam, altDepth)), digits=0))
      } else
        vdf$DP <- vdf$REFDP <- vdf$ALTDP <- integer()
      vdf <- as.data.frame(vdf)
      
      if (nrow(vdf) > 0) {
        varlocs <- paste0("<a href=http://genome.ucsc.edu/cgi-bin/hgTracks?org=", UCSCorg,
                          "&db=", genomeBuild, 
                          "&position=", vdf$CHR, ":", vdf$POS, " target=\"ucsc\">",
                          vdf$CHR, ":", vdf$POS,
                          "</a>")
        tempOMIM <- rep(NA_character_, nrow(vdf))
        tempOMIM[!is.na(vdf$OMIM)] <- sapply(strsplit(vdf$OMIM[!is.na(vdf$OMIM)], " *, *"), function(x) {
          z <- paste0("<a href=http://www.omim.org/entry/", x, " target=\"omim\">", x, "</a>")
          a <- paste(z, collapse=", ")
          a
        })
      
        tempdbSNP <- rep(NA_character_, nrow(vdf))
        tempdbSNP[!is.na(vdf$dbSNP)] <- sapply(strsplit(gsub("rs", "", vdf$dbSNP[!is.na(vdf$dbSNP)]), " *, *"),
                                               function(x) {
          z <- paste0("<a href=http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=", x,
                      " target=\"dbsnp\">", paste0("rs", x), "</a>")
          a <- paste(z, collapse=", ")
          a
        })
      
        vdf[["POSITION"]] <- varlocs
        vdf[["OMIM"]] <- tempOMIM
        vdf[["dbSNP"]] <- tempdbSNP
      }

      vdf <- vdf[, -match(c("CHR", "POS"), colnames(vdf))]
      rownames(vdf) <- NULL

      }) ## withProgress
      vdf 
    })

    output$mytabs <- renderUI({
      tabPanelList <- list(tabPanel("Genome", tableOutput('tableGenome'), value="genome"),
                           tabPanel("Gene", tableOutput('tableGene'), value="gene"),
                           tabPanel("Transcript", tableOutput('tableTranscript'), value="transcript"),
                           tabPanel("Protein", htmlOutput('tableProtein'), value="protein"))

      if ("MafDb" %in% annotationObjClasses)
        tabPanelList[[length(tabPanelList)+1]] <- tabPanel("MAF", tableOutput('tableMAF'), value="maf")

      if ("GScores" %in% annotationObjClasses || "GenePhylostrataDb" %in% annotationObjClasses)
        tabPanelList[[length(tabPanelList)+1]] <- tabPanel("Conservation", tableOutput('tableConservation'), value="conservation")

      ## if (all(!is.na(param(vfResultsObj)$spliceSiteMatricesFilenames)))
      ##   tabPanelList[[length(tabPanelList)+1]] <- tabPanel("CrypSplice", tableOutput('tableCrypSplice'), value="cryp")

      tabPanelList[[length(tabPanelList)+1]] <- tabPanelAbout()

      tabPanelList[[length(tabPanelList)+1]] <- "tsp"

      names(tabPanelList) <- c(rep("", length(tabPanelList)-1), "id")

      do.call(tabsetPanel, tabPanelList)
    })

    output$tableGenome <- renderTable({
      filteredVariantsReact()[, c("VarID", "POSITION", "dbSNP", "TYPE", "HGVSg", "DP", "REFDP", "ALTDP")]
    }, NA.string="NA",  sanitize.text.function=function(x){x})

    output$tableGene <- renderTable({
      filteredVariantsReact()[, c("VarID", "POSITION", "GENE", "LOCATION", "HGVSc", "OMIM")]
    }, NA.string="NA",  sanitize.text.function=function(x){x})

    output$tableTranscript <- renderTable({
      filteredVariantsReact()[, c("VarID", "POSITION", "GENE", "TXID", "LOCATION", "LOCSTART", "LOCEND", "cDNALOC", "HGVSc", "CUREF", "CUALT", "CUFC")]
    }, NA.string="NA",  sanitize.text.function=function(x){x})

    output$tableProtein <- renderTable({
      selcols <- c("VarID", "GENE", "HGVSp", "CONSEQUENCE", "AAchange", "AAchangeType")
      if ("PolyPhenDb" %in% annotationObjClasses)
        selcols <- c(selcols, "PolyPhen2")
      if ("PROVEANDb" %in% annotationObjClasses)
        selcols <- c(selcols, "PROVEAN")

      filteredVariantsReact()[, selcols]
    }, NA.string="NA",  sanitize.text.function=function(x){x})

    output$tableMAF <- renderTable({
      selcols <- selcolsnames <- c("VarID", "dbSNP")
      if ("MafDb" %in% annotationObjClasses) {
        selcols <- c(selcols, "maxMAF", names(MAFpop(vfResultsObj))[MAFpop(vfResultsObj)])
        selcolsnames <- c(selcolsnames, "Max", gsub("_", " ", names(MAFpop(vfResultsObj))[MAFpop(vfResultsObj)]))
      }

      fv <- filteredVariantsReact()
      fv <- fv[, selcols]
      colnames(fv) <- selcolsnames

      fv
    }, NA.string="NA",  sanitize.text.function=function(x){x})

    output$tableConservation <- renderTable({
      selcols <- c("VarID", "POSITION", "GENE")
      if ("GScores" %in% annotationObjClasses && !is.na(mtPhastConsDb)) {
        cname <- type(get(param(vfResultsObj)$otherAnnotations[mtPhastConsDb]))
        selcols <- c(selcols, cname)
      }
      if ("GenePhylostrataDb" %in% annotationObjClasses)
        selcols <- c(selcols, "GenePhylostratum", "GenePhylostratumTaxID")

      fv <- filteredVariantsReact()[, selcols]

      if ("GenePhylostrataDb" %in% annotationObjClasses) {
        fv$GenePhylostratum[!is.na(fv$GenePhylostratum)] <- sprintf("<a href=\"http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=%s&lvl=5\" target='ncbitaxonomybrowser'>%s</a>", fv$GenePhylostratumTaxID[!is.na(fv$GenePhylostratum)], fv$GenePhylostratum[!is.na(fv$GenePhylostratum)])
        fv <- fv[, -match("GenePhylostratumTaxID", colnames(fv))]
      }

      fv
    }, NA.string="NA", sanitize.text.function=function(x){x})

    output$tableCrypSplice <- renderTable({
      selcols <- selcolsnames <- c("VarID", "POSITION")

      ## if (all(!is.na(param(vfResultsObj)$spliceSiteMatricesFilenames))) {
      ##   selcols <- c(selcols, "SCORE5ssREF", "SCORE5ssALT", "SCORE5ssPOS", "SCORE3ssREF", "SCORE3ssALT", "SCORE3ssPOS")
      ##   selcolsnames <- c(selcolsnames, "5'ss Ref", "5'ss Alt", "5'ss Pos", "3'ss Ref", "3'ss Alt", "3'ss Pos")
      ## }

      fv <- filteredVariantsReact()[, selcols]
      colnames(fv) <- selcolsnames

      fv
    }, NA.string="NA",  sanitize.text.function=function(x){x})

    ## this provides a tab-separated value file with the selected vfResultsObj
    output$downloadData <- downloadHandler(
      filename = function() { "vfResultsObj.tsv" },
      content = function(file) {
        ## presence in dbSNP
        ## if (input$dbSNPpresentFlag)
        ##   dbSNPpresent(vfResultsObj) <- input$dbSNPpresent
        ## else
        ##   dbSNPpresent(vfResultsObj) <- NA_character_

        ## presence in OMIM
        ## if (input$OMIMpresentFlag)
        ##   OMIMpresent(vfResultsObj) <- input$OMIMpresent
        ## else
        ##   OMIMpresent(vfResultsObj) <- NA_character_

        ## type of variant
        ## varTypMask <- variantType(vfResultsObj)
        ## for (i in names(varTypMask))
        ##   varTypMask[i] <- input[[i]]
        ## variantType(vfResultsObj) <- varTypMask
 
        ## variant location
        ## locMask <- variantLocation(vfResultsObj)
        ## for (i in names(locMask))
        ##   locMask[i] <- input[[i]]
        ## variantLocation(vfResultsObj) <- locMask

        ## variant consequence
        conMask <- variantConsequence(vfResultsObj)
        for (i in names(conMask))
          conMask[i] <- input[[i]]
        variantConsequence(vfResultsObj) <- conMask

        ## type of amino acid change
        ## aaChangeType(vfResultsObj) <- input$aaChangeType

        ## minimum allele frequency
        if (!is.na(mtMafDb)) {
          mafMask <- MAFpop(vfResultsObj)
          for (i in names(mafMask)) ## unlisting a reactivevalues object does not work :(
            mafMask[i] <- input[[i]]
          MAFpop(vfResultsObj) <- mafMask
          maxMAF(vfResultsObj) <- input$maxMAF
          naMAF(vfResultsObj) <- input$naMAF
        }

        ## nucleotide conservation
        if (!is.na(mtPhastConsDb)) {
          if (input$minPhastConsFlag)
            minPhastCons(vfResultsObj) <- input$minPhastCons
          else
            minPhastCons(vfResultsObj) <- NA_real_
        }

        ## gene conservation
        if (!is.na(mtGenePhylostrataDb)) {
          if (input$minPhylostratumFlag)
            minPhylostratum(vfResultsObj) <- input$minPhylostratum
          else
            minPhylostratum(vfResultsObj) <- NA_integer_
        }

        ## cryptic splice sites
        ## if (all(!is.na(crypsplice))) {
        ##   if (input$minScore5ssFlag)
        ##     minScore5ss(vfResultsObj) <- as.numeric(input$minScore5ss)
        ##   else
        ##     minScore5ss(vfResultsObj) <- NA_real_
        ##
        ##   if (input$minScore3ssFlag)
        ##     minScore3ss(vfResultsObj) <- as.numeric(input$minScore3ss)
        ##   else
        ##     minScore3ss(vfResultsObj) <- NA_real_
        ## }

        ## codon-usage fold-change
        minCUFC(vfResultsObj) <- as.numeric(input$minCUFC)

        ## apply filters
        fvxsam <- filteredVariants(vfResultsObj)
        fv <- fvxsam[[1]]

        ## build data frame with the annotated and filtered variants
        vdf <- DataFrame(VarID=fv$VARID, CHR=seqnames(fv),
                         POS=start(fv), POSITION=start(fv),
                         mcols(fv))
        vdf$REF <- ref(fv)
        vdf$ALT <- alt(fv)
        vdf <- as.data.frame(vdf)

        write.table(vdf, file, row.names=FALSE, col.names=TRUE, quote=FALSE, sep="\t")
      },
      contentType = "text/plain"
    )

    ## this should generate a report to facilitate reproducing the analysis with R code
    output$generateReport <- downloadHandler(
      filename = function() { "report.txt" },
      content = function(file) {
        ## presence in dbSNP
        ## if (input$dbSNPpresentFlag)
        ##   dbSNPpresent(vfResultsObj) <- input$dbSNPpresent
        ## else
        ##   dbSNPpresent(vfResultsObj) <- NA_character_

        ## presence in OMIM
        ## if (input$OMIMpresentFlag)
        ##   OMIMpresent(vfResultsObj) <- input$OMIMpresent
        ## else
        ##   OMIMpresent(vfResultsObj) <- NA_character_

        ## type of variant
        ## varTypMask <- variantType(vfResultsObj)
        ## for (i in names(varTypMask))
        ##   varTypMask[i] <- input[[i]]
        ## variantType(vfResultsObj) <- varTypMask

        ## variant location
        ## locMask <- variantLocation(vfResultsObj)
        ## for (i in names(locMask))
        ##   locMask[i] <- input[[i]]
        ## variantLocation(vfResultsObj) <- locMask

        ## variant consequence
        conMask <- variantConsequence(vfResultsObj)
        for (i in names(conMask))
          conMask[i] <- input[[i]]
        variantConsequence(vfResultsObj) <- conMask

        ## type of amino acid change
        ## aaChangeType(vfResultsObj) <- input$aaChangeType

        ## minimum allele frequency
        if (!is.na(mtMafDb)) {
          mafMask <- MAFpop(vfResultsObj)
          for (i in names(mafMask)) ## unlisting a reactivevalues object does not work :(
            mafMask[i] <- input[[i]]
          MAFpop(vfResultsObj) <- mafMask
          maxMAF(vfResultsObj) <- input$maxMAF
          naMAF(vfResultsObj) <- input$naMAF
        }

        ## nucleotide conservation
        if (!is.na(mtPhastConsDb)) {
          if (input$minPhastConsFlag)
            minPhastCons(vfResultsObj) <- input$minPhastCons
          else
            minPhastCons(vfResultsObj) <- NA_real_
        }

        ## gene conservation
        if (!is.na(mtGenePhylostrataDb)) {
          if (input$minPhylostratumFlag)
            minPhylostratum(vfResultsObj) <- input$minPhylostratum
          else
            minPhylostratum(vfResultsObj) <- NA_integer_
        }

        ## cryptic splice sites
        ## if (all(!is.na(crypsplice))) {
        ##   if (input$minScore5ssFlag)
        ##     minScore5ss(vfResultsObj) <- as.numeric(input$minScore5ss)
        ##   else
        ##     minScore5ss(vfResultsObj) <- NA_real_
        ##
        ##   if (input$minScore3ssFlag)
        ##     minScore3ss(vfResultsObj) <- as.numeric(input$minScore3ss)
        ##   else
        ##     minScore3ss(vfResultsObj) <- NA_real_
        ## }

        ## codon-usage fold-change
        minCUFC(vfResultsObj) <- as.numeric(input$minCUFC)

        sink(file)
        cat("R script and output to get the resulting filtered variants with VariantFiltering\n")
        cat("================================================================================\n\n")

        cat("> library(VariantFiltering)\n\n")
        cat(sprintf("> %s <- %s\n", gettext(vfResultsObj@callObj)[2], param(vfResultsObj)$callStr))
        cat(sprintf("> res <- %s\n", deparse(vfResultsObj@callObj)))
        ## if (!is.na(dbSNPpresent(vfResultsObj)))
        ##   cat(sprintf("> dbSNPpresent(res) <- \"%s\"\n", dbSNPpresent(vfResultsObj)))
        ## if (!is.na(OMIMpresent(vfResultsObj)))
        ##   cat(sprintf("> OMIMpresent(res) <- \"%s\"\n", OMIMpresent(vfResultsObj)))
        ## if (!all(variantType(vfResultsObj)))
        ##     cat(sprintf("> variantType(res) <- c(%s)\n", paste(paste(names(variantType(vfResultsObj)), as.character(variantType(vfResultsObj)), sep="="), collapse=", ")))
        ## if (!all(variantLocation(vfResultsObj)))
        ##     cat(sprintf("> variantLocation(res) <- c(%s)\n", paste(paste(names(variantLocation(vfResultsObj)), as.character(variantLocation(vfResultsObj)), sep="="), collapse=", ")))
        if (!all(variantConsequence(vfResultsObj)))
            cat(sprintf("> variantConsequence(res) <- c(%s)\n", paste(paste(names(variantConsequence(vfResultsObj)), as.character(variantConsequence(vfResultsObj)), sep="="), collapse=", ")))
        ## if (aaChangeType(vfResultsObj) != "Any")
        ##   cat(sprintf("> aaChangeType(res) <- \"%s\"\n", aaChangeType(vfResultsObj)))
        if (minCUFC(vfResultsObj) > 0)
          cat(sprintf("> minCUFC(res) <- %.2f\n", minCUFC(vfResultsObj)))
        if (!is.na(mtMafDb)) {
          if (!all(MAFpop(vfResultsObj))) ## default values need not to be set
            cat(sprintf("> MAFpop(res) <- c(%s)\n", paste(paste(names(MAFpop(vfResultsObj)), as.character(MAFpop(vfResultsObj)), sep="="), collapse=", ")))
          if (!naMAF(vfResultsObj)) ## default values need not to be set
            cat(sprintf("> naMAF(res) <- %s\n", ifelse(input$naMAF, "TRUE", "FALSE")))
          cat(sprintf("> maxMAF(res) <- %.3f\n", input$maxMAF))
        }
        if (!is.na(mtPhastConsDb)) {
          if (!is.na(minPhastCons(vfResultsObj)))
            cat(sprintf("> minPhastCons(res) <- %.1f\n", minPhastCons(vfResultsObj)))
        }
        if (!is.na(mtGenePhylostrataDb)) {
          if (!is.na(minPhylostratum(vfResultsObj)))
            cat(sprintf("> minPhylostratum(res) <- %d\n", minPhylostratum(vfResultsObj)))
        }
        ## if (all(!is.na(crypsplice))) {
        ##   if (!is.na(minScore5ss(vfResultsObj)))
        ##     cat(sprintf("> minScore5ss(res) <- %.2f\n", minScore5ss(vfResultsObj)))
        ##   if (!is.na(minScore3ss(vfResultsObj)))
        ##     cat(sprintf("> minScore3ss(res) <- %.2f\n", minScore3ss(vfResultsObj)))
        ## }
        cat("> res\n")
        print(vfResultsObj)
        cat("> reportVariants(res, type=\"tsv\", file=\"variants.tsv\")\n")
        cat("> sessionInfo()\n")
        print(sessionInfo())
        sink()
      },
      contentType = "text/plain"
    )

  }

  runApp(app)
})