## ----LoadFunctions, echo=FALSE, message=FALSE, warning=FALSE, results='hide'----
library(knitr)
opts_chunk$set(error = FALSE)
library(EventPointer)
library(dplyr)
library(kableExtra)

## ----style, echo = FALSE, results = 'asis'------------------------------------
##BiocStyle::markdown()

## ----eval=FALSE---------------------------------------------------------------
# 
# library(BiocManager)
# 
# if (!requireNamespace("BiocManager", quietly=TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("EventPointer")

## -----------------------------------------------------------------------------
path_bam_files = system.file("extdata/bams", package = "EventPointer")
path_file_1 = dir(path_bam_files,pattern = "*.bam$",full.names = TRUE)[1]

file_tmp <- BamFile(path_file_1)
si <- seqinfo(file_tmp)
sl <- rep(seqlevels(si),2)
st <- c(rep(c("+"), length(si)),rep(c("-"), length(si)))
ranges_info <- GRanges(sl, IRanges(1, seqlengths(si)[sl]), st)

## ----echo=FALSE---------------------------------------------------------------
ranges_info = as.data.frame(ranges_info)

ranges_info$index = 1:nrow(ranges_info)
ranges_info = ranges_info[,c(6,1:5)]
levels(ranges_info$strand) = c("p","n","o")
# ranges_info
kable(ranges_info) %>%
  kable_paper() %>%
  scroll_box(width = "500px", height = "200px")

## ----eval=FALSE, warning=FALSE, collapse=TRUE---------------------------------
# 
# PathSamplesAbundance <- system.file("extdata/bams", package = "EventPointer")
# PathTranscriptomeGTF <- list.files(PathSamplesAbundance,"*.gtf",full.names = T)
# 
# PathSGResult <- tempdir()
# 
# EventsDetection_BAM(PathSamplesAbundance,
#                     PathTranscriptomeGTF,
#                     region = 16,
#                     min_junction_count = 2,
#                     max_complexity = 30,
#                     min_n_sample = NULL,
#                     min_anchor = 6,
#                     nboot = 20,
#                     lambda = NULL,
#                     cores = 1,
#                     PathSGResult = PathSGResult)
# 

## ----eval=FALSE---------------------------------------------------------------
# PathSamplesAbundance <- system.file("extdata/bams", package = "EventPointer")
# PathTranscriptomeGTF <- list.files(PathSamplesAbundance,"*.gtf",full.names = T)
# PathSGResult <- tempdir()
# 
# SgFC <- PrepareBam_EP(PathSamplesAbundance = PathSamplesAbundance,
#                         PathTranscriptomeGTF = PathTranscriptomeGTF,
#                         region = 16,
#                         min_junction_count = 2,
#                         max_complexity = 30,
#                         min_n_sample = NULL,
#                         min_anchor = 6,
#                         cores = 1,
#                         PathSGResult = PathSGResult,
#                         verbose = FALSE)

## ----eval=TRUE----------------------------------------------------------------
load(system.file("extdata/SgFC.RData", package = "EventPointer"))
SgFC

## ----eval=FALSE---------------------------------------------------------------
# EventsDetection_EPBAM<-EventDetection(SgFC, cores=1)

## ----echo=FALSE---------------------------------------------------------------
load(system.file("extdata/EventsDetection_EPBAM.RData", package = "EventPointer"))

## ----eval=FALSE---------------------------------------------------------------
# index_gene <- 1
# index_event <- 1
# length(EventsDetection_EPBAM) #number of genes
# length(EventsDetection_EPBAM[[index_gene]]) #number of events in the corresponding gene
# names(EventsDetection_EPBAM[[index_gene]][[index_event]]) #the different information stored in the object.

