inst/doc/EventPointer.R

In EventPointer: An effective identification of alternative splicing events using junction arrays and RNA-Seq data

## ----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=TRUE, 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=TRUE, 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()
PathGTF <- tempdir()

# Run the function
EventXtrans <- EventsGTFfromTrancriptomeGTF(inputFile = inputFile,
                                           Transcriptome = Transcriptome,
                                           Pathtxt=Pathtxt,
                                           PathGTF=PathGTF)


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

names(EventXtrans)


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

#first: load data from kallisto output

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","NULL",NA))
for (n in 2:length(dirtoload)){
  RNASeq[,n+1] <- read.delim(dirtoload[n],sep = '\t', colClasses = c('NULL','NULL','NULL','NULL',NA))
}
rownames(RNASeq)<-RNASeq[,1]
RNASeq<-RNASeq[,-1]
colnames(RNASeq) <- filesnames


#second: the annotation of the transcript names must be equal in RNASeq and PathsxTranscript

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

#third: Obtain values of PSI
PSIss <- GetPSI_FromTranRef(PathsxTranscript = EventXtrans,Samples = RNASeq,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")

## ----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,
#                          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=TRUE, 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=TRUE---------------------------------------------
  # 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]]

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