######################
#####For Gene and TE subfamily expression shiny app####
#Load source data OUTSIDE  of ui as this only needs to be run once
libs <- .libPaths()
.libPaths(c("Rlibs", libs))
library(shiny)
library(RColorBrewer)
library(gplots)

load("Brainspan_gene_merged_stages.RData")
gexp<-tmp2
load("Brainspan_TE_subfam_merged_stages.RData")
texp<-tmp2
#Source helper functions -----
source("helpers.R")

#Shiny app ----
ui <- fluidPage(
  titlePanel("Brain TExplorer: Gene and TE Expression in the Developing Human Brain"),
  h4("If you are using this app please cite: Playfoot et al., 2021"),
  sidebarLayout(
    sidebarPanel(
      helpText("Type human gene name and/or TE subfamily name (case sensitive)"),
      textInput("gene", label="Gene Name", value="TRIM28"),
      selectInput("te", label="TE Subfamily Name", choices=texp[[1]]$ID,selected="LTR7C"),
      selectInput("region",label="Region for correlation",choices= c("AMY", "CB",  "DFC", "HIP", "ITC", "MFC", "THA", "OFC", "IPC", "A1C", "M1C", "S1C", "V1C", "STR", "STC", "VFC")),
      submitButton("Submit"),
      helpText("Note: TE subfamily expression is determined from multimapped RNA-seq reads and does not exclude reads resulting from 'pervasive transcription'. Only data from Brainspan is shown."),
      helpText("Stage Info: 2B (10-12pcw), 3A (13-15pcw), 3B-5 (16-38pcw), 6-9 (birth - 11 years), 10 (12-19 years), 11 (20-60+ years)."),
      helpText("Region Info: AMY (amygdala), CB (cerebellar cortex), DFC (dorsolateral prefrontal cortex), HIP (hippocampus), ITC (inferior temporal cortex), MFC (medial prefrontal cortex), THA (mediodorsal nucleus of the thalamus), OFC (orbital prefrontal cortex), IPC (posterior inferior parietal cortex), A1C (primary auditory (A1) cortex), M1C (primary motor (M1) cortex), S1C (primary somatosensory (S1) cortex), V1C (primary visual (V1) cortex), STR (striatum), STC (superior temporal cortex), VFC (ventrolateral prefrontal cortex)."),
    ),
    mainPanel(
      helpText("Note: Heatmap scale is given in Log2CPM and changes for each query. Sample information, unmerged stage and regional expression data in tabular format can be found in Playfoot et al., 2021 Supplemental Materials."),
      downloadButton("downloadgene", "Download Gene Plot"), 
      plotOutput("gene"),
      downloadButton("downloadte", "Download TE Plot"),
      plotOutput("te"),
      downloadButton("downloadregion", "Download Correlation Plot"),
      plotOutput("region")
    )
  )
)
# Define server logic ----
server <- function(input, output) {
  geneplot<-function(){
    ssubfunc<-function(gexp){subset(gexp,gexp$symbol==input$gene)}
    gexp<-lapply(gexp, ssubfunc)
    
    #Collapse all lists into dataframe
    stam<-do.call(rbind,gexp)
    rownames(stam)<-c("AMY", "CB",  "DFC", "HIP", "ITC", "MFC", "THA", "OFC", "IPC", "A1C", "M1C", "S1C", "V1C", "STR", "STC", "VFC")
    stamfin<-stam[,c(6:11)]
    #Change CPM values to log2CPM
    stamfin<-log2(stamfin)
    #Replace infinite value with NA
    stamfin[stamfin=="-Inf"]<-NA
    #Generating heatmap for gene of interest
    stamfin<-as.matrix(stamfin)
    heatmap.2(stamfin,Colv=NA,Rowv=NA,scale="none",srtCol=45,margins=c(8,4),trace="none",main=input$gene,col=rev(brewer.pal(11, "RdYlBu")),density.info="none")
  }
  output$gene<-renderPlot({
     geneplot()
    })
  teplot<-function(){
    tsubfunc<-function(texp){subset(texp,texp$ID==input$te)}
    texp<-lapply(texp, tsubfunc)
    
    #Collapse all lists into dataframe
    tstam<-do.call(rbind,texp)
    rownames(tstam)<-c("AMY", "CB",  "DFC", "HIP", "ITC", "MFC", "THA", "OFC", "IPC", "A1C", "M1C", "S1C", "V1C", "STR", "STC", "VFC")
    tstamfin<-tstam[,c(3:8)]
    #Change CPM values to log2CPM
    tstamfin<-log2(tstamfin)
    tstamfin[tstamfin=="-Inf"]<-NA
    #Generating heatmap for gene of interest
    tstamfin<-as.matrix(tstamfin)
    heatmap.2(tstamfin,Colv=NA,Rowv=NA,scale="none",trace="none",srtCol=45,margins=c(8,4),main=input$te,col=rev(brewer.pal(11, "RdYlBu")),density.info="none")
}
  output$te<-renderPlot({
    teplot()
  })
  
  regionplot<-function(){
    #Data prep. for selected gene
    ssubfunc<-function(gexp){subset(gexp,gexp$symbol==input$gene)}
    gexp<-lapply(gexp, ssubfunc)
    stam<-do.call(rbind,gexp)
    rownames(stam)<-c("AMY", "CB",  "DFC", "HIP", "ITC", "MFC", "THA", "OFC", "IPC", "A1C", "M1C", "S1C", "V1C", "STR", "STC", "VFC")
    stamfin<-stam[,c(6:11)]
    #stamfin<-log2(stamfin)
    #Replace infinite value with NA
    stamfin[stamfin=="-Inf"]<-NA
    cnames<-colnames(stamfin)
    
    #Data prep. for selected TE
    tsubfunc<-function(texp){subset(texp,texp$ID==input$te)}
    texp<-lapply(texp, tsubfunc)
    
    #Collapse all lists into dataframe
    tstam<-do.call(rbind,texp)
    rownames(tstam)<-c("AMY", "CB",  "DFC", "HIP", "ITC", "MFC", "THA", "OFC", "IPC", "A1C", "M1C", "S1C", "V1C", "STR", "STC", "VFC")
    tstamfin<-tstam[,c(3:8)]
    #Change CPM values to log2CPM
    #tstamfin<-log2(tstamfin)
    tstamfin[tstamfin=="-Inf"]<-NA
    #Get pearson correlation with p-value for each query
    ctest<-cor.test(as.numeric(stamfin[input$region,]),as.numeric(tstamfin[input$region,]),use="pairwise.complete.obs")
    pval<-round(ctest$p.value,digits=4)
    est<-round(ctest$estimate,digits=2)
    
    #Plot line plots for specific regions
    par(mar=c(10.1,5.1,5.1,5.1))
    plot(as.numeric(tstamfin[input$region,]), type="o",xlab=NA,ylab="SubFamily Expression (CPM)",col.lab="deepskyblue4",cex.lab=1.5,cex.axis=1.5,lty=1,lwd=2, las=2,xaxt="n",col="deepskyblue4")
    axis(1, at=1:6,labels=cnames, cex.lab=1.5, cex.axis=1.5,cex.sub=1.5,las=2)
    par(new=TRUE)
    plot(as.numeric(stamfin[input$region,]), type="o",ylab=NA,xlab=NA,lty=1,lwd=2,xaxt="n",yaxt="n",col="firebrick")
    axis(4,cex.lab=1.5,cex.axis=1.5,cex.sub=1.5,las=2)
    mtext("Gene Expression (CPM)", col="firebrick",side=4, line=4, cex=1.5,cex.lab=1.5,cex.axis=1.5,cex.sub=1.5,las=0)
    #mtext("Developmental Stage", side=1, line=5, cex=1.5,cex.lab=1.5,cex.axis=1.5,cex.sub=1.5,las=0)
    abline(v=4)
    title(paste0(input$gene," vs ",input$te," in ",input$region," Expression Corr=",round(ctest$estimate,digits=2), " p-value=",round(ctest$p.value,digits=4)))
  }
  output$region<-renderPlot({
    regionplot()
    })
  
  ###Define plots again for download.
  # downloadHandler contains 2 arguments as functions, namely filename, content
  output$downloadgene <- downloadHandler(
    filename =  function() {
      paste(input$gene,"_heatmap.pdf",sep="")
    },
    # content is a function with argument file. content writes the plot to the device
    content = function(file) {
      pdf(file) # open the pdf device
      geneplot()
      dev.off()  # turn the device off
      
    } 
  )
  output$downloadte <- downloadHandler(
    filename =  function() {
      paste(input$te,"_heatmap.pdf",sep="")
    },
    # content is a function with argument file. content writes the plot to the device
    content = function(file) {
      pdf(file) # open the pdf device
      teplot()
      dev.off()  # turn the device off
      
    } 
  )  
  output$downloadregion <- downloadHandler(
    filename =  function() {
      paste(input$gene,input$te,"in",input$region,"lineplot.pdf",sep="_")
    },
    # content is a function with argument file. content writes the plot to the device
    content = function(file) {
      pdf(file) # open the pdf device
      regionplot()
      dev.off()  # turn the device off
      
    } 
  )  
  
}

# Run the app ----
shinyApp(ui, server)

