## ----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")

## ----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)
#  
#  

## ----PrepareBam, eval=FALSE, collapse=TRUE------------------------------------
#  # Obtain the samples and directory for .bam files
#  
#  # the object si contains example sample information from the SGSeq R package
#  # use ?si to see the corresponding documentation
#  
#     BamInfo<-si
#     Samples<-BamInfo[,2]
#     PathToSamples <- system.file("extdata/bams", package = "SGSeq")
#     PathToGTF<-paste(system.file("extdata",package="EventPointer"),"/FBXO31.gtf",sep="")
#  
#    # Run PrepareBam function
#     SG_RNASeq<-PrepareBam_EP(Samples=Samples,
#                              SamplePath=PathToSamples,
#                              Ref_Transc="GTF",
#                              fileTransc=PathToGTF,
#                              cores=1)

## ----EventDetection, eval=TRUE------------------------------------------------
  # Run EventDetection function
   data(SG_RNASeq)
   TxtPath<-tempdir()
   AllEvents_RNASeq<-EventDetection(SG_RNASeq,cores=1,Path=TxtPath)
   

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

## ----EP_RNASeq, eval=TRUE-----------------------------------------------------
   Dmatrix<-matrix(c(1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1),ncol=2,byrow=FALSE)
   Cmatrix<-t(t(c(0,1)))
   Events <- EventPointer_RNASeq(AllEvents_RNASeq,Dmatrix,Cmatrix,Statistic="LogFC",PSI=TRUE)

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

## ----RNAS_IGV, eval=TRUE, collapse=TRUE---------------------------------------

   # IGV Visualization
   EventsTxt<-paste(system.file("extdata",package="EventPointer"),"/EventsFound_RNASeq.txt",sep="")
   PathGTF<-tempdir()
   EventPointer_RNASeq_IGV(Events,SG_RNASeq,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="")
Transcriptome <- "Gencode24_2genes"
Pathtxt <- tempdir()

## ----GTFfGTF_2, eval=FALSE, warning=FALSE, collapse=TRUE----------------------
#  # Run the function
#  EventXtrans <- EventDetection_transcriptome(inputFile = inputFile,
#                                               Transcriptome = Transcriptome,
#                                               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

## ----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_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 <- 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,coef = 1,number = 5)


## ----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=TRUE, warning=FALSE, collapse=TRUE-------------------

data("TxD")

#check same annotation for transcripts:
EventXtrans$transcritnames[1]
rownames(TxD)[1]


## as si 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$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_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()

