################This files documents the scripts for AluAluCNVpredictor
################examples of input files can be accessed at https://github.com/BCM-Lupskilab/AluAluCNVpredictor

####server.R
library(shiny)
source("prediction.R")
source("prediction_hg38.R")

riskscore <- read.delim("data/Gene_level_prediction.txt",header = T)

shinyServer(function(input, output){

##query the gene-level prediction results
  allAlu <- reactive({
          as.character(riskscore[riskscore$genename==toupper(input$geneID),2])
  	}) 
  CNVAlu <- reactive({
          as.character(riskscore[riskscore$genename==toupper(input$geneID),3])
  	}) 
  MIMID<- reactive({
            as.character(riskscore[riskscore$genename==toupper(input$geneID),8])
          })     
  MIMscore <- reactive({
          round(as.numeric(riskscore[riskscore$genename==toupper(input$geneID),6]),digits = 3)
  	}) 
  MIMrank <- reactive({
            as.numeric(riskscore[riskscore$genename==toupper(input$geneID),5])
          })
  allscore <- reactive({
          round(as.numeric(riskscore[riskscore$genename==toupper(input$geneID),4]), digits= 3)
  	})
  allrank <- reactive({
          as.numeric(riskscore[riskscore$genename==toupper(input$geneID),11])
  	})

  output$textallAlu <- renderText({ 
          paste("Count of total intersecting Alu pairs: ", allAlu(), sep="") 
  	}) 
  output$textCNVAlu <- renderText({ 
          paste("Count of predicted CNV-Alu pairs: ", CNVAlu(), sep="") 
  	}) 
  output$textMIMID <- renderText({ 
            paste("MIM number: ", MIMID(), sep="") 
          })      
  output$textMIMscore <- renderText({ 
            paste("AAMR risk score for MIM genes: ", MIMscore(), sep="") 
          })      
  output$textMIMrnk <- renderText({
            paste("Rank in MIM genes: ",MIMrank(), "/12,074", sep="")
          })      
  output$textallscore <- renderText({ 
          paste("AAMR risk score for available RefSeq genes: ", allscore(), sep="") 
  	})         
  output$textallrank <- renderText({
          paste("Rank in available RefSeq genes: ",allrank(), "/23,637", sep="")
  	})    
  output$hist1 <- renderPlot({
          par(mar=c(2.5,2.5,1,1)+0.1, mgp = c(1.5, 0.5, 0))
          hist(unique(riskscore[,c(1,6)])[,2],breaks=100, xlab="AAMR risk scores for OMIM genes", main=NULL)
          abline(v=MIMscore(),col="red")
	})      
  output$hist2 <- renderPlot({
          par(mar=c(2.5,2.5,1,1)+0.1, mgp = c(1.5, 0.5, 0))
          hist(unique(riskscore[,c(1,4)])[,2],breaks=100, xlab="AAMR risk scores for RefSeq genes", main=NULL)
          abline(v=allscore(),col="red")
          })      
       
 ##query the interval-level prediction results
  my_data <- eventReactive(input$go,{
          if(input$hg=="GRCh37_hg19"){
          		df <-  prediction_hg19(input$GI5,input$GI3)
          	}else {
          		df <- prediction_hg38(input$GI5,input$GI3)
          		}
          })
  output$table1 <- renderTable({
	my_data()       
     	}, include.rownames = FALSE)
  })

####ui.R