## ----eval=FALSE---------------------------------------------------------------
# event_list <- list()
# index <- 1
# for (geneEvent in EventsDetection_EPBAM) {
#   for (event in geneEvent) {
#     event_list[[index]] <- event$Info
#     index <- index + 1
#   }
# }
# totalEventTable <- do.call(rbind, event_list)
# write.csv(totalEventTable, file = paste0(PathSGResult, "/TotalEventsFound.csv"), row.names = FALSE)
# 

## ----eval=FALSE---------------------------------------------------------------
# cores <- 1
# lambda <- NULL
# nboot <- 20
# PSI_boots <- getPSI_RNASeq_boot(EventsDetection_EPBAM, lambda = lambda,
#                                 cores = cores, nboot = nboot)

## -----------------------------------------------------------------------------
data(PSI_boots)
Design <- cbind(rep(1, 9), rep(c(1, 0, 0), 3), rep(c(0, 1, 0), 3))
Contrast <- cbind(c(0, 1, 0), c(0, 0, 1))

PathSGResult <- tempdir()
result = EventPointerStats_BAM(PSI_boots,
                               Design,
                               Contrast,
                               Threshold = 0,
                               nbootstraps = 1000,
                               cores = 1,
                               ram = 0.1,
                               pathResult = PathSGResult)

## ----eval=FALSE, collapse=TRUE------------------------------------------------
# 
# EventsTxt<-paste(system.file("extdata",package="EventPointer"),"/TotalEventsFound.csv",sep="")
# load(paste(system.file("extdata",package="EventPointer"),"/SgFC.RData",sep=""))
# SG_RNASeq <- SgFC
# PathGTF <- tempdir()
# 
# EventPointerBAM_IGV(SG_RNASeq = SgFC,
#                     EventsCSV =  EventsTxt,
#                     PathGTF =  "./")

## ----GTFfGTF, eval=TRUE, warning=FALSE, collapse=TRUE-------------------------


# Set input variables
PathFiles<-system.file("extdata",package="EventPointer")
inputFile <- paste(PathFiles,"/gencode.v24.ann_2genes.gtf",sep="")
Pathtxt <- tempdir()

## ----GTFfGTF_2, eval=FALSE, warning=FALSE, collapse=TRUE----------------------
# # Run the function
# EventXtrans <- EventDetection_transcriptome(inputFile = inputFile,
#                                              Pathtxt=Pathtxt,
#                                              cores=1)
# 

## ----GTFfGTFnames, eval=FALSE, warning=FALSE, collapse=TRUE-------------------
# 
# names(EventXtrans)
# 

## ----getbootstrpdata, eval=TRUE, warning=FALSE, collapse=TRUE-----------------

PathSamples<-system.file("extdata",package="EventPointer")
PathSamples <- paste0(PathSamples,"/output")
PathSamples <- dir(PathSamples,full.names = TRUE)

data_exp <- getbootstrapdata(PathSamples = PathSamples,type = "kallisto")


## ----getdata_nb, eval=TRUE, warning=FALSE, collapse=TRUE----------------------

# this code chunk is an example of how to load the data from kallisto output.
# the expression of the isoforms are counts

PathFiles <- system.file("extdata",package="EventPointer")
filesnames <- dir(paste0(PathFiles,"/output"))
PathFiles <- dir(paste0(PathFiles,"/output"),full.names = TRUE)
dirtoload <- paste0(PathFiles,"/","abundance.tsv")
RNASeq <- read.delim(dirtoload[1],sep = "\t", colClasses = c(NA,"NULL","NULL",NA,"NULL"))
for (n in 2:length(dirtoload)){
  RNASeq[,n+1] <- read.delim(dirtoload[n],sep = '\t', colClasses = c('NULL','NULL','NULL',NA,'NULL'))
}
rownames(RNASeq)<-RNASeq[,1]
RNASeq<-RNASeq[,-1]
colnames(RNASeq) <- filesnames



## ----check_ann_nb, eval=TRUE, warning=FALSE, collapse=TRUE--------------------
rownames(RNASeq)[1:5]
EventXtrans$transcritnames[1:5]

