#1  packages installation

install.packages('dplyr')
library(dplyr)
install.packages('tibble')
library(tibble)
install.packages('GenomicRanges')
library(GenomicRanges)
# install.packages("writexl")
# library("writexl")
install.packages('IDPmisc')
library(IDPmisc)
install.packages('ggplot2')
library(ggplot2)
install.packages('readbitmap')
library(readbitmap)


#2  it is possible to choose between test data or the entire dataset (CDS and UTR regions need to be further downloaded)

  #OPTION A
  #test data to run this program without any further files downloading; 
  #only chromosome 1 regions; CDS,UT5 and UT3 files do not need to be downloaded for this part because the BED files 
  #with chromosome 1 regions are already stored in this directory
  regions<-read.csv('test_regions.csv')
  UCSC_outputs<-c('test_cds.bed','test_utr5.bed','test_utr3.bed')

  #OPTION B
  #entire dataset of SV regions; CDS, UTR5 and UTR3 BED files are need to be downloaded from UCSC Table Browser for this analysis
  #upload deletion positions from 1000 GP
  regions <- read.csv('NIHMS64772-supplement-Table_S6.csv')
  #define CDS, UTR5 and UTR3 files names downloaded from UCSC Table Browser 
  UCSC_outputs<-c('place bed file containing CDS regions','place bed filecontaining UTR5 regions','place bed file containing UTR3 regions')


#3  GET OUTPUT 
#find region overlaps and get genes and RVIS score for each variant that deleted at least one exon
output <- findRegions(regions,UCSC_outputs)




##function findRegions: input- regions...deletion positions from 1000 GP
#                              GS...gene symbols for gene annotation downloaded from DAVID
#                     output- regions(data.frame)...1col=chromosome,2col=START,3col=END,4col=regions within particular deletion (UTR;CDS),
#                                                   5col=gene within particular deletion region, 6col=RVIS score for each site that deleted at least one exon
findRegions<-function(regions,UCSC_outputs){
             regions <- regions[, 1:3]
             df <- data.frame(
                   seqnames = regions[, 1],
                   start = regions[, 2],
                   end = regions[, 3],
                   stringsAsFactors = FALSE)
             GR <- makeGRangesFromDataFrame(df)
             regions_type<-c(UCSC_outputs)
             regions<-data.frame(append(regions, data.frame(region = "", gene = "",RVIS = NA,stringsAsFactors = FALSE), after =  5),stringsAsFactors = FALSE)
    
             for (i in 1:3){
                 reg<-c('CDS','UTR5','UTR3')
                 RT<-read.table(regions_type[i],header = FALSE, sep="\t",stringsAsFactors=FALSE, quote="")
                 NAMES<-RT$V4
                 NAMES<-gsub("\\..*","",NAMES)
                 RT<-RT[,1:3]
                 df_ex<-data.frame(seqnames = RT[,1],
                                   start = RT[,2],
                                   end = RT[,3],
                                   stringsAsFactors = FALSE)
                 GR_ex<-makeGRangesFromDataFrame(df_ex)
                 for (j in 1:length(GR)){
                     y<-c()
                     y<-findOverlaps(GR_ex, GR[j,], type='within') #function finding CDS, UTR5 or UTR3 regions within SV regions.
                   
                     if(length(y)>0){
                       gene_id<-c()
                       IDs<-c()
                       gene_id<-NAMES[y@from] #getting gene names of overlapping regions found by function findOverlaps
                       gene_id<-unique(unlist(strsplit(gene_id, " "))[])
                       gene_id<-toupper(RefSeqToSymbol(gene_id)) #RefSeq gene annotations are transformed into gene symbols
                       IDs<-c(regions[j,5],gene_id)
                       IDs<-unique(unlist(strsplit(IDs, ","))[])
                       IDs<-paste(IDs,collapse = ',')
                       regions[j,5]<-IDs
                       pom<-c(regions[j,4],reg[i])
                       pom<-pom[pom != ""] #eliminate an empty string in vector 
                       pom<-paste(pom,collapse = ',')
                       regions[j,4]<-pom
                       RVIS_score<-RVIS(gene_id) #assigning minimum RVIS score for a particular deletion site
                       if (is.na(RVIS_score) & is.na(regions[j,6])){ 
                          regions[j,6]<-NA #to avoid INF value in min() function
                          }
                       else {
                       regions[j,6]<-min(regions[j,6],RVIS_score,na.rm = T) #minimum RVIS score per site
                       }
                       } 
                    }
             }
             regions <- do.call(data.frame, lapply(regions, function(x) replace(x, is.infinite(x), NA)))
             return(regions)
}




#function RefSeqToSymbol: converts RefSeq gene annotation into symbol:
#                      input- vector(character)...vector containing RefSeq gene annotations
#                     output- symbols(character)...vector containing gene symbols
RefSeqToSymbol<-function(vector){
                GS<-read.table('gene_SYMBOL.bed', header = T, sep = "\t", stringsAsFactors = FALSE, quote = "")
                mat<-GS %>%
                mutate_all(list(~ . %in% vector)) %>% #find gene positions in GS data frame 
                as.matrix %>% 
                which(., arr.ind = TRUE)
                symbols<-GS[mat[,1], 'To'] #get symbols for defined positions in GS data frame 
                return(symbols)
}

#function RVIS: assign minimum RVIS score per deletion site with annotated genes 
#               input- vector(character)...vector containing gene IDs in gene symbol format
#              output- RVIS (integer)...minimum RVIS score of inputted genes 
RVIS<-function(vector){
      RVIS_score<-read.table('RVIS.bed',header = T, sep="\t",stringsAsFactors=FALSE, quote="") #upload table containing RVIS scores
      mat<-RVIS_score %>%  #find positions containing inputted gene names in RVIS data frame 
      mutate_all(list(~ . %in% vector)) %>%
          as.matrix %>% 
          which(., arr.ind = TRUE)
      pom<-round(RVIS_score[mat[,1], 'X.ALL_0.1.']) #get RVIS score for each gene via its position in RVIS data frame 
      pom<-na.omit(pom) 
      if(isEmpty(pom)==T){
        RVIS<-NA  #to avoid an INF value produced by min() function
      } else {
      RVIS<-min(pom,na.rm = T)
      }
      return(RVIS)
}