library(shiny)
shinyUI(fluidPage(
  h1(“AluAluCNVpredictor", style="color:darkcyan"),  #darkcyan
  p("Lupski lab",style="color:dimgray;font-size: 20px"),
  hr(),
  p("Gene-level prediction",style="color:darkcyan;font-size: 15px"),
  wellPanel(fluidRow(
          column(2,textInput('geneID', 'Gene Name', value = "")),
          column(3,br(),
                 textOutput('textallAlu'),
                 textOutput('textCNVAlu'),
                 br(),
                 textOutput('textMIMID'),
                 textOutput('textMIMscore'),
                 textOutput('textMIMrnk'),
                 br(),
                 textOutput('textallscore'),
                 textOutput('textallrank')),
          column(3, plotOutput('hist1', height="200px")),
          column(3, plotOutput('hist2', height="200px")))),
   p("Interval-level prediction",style="color:darkcyan;font-size: 15px"),
   sidebarLayout(
          sidebarPanel(
                 selectInput('hg',"Reference Genome Assembly:", 
                                    choices=c("GRCh37_hg19","GRCh38_hg38")),
                 helpText("Query any predicted CNV-Alu pairs intersecting your regions of interest"),
                 helpText("format:chr:start-end"),
                 tags$style(type="text/css", "textarea {width:100%}") ,
                 tags$textarea(id="GI5", rows=5, placeholder = c("Enter the 5' Genomic Interval:")),
                 tags$style(type="text/css", "textarea {width:100%}") ,
                 tags$textarea(id="GI3", rows=5, placeholder = c("Enter the 3' Genomic Interval:")),
                 actionButton('go', 'Search'),width = 3) ,
          mainPanel(
                 tableOutput('table1')  )
         )))

#####prediction.R
prediction_hg19 <- function(up,down){
  setwd("/var/shiny-server/www/xiaofei/Alu_app")
  up2 <-  unlist(strsplit(gsub(" ", "", up, fixed = TRUE),"\n|:|-"))
  up2 <-  as.data.frame(matrix(up2, ncol =3, byrow = T))
  up2[,2]<- as.numeric(as.character(up2[,2]))-1
  up2[,4] <- 1:nrow(up2)
  down2 <- unlist(strsplit(gsub(" ", "", down, fixed = TRUE),"\n|:|-"))
  down2 <-  as.data.frame(matrix(down2, ncol =3, byrow = T))
  down2[,2]<- as.numeric(as.character(down2[,2]))-1
  down2[,4] <- 1:nrow(down2)
  write.table(up2, "data/test_up", col.names = F,row.names = F, quote = F, sep = "\t")
  write.table(down2, "data/test_down", col.names = F,row.names = F, quote = F, sep = "\t")
  
  command <- c("intersectBed -wa -wb -a data/predictedTRUE_up -b data/test_up|cut -f4,8 > data/PredictedAlupair")
  try(system(command))
  command <- c("intersectBed -wa -wb -a data/predictedTRUE_down -b data/test_down|cut -f4,8 >> data/PredictedAlupair")
  try(system(command))
  command <- c("cat data/PredictedAlupair|sort |uniq -c|grep  '^ *2 ' >data/PredictInPair")
  try(system(command))
  
  if(file.info("data/PredictInPair")$size>0){
  	pre_rslt <- read.table("data/PredictInPair", sep = "\t",header = F)
 	pre_rslt[,1] <- gsub(".*\\s","",pre_rslt[,1])
  	write.table(pre_rslt,"data/temp", quote = F,col.names = F,row.names = F,sep = "\t")
  	command <- c("sort -n -k1 data/temp >data/temp6 ")
  	try(system(command))
  	command <- c("awk 'NR==FNR {h[$1]; next} $4 in h' data/temp data/predictedTRUE_up |sort -n -k4>data/temp2")
  	try(system(command))
  	command <- c("awk 'NR==FNR {h[$1]; next} $4 in h' data/temp data/predictedTRUE_down |sort -n -k4>data/temp3")
  	try(system(command))
  	command <- c("join -1 1 -2 4 -a 1 data/temp6 data/temp2 >data/temp4")
  	try(system(command))
  	command <- c('join -1 1 -2 4 -a 1 data/temp4 data/temp3 >data/temp5')
 	try(system(command))
  	df <- read.table("data/temp5", header=F, sep = " " )
  	df2 <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
  	df2[,1]<- as.data.frame(df[,2])
  	df2[,2] <- as.data.frame(paste(df[,3],":",df[,4]+1,"-",df[,5], " | ", df[,6],":",df[,7]+1,"-",df[,8],sep=""))
 	colnames(df2) <- c("No.","Predicted_Alu_pair")
  	df3 <- data.frame(up2[,4],as.data.frame(paste(up2[,1],":",up2[,2]+1,"-",up2[,3], sep="")),
                    as.data.frame(paste(down2[,1],":",down2[,2]+1,"-",down2[,3], sep="")) )
  	colnames(df3) <-c("No.","5'_Genomic_Interval","3'_Genomic_Interval")
  	df4 <- merge(df3,df2, all = T)
 	command <- c("rm data/temp* data/test* data/Predict*")
  	try(system(command))
  	return(df4)
  }else{
          df3 <- data.frame(up2[,4],as.data.frame(paste(up2[,1],":",up2[,2]+1,"-",up2[,3], sep="")),
                            as.data.frame(paste(down2[,1],":",down2[,2]+1,"-",down2[,3], sep="")), NA )
          colnames(df3) <-c("No.","5'_Genomic_Interval","3'_Genomic_Interval", "Predicted_Alu_pair")
          command <- c("rm data/test* data/Predict*")
          try(system(command))
          return(df3)
  }
 }

####prediction_hg38.R
prediction_hg38 <- function(up,down){
        setwd("/var/shiny-server/www/xiaofei/Alu_app")
        up2 <-  unlist(strsplit(gsub(" ", "", up, fixed = TRUE),"\n|:|-"))
        up2 <-  as.data.frame(matrix(up2, ncol =3, byrow = T))
        up2[,2]<- as.numeric(as.character(up2[,2]))-1
        up2[,4] <- 1:nrow(up2)
        down2 <- unlist(strsplit(gsub(" ", "", down, fixed = TRUE),"\n|:|-"))
        down2 <-  as.data.frame(matrix(down2, ncol =3, byrow = T))
        down2[,2]<- as.numeric(as.character(down2[,2]))-1
        down2[,4] <- 1:nrow(down2)
        write.table(up2, "data/test_up", col.names = F,row.names = F, quote = F, sep = "\t")
        write.table(down2, "data/test_down", col.names = F,row.names = F, quote = F, sep = "\t")
        command <- c("intersectBed -wa -wb -a data/predictedTRUE_up_hg38 -b data/test_up|cut -f4,8 > data/PredictedAlupair")
        try(system(command))
        command <- c("intersectBed -wa -wb -a data/predictedTRUE_down_hg38 -b data/test_down|cut -f4,8 >> data/PredictedAlupair")
        try(system(command))
        command <- c("cat data/PredictedAlupair|sort |uniq -c|grep  '^ *2 ' >data/PredictInPair")
        try(system(command))
        if(file.info("data/PredictInPair")$size>0){
                pre_rslt <- read.table("data/PredictInPair", sep = "\t",header = F)
                pre_rslt[,1] <- gsub(".*\\s","",pre_rslt[,1])
                write.table(pre_rslt,"data/temp", quote = F,col.names = F,row.names = F,sep = "\t")             
                command <- c("sort -n -k1 data/temp >data/temp6 ")
                try(system(command))
                command <- c("awk 'NR==FNR {h[$1]; next} $4 in h' data/temp data/predictedTRUE_up_hg38 |sort -n -k4>data/temp2")
                try(system(command))            
                command <- c("awk 'NR==FNR {h[$1]; next} $4 in h' data/temp data/predictedTRUE_down_hg38 |sort -n -k4>data/temp3")
                try(system(command))             
                command <- c("join -1 1 -2 4 -a 1 data/temp6 data/temp2 >data/temp4")
                try(system(command))           
                command <- c('join -1 1 -2 4 -a 1 data/temp4 data/temp3 >data/temp5')
                try(system(command))          
                df <- read.table("data/temp5", header=F, sep = " " )
                df2 <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
                df2[,1]<- as.data.frame(df[,2])
                df2[,2] <- as.data.frame(paste(df[,3],":",df[,4]+1,"-",df[,5], " | ", df[,6],":",df[,7]+1,"-",df[,8],sep=""))    
                colnames(df2) <- c("No.","Predicted_Alu_pair")
                df3 <- data.frame(up2[,4],as.data.frame(paste(up2[,1],":",up2[,2]+1,"-",up2[,3], sep="")),
                                  as.data.frame(paste(down2[,1],":",down2[,2]+1,"-",down2[,3], sep="")) )
                colnames(df3) <-c("No.","5'_Genomic_Interval","3'_Genomic_Interval")
                df4 <- merge(df3,df2, all = T)
                command <- c("rm data/temp* data/test* data/Predict*")
                try(system(command))
                
                return(df4)
        }else{
                df3 <- data.frame(up2[,4],as.data.frame(paste(up2[,1],":",up2[,2]+1,"-",up2[,3], sep="")),
                                  as.data.frame(paste(down2[,1],":",down2[,2]+1,"-",down2[,3], sep="")), NA )
                colnames(df3) <-c("No.","5'_Genomic_Interval","3'_Genomic_Interval", "Predicted_Alu_pair")
                command <- c("rm data/test* data/Predict*")
                try(system(command))
                return(df3)
        }
 }