## ----PSI_Statistic, eval=TRUE, warning=FALSE, collapse=TRUE-------------------

#change rownames of RNASeq variable
rownames(RNASeq) <- sapply(strsplit(rownames(RNASeq),"\\|"),function(X) return(X[1]))
RNASeq<-as.matrix(RNASeq) #must be a matrix variable


PSIss_nb <- GetPSI_FromTranRef(PathsxTranscript = EventXtrans,
                            Samples = RNASeq,
                            Bootstrap = FALSE,
                            Filter = FALSE)

PSI <- PSIss_nb$PSI
Expression <- PSIss_nb$ExpEvs


## ----check_ann, eval=TRUE, warning=FALSE, collapse=TRUE-----------------------
rownames(data_exp[[1]])[1:5]
EventXtrans$transcritnames[1:5]


## ----psi_with_bootstrap, eval=TRUE, warning=FALSE, collapse=TRUE--------------

#change rownames of the first element of teh list data_exp
rownames(data_exp[[1]]) <- sapply(strsplit(rownames(data_exp[[1]]),"\\|"),function(X) return(X[1]))

PSIss <- GetPSI_FromTranRef(PathsxTranscript = EventXtrans,
                               Samples = data_exp,
                               Bootstrap = TRUE,
                               Filter = FALSE)

PSI <- PSIss$PSI
Expression <- PSIss$ExpEvs

## ----ep_bootstrp_statistic, eval=TRUE, warning=FALSE, collapse=TRUE-----------

Dmatrix <- cbind(1,rep(c(0,1),each=2))
Cmatrix <- matrix(c(0,1),nrow=2)

Fit_bots <- EventPointer_Bootstraps(PSI = PSI,
                               Design = Dmatrix,
                               Contrast = Cmatrix,
                               cores = 1,
                               ram = 1,
                               nbootstraps = 10,
                               UsePseudoAligBootstrap = TRUE)


## ----ResultTable, eval=TRUE, warning=FALSE, collapse=TRUE---------------------

ResulTable(EP_Result = Fit_bots,coef = 1,number = 5)


## ----PSI_Statistic2, eval=TRUE, warning=FALSE, collapse=TRUE------------------

# Design and contrast matrix:

Design <- matrix(c(1,1,1,1,0,0,1,1),nrow=4)
Contrast <- matrix(c(0,1),nrow=2)

# Statistical analysis:

Fit <- EventPointer_RNASeq_TranRef(Count_Matrix = Expression,Statistic = "LogFC",Design = Design, Contrast = Contrast)


## ----EP_TranRef_Res_Table, echo=FALSE-----------------------------------------
kable(Fit,digits=5,row.names=FALSE,align="c",caption = "Table 3: PSI_Statistic results")

## ----ep_tranref_igv, eval=TRUE, warning=FALSE, collapse=TRUE------------------

SG_List <- EventXtrans$SG_List
PathEventsTxt<-system.file('extdata',package='EventPointer')
PathEventsTxt <- paste0(PathEventsTxt,"/EventsFound_Gencode24_2genes.txt")
PathGTF <- tempdir()

EventPointer_RNASeq_TranRef_IGV(SG_List = SG_List,pathtoeventstable = PathEventsTxt,PathGTF = PathGTF)


## ----biomart, eval=FALSE, warning=FALSE, collapse=TRUE------------------------
# 
# library(biomaRt)
# mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", host = "mar2016.archive.ensembl.org")
# mart<-useDataset("hsapiens_gene_ensembl",mart)
# 
# mistranscritos <- EventXtrans$transcritnames
# head(mistranscritos)
# 
# ## we need to remove the ".x":
# mistranscritos <- gsub("\\..*","",mistranscritos)
# 
# Dominios <- getBM(attributes = c("ensembl_transcript_id","interpro","interpro_description"),
#                   filters = "ensembl_transcript_id",
#                   values = mistranscritos,
#                   mart=mart)
# 
# #we build the isoform x protein domain matrix
# library(Matrix)
# ii <- match(Dominios$ensembl_transcript_id,mistranscritos)
# misDominios <- unique(Dominios[,2])
# jj <- match(Dominios[,2],misDominios)
# 
# TxD <- sparseMatrix(i=ii,j = jj,dims = c(length(mistranscritos),length(misDominios)),x = 1)
# rownames(TxD) <- mistranscritos
# colnames(TxD) <- misDominios
# 

