Query Uniprot and plot proteins

NAMESPACE

@@ -73,6 +73,7 @@ export(prepareFirehoseArchives)
 export(psichomics)
 export(queryEnsembl)
 export(queryFirehoseData)
+export(queryUniprot)
 export(readAnnotation)
 export(renameDuplicated)
 export(rm.null)

R/plots_information.R

@@ -24,25 +24,108 @@
 queryEnsembl <- function(path, query, grch37 = TRUE) {
     url <- paste0("http://", if(grch37) "grch37.", "rest.ensembl.org")
     resp <- GET(url, path=path, query=query)
-    stop_for_status(resp)
+    warn_for_status(resp)
     r <- content(resp, "text", encoding = "UTF8")
     return(fromJSON(r))
 }
 
+#' Query the Uniprot REST API
+#' 
+#' @param protein Character: protein to query
+#' @param format Character: format of the response
+#' 
+#' @importFrom httr GET
+#' @importFrom jsonlite fromJSON
+#' 
+#' @return Parsed response
+#' @export
+#' 
+#' @examples 
+#' query <- list(feature = "gene")
+#' queryUniprot(protein, format)
+queryUniprot <- function(protein, format="xml") {
+    url <- "http://www.uniprot.org"
+    path <- paste0("uniprot/", protein, ".", format)
+    resp <- GET(url, path=path)
+    warn_for_status(resp)
+    r <- content(resp, "text", encoding = "UTF8")
+    return(r)
+}
+
 infoUI <- function(id) {
     ns <- NS(id)
     tagList(
         uiOutput(ns("info")),
-        plotOutput(ns("transcripts"))
+        plotOutput(ns("plotTranscripts")),
+        uiOutput(ns("selectProtein")),
+        highchartOutput(ns("plotProtein"))
     )
 }
 
 noinfo <- function(output) {
-    output$info <- renderUI({
-        "No information available. Try another event."
-    })
+    output$info <- renderUI("No information available. Try another event.")
+    output$plotTranscripts <- renderPlot(NULL)
+}
+
+#' Plot a Uniprot protein
+#' 
+#' @param feature List of XML nodes: features of the Uniprot protein
+#' @param proteinLength Integer: protein length
+#' 
+#' @return Highcharter object
+proteinHighcharts <- function(feature, proteinLength) {
+    hc <- highchart() %>%
+        hc_chart(type="area", zoomType="x") %>%
+        hc_xAxis(title=list(text="Position (aminoacids)"), min=0,
+                 max=proteinLength) %>%
+        hc_yAxis(visible=FALSE) %>%
+        hc_tooltip(pointFormat="<b>{series.name} {point.id}</b>
+                   <br>{point.variant}{point.description}")
+
+    # The diverse types of features available
+    types <- unique(sapply(sapply(feature, xmlAttrs), "[[", 1))
 
-    output$transcripts <- renderPlot(NULL)
+    featureList <- NULL
+    # Reverse elements from features so the first ones (smaller Y) are above
+    for (feat in rev(feature)) {
+        attrs <- xmlAttrs(feat)
+        type  <- attrs[[1]]
+        description <- attrs[[2]]
+        id <- tryCatch(attrs[[3]], error=function(e) NULL)
+        location  <- feat[[match("location", names(feat))]]
+
+        # Get original and variant aminoacid
+        variant <- match("variation", names(feat))
+        if (!is.na(variant) && !is.null(variant)) {
+            original  <- xmlToList(feat[[match("original", names(feat))]])
+            variation <- xmlToList(feat[[variant]])
+            variant <- sprintf("%s>%s: ", original, variation)
+        } else {
+            variant <- NULL
+        }
+
+        start <- as.numeric(xmlAttrs(location[[1]]))
+        # If there's no stop position, simply sum 1 to the start position
+        stop  <- tryCatch(as.numeric(xmlAttrs(location[[2]])),
+                          error=function(e) {return(start+1)})
+        y     <- match(type, types)
+
+        # Create a list with two points based on this region
+        temp <- list(
+            NULL, 
+            list(x=start, y=y, description=description, id=id, variant=variant), 
+            list(x=stop,  y=y, description=description, id=id, variant=variant))
+
+        # Either make a new list or append to existing
+        if (is.null(featureList[type])) {
+            featureList[[type]] <- temp[2:3]
+        } else {
+            featureList[[type]] <- c(featureList[[type]], temp)
+        }
+    }
+    for (type in names(featureList))
+        hc <- hc %>% hc_add_series(name=type, data=featureList[[type]])
+    return(hc)
 }
 
 #' @importFrom Sushi plotGenes zoomsregion labelgenome
@@ -63,11 +146,12 @@ infoServer <- function(input, output, session) {
 
         species <- tolower(getSpecies())
         assembly <- getAssemblyVersion()
-        path <- paste0("lookup/symbol/", species, "/", gene)
         grch37 <- assembly == "hg19"
+
+        path <- paste0("lookup/symbol/", species, "/", gene)
         info <- queryEnsembl(path, list(expand=1), grch37=grch37)
 
-        output$transcripts <- renderPlot({
+        output$plotTranscripts <- renderPlot({
             transcripts <- data.frame()
 
             for (i in 1:nrow(info$Transcript)) {
@@ -104,10 +188,10 @@ infoServer <- function(input, output, session) {
                 h2(info$display_name, tags$small(info$id)),
                 sprintf("Species: %s (assembly %s)", species, assembly),
                 br(), sprintf("Chromosome %s: %s-%s (%s strand)",
-                        info$seq_region_name, 
-                        format(start, big.mark=",", scientific=FALSE),
-                        format(end, big.mark=",", scientific=FALSE),
-                        ifelse(info$strand == -1,"reverse", "forward")),
+                              info$seq_region_name, 
+                              format(start, big.mark=",", scientific=FALSE),
+                              format(end, big.mark=",", scientific=FALSE),
+                              ifelse(info$strand == -1,"reverse", "forward")),
                 br(), sprintf("%s (%s)", info$description,
                               info$biotype), br(),
                 tags$a("Ensembl", icon("external-link"), target="_blank",
@@ -124,6 +208,42 @@ infoServer <- function(input, output, session) {
                                "Zoom to splicing event"="event"))
             )
         })
+
+        output$selectProtein <- renderUI({
+            proteins <- info$Transcript$Translation$id
+            proteins <- proteins[!is.na(proteins)]
+            selectizeInput(ns("protein"), label="Select protein",
+                           choices=proteins)
+        })
+
+        observe({
+            # Convert from ENSEMBL to Uniprot/SWISSPROT
+            ensembl <- input$protein
+
+            if (is.null(ensembl)) return(NULL)
+            print("Looking for ENSEMBL protein in Uniprot...")
+            uniprot <- queryEnsembl(paste0("xrefs/id/", ensembl),
+                                    list("content-type"="application/json"), 
+                                    grch37=grch37)
+            uniprot <- uniprot[grepl("SWISSPROT", uniprot$dbname), ]
+
+            if (nrow(uniprot) == 0)
+                print("No protein from Uniprot :(")
+            else {
+                protein <- uniprot$primary_id[1]
+                resp <- queryUniprot(protein, "xml")
+
+                doc <- xmlTreeParse(resp)
+                root <- xmlRoot(doc)[[1]]
+                feature <- getNodeSet(root, "//feature")
+                proteinLength <- as.numeric(
+                    xmlAttrs(getNodeSet(root, "//sequence")[[1]])[[1]])
+
+                hc <- proteinHighcharts(feature, proteinLength)
+                output$plotProtein <- renderHighchart(hc)
+                print("Uniprot protein found")
+            }
+        })
     })
 }