## ----pd_enrichmetn, eval=FALSE, warning=FALSE, collapse=TRUE------------------
# 
# data("TxD")
# 
# #check same annotation for transcripts:
# EventXtrans$transcritnames[1]
# # [1] "ENST00000611540.4"
# rownames(TxD)[1]
# # [1] "ENST00000611540"
# 
# 
# ## as is not the same, we change EventXtrans$transcritnames annotation
# transcriptnames <- EventXtrans$transcritnames
# transcriptnames <- gsub("\\..*","",transcriptnames)
# EventXtrans$transcritnames <- transcriptnames
# 
# Result_PDEA <- Protein_Domain_Enrichment(PathsxTranscript = EventXtrans,
#                                          TxD = TxD,
#                                          Diff_PSI = Fit_bots$deltaPSI)
# 
# 

## ----pathprimer3, eval=FALSE, warning=FALSE, collapse=TRUE--------------------
# 
# Primer3Path <- Sys.which("primer3_core")
# 

## ----FindPrimers, eval=FALSE, warning=FALSE, collapse=TRUE--------------------
# 
# data("EventXtrans")
# #From the output of EventsGTFfromTranscriptomeGTF we take the splicing graph information
# SG_list <- EventXtrans$SG_List
# #SG_list contains the information of the splicing graphs for each gene
# #Let's supone we want to design primers for the event 1 of the gene ENSG00000254709.7
# 
# #We take the splicing graph information of the required gene
# SG <- SG_list$ENSG00000254709.7
# 
# #We point the event number
# EventNum <- 1
# 
# #Define rest of variables:
# Primer3Path <- Sys.which("primer3_core")
# Dir <- "C:\\PROGRA~2\\primer3\\"
# 
# 
# MyPrimers_taqman <- FindPrimers(SG = SG,
#                                 EventNum = EventNum,
#                                 Primer3Path = Primer3Path,
#                                 Dir = Dir,
#                                 mygenomesequence = BSgenome.Hsapiens.UCSC.hg38::Hsapiens,
#                                 taqman = 1,
#                                 nProbes=1,
#                                 nPrimerstwo=4,
#                                 ncommonForward=4,
#                                 ncommonReverse=4,
#                                 nExons=10,
#                                 nPrimers =5,
#                                 maxLength = 1200)
# 

## ----EP_DesiCP_Res_Table, eval=TRUE, warning=FALSE, collapse=TRUE,echo=FALSE----
data("MyPrimers")
kable(MyPrimers[1:5,],digits=5,row.names=FALSE,align="c",caption = "Table 4: Data.frame output of FindPrimers for conventional PCR") %>%
  kable_styling() %>%
  scroll_box(width ="660px")


## ----EP_DesiTP_Res_Table, eval=TRUE, warning=FALSE, collapse=TRUE,echo=FALSE----
data("MyPrimers_taqman")
kable(MyPrimers_taqman[1:5,],digits=5,row.names=FALSE,align="c",caption = "Table 5: Data.frame output of FindPrimers for conventional PCR") %>%
  kable_styling() %>%
  scroll_box(width ="660px")


## ----CDFGTF, eval=FALSE, warning=FALSE, collapse=TRUE-------------------------
# 
# # Set input variables
#    PathFiles<-system.file("extdata",package="EventPointer")
#    DONSON_GTF<-paste(PathFiles,"/DONSON.gtf",sep="")
#    PSRProbes<-paste(PathFiles,"/PSR_Probes.txt",sep="")
#    JunctionProbes<-paste(PathFiles,"/Junction_Probes.txt",sep="")
#    Directory<-tempdir()
#    array<-"HTA-2_0"
# 
# # Run the function
# 
#    CDFfromGTF(input="AffyGTF",inputFile=DONSON_GTF,
#               PSR=PSRProbes,Junc=JunctionProbes,
#               PathCDF=Directory,microarray=array)
# 

## ----aroma, eval=FALSE--------------------------------------------------------
# 
# # Simple example of Aroma.Affymetrix Preprocessing Pipeline
# 
# verbose <- Arguments$getVerbose(-8);
# timestampOn(verbose);
# projectName <- "Experiment"
# cdfGFile <- "EP_HTA-2_0,r"
# cdfG <- AffymetrixCdfFile$byChipType(cdfGFile)
# cs <- AffymetrixCelSet$byName(projectName, cdf=cdfG)
# bc <- NormExpBackgroundCorrection(cs, method="mle", tag=c("*","r11"));
# csBC <- process(bc,verbose=verbose,ram=20);
# qn <- QuantileNormalization(csBC, typesToUpdate="pm");
# csN <- process(qn,verbose=verbose,ram=20);
# plmEx <- ExonRmaPlm(csN, mergeGroups=FALSE)
# fit(plmEx, verbose=verbose)
# cesEx <- getChipEffectSet(plmEx)
# ExFit <- extractDataFrame(cesEx, addNames = TRUE)

## ----EP_arrays, eval=TRUE-----------------------------------------------------

   data(ArraysData)

   Dmatrix<-matrix(c(1,1,1,1,0,0,1,1),nrow=4,ncol=2,byrow=FALSE)
   Cmatrix<-t(t(c(0,1)))
   EventsFound<-paste(system.file("extdata",package="EventPointer"),"/EventsFound.txt",sep="")
   
   Events<-EventPointer(Design=Dmatrix,
                      Contrast=Cmatrix,
                      ExFit=ArraysData,
                      Eventstxt=EventsFound,
                      Filter=FALSE,
                      Qn=0.25,
                      Statistic="LogFC",
                      PSI=TRUE)

## ----EP_Arrays_Res_Table, echo=FALSE------------------------------------------
kable(Events[1:5,],digits=5,row.names=TRUE,align="c",caption = "Table 1: EventPointer Arrays results")

## ----Arrays_IGV, eval=FALSE, collapse=TRUE------------------------------------
# 
# # Set Input Variables
# 
#    DONSON_GTF<-paste(PathFiles,"/DONSON.gtf",sep="")
#    PSRProbes<-paste(PathFiles,"/PSR_Probes.txt",sep="")
#    JunctionProbes<-paste(PathFiles,"/Junction_Probes.txt",sep="")
#    Directory<-tempdir()
#    EventsFound<-paste(system.file("extdata",package="EventPointer"),"/EventsFound.txt",sep="")
#    array<-"HTA-2_0"
# 
# 
# # Generate Visualization files
# 
#  EventPointer_IGV(Events[1,,drop=FALSE],"AffyGTF",DONSON_GTF,PSRProbes,JunctionProbes,Directory,EventsFound,array)
# 
# 

## ----CDFGTF_MP, eval=FALSE, warning=FALSE, collapse=TRUE----------------------
# 
# # Set input variables
#    PathFiles<-system.file("extdata",package="EventPointer")
#    DONSON_GTF<-paste(PathFiles,"/DONSON.gtf",sep="")
#    PSRProbes<-paste(PathFiles,"/PSR_Probes.txt",sep="")
#    JunctionProbes<-paste(PathFiles,"/Junction_Probes.txt",sep="")
#    Directory<-tempdir()
#    array<-"HTA-2_0"
# 
# # Run the function
# 
#    CDFfromGTF_Multipath(input="AffyGTF",inputFile=DONSON_GTF,
#               PSR=PSRProbes,Junc=JunctionProbes,
#               PathCDF=Directory,microarray=array,paths=3)
# 

## ----EventDetection_MP, eval=FALSE--------------------------------------------
#   # Run EventDetection function
#    data(SG_RNASeq)
#    TxtPath<-tempdir()
#    AllEvents_RNASeq_MP<-EventDetectionMultipath(SG_RNASeq,cores=1,Path=TxtPath,paths=3)

## ----ListofLists_MP, eval=FALSE-----------------------------------------------
# Events[[i]][[j]]

## ----reclassification, eval=TRUE----------------------------------------------
#load splicing graph
data("SG_reclassify")

#load table with info of the events
PathFiles<-system.file("extdata",package="EventPointer")
inputFile <- paste(PathFiles,"/Events_found_class.txt",sep="")
EventTable <- read.delim(file=inputFile)
#this table has the information of 5 complex events.

EventTable_new <- Events_ReClassification(EventTable = EventTable,
                        SplicingGraph = SG_reclassify)

## ----PSI_ADV, eval=TRUE, collapse=TRUE----------------------------------------

# Microarrays (two paths)
data(ArraysData)
PSI_Arrays_list<-EventPointer:::getPSI(ArraysData)
PSI_Arrays <- PSI_Arrays_list$PSI
Residuals_Arrays <- PSI_Arrays_list$Residuals

# Microarrays (Multi-Path)
data(ArrayDatamultipath)
PSI_MP_Arrays_list <- EventPointer:::getPSImultipath(ArrayDatamultipath)
PSI_multipath_Arrays <- PSI_MP_Arrays_list$PSI
Residuals_multipath_Arrays <- PSI_MP_Arrays_list$Residuals

# RNASeq (two paths)
data(AllEvents_RNASeq)
PSI_RNASeq_list<-EventPointer:::getPSI_RNASeq(AllEvents_RNASeq)
PSI_RNASeq <- PSI_RNASeq_list$PSI
Residuals_RNASeq <- PSI_RNASeq_list$Residuals

# RNASeq (Multi-Path)
data(AllEvents_RNASeq_MP)
PSI_MP_RNASeq_list <- EventPointer:::getPSI_RNASeq_MultiPath(AllEvents_RNASeq_MP)
PSI_multipath_RNASeq <- PSI_MP_RNASeq_list$PSI
Residuals_multipath_RNASeq <- PSI_MP_RNASeq_list$Residuals


## ----PSI_ADV2, eval=TRUE, collapse=TRUE---------------------------------------

Dmatrix<-matrix(c(1,1,1,1,0,0,1,1),nrow=4,ncol=2,byrow=FALSE)
Cmatrix<-t(c(0,1))

table <- PSI_Statistic(PSI = PSI_Arrays,Design = Dmatrix,Contrast = Cmatrix,nboot = 20)


## ----PSI_ADV3, eval=TRUE, collapse=TRUE---------------------------------------

Dmatrix<-matrix(c(1,1,1,1,0,0,1,1),nrow=4,ncol=2,byrow=FALSE)
Cmatrix<-t(t(c(0,1)))

Ress <- vector("list", length = ncol(Cmatrix))

fitresiduals <- limma::lmFit(Residuals_Arrays,design = Dmatrix)
fitresiduals2 <- limma::contrasts.fit(fitresiduals, Cmatrix)
fitresiduals2 <- limma::eBayes(fitresiduals2)

# repeated if there is more than one contrast
for (jj in 1:ncol(Cmatrix)) {
  TopPSI <- limma::topTable(fitresiduals2, coef = jj, number = Inf)[, 1, drop = FALSE]
  colnames(TopPSI)<-"Residuals"
  Ress[[jj]] <- TopPSI
}



## -----------------------------------------------------------------------------
sessionInfo()

