aroma.affymetrix pre-processing pipeline

For microarray data, a pre-processing pipeline must be applied to the .cel files of the experiment. Usually this pre-processing is done using the aroma.affymetrix R package. This procedure normalizes and summarizes the expression of the different probesets into single values.

The aroma.affymetrix R package provides users different functions to work with affymetrix arrays. The functions are used to extract all the information contained in the .cel files and to do all the required pre-processing steps such as background correction, normalization and summarization. The package requires a certain folder structure in order to work correctly. For more information about aroma.affymetrix visit the webpage:aroma project

The following code chunk displays the pipeline used to get the results required by EventPointer after the pre-processing using aroma.affymetrix. The code is not intended to be a runable example, but just to show users the settings and functions that can be used. In order for users to have a runable example, the corrrespoding folder structure and files are required.

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

EventPointer function

EventPointer is the main function of the algorithm

The function requires the following parameters:

• Design : The design matrix for the experiment.
• Contrast: The contrast matrix for the experiment.
• ExFit: [aroma.affymetrix] pre-processed variable after using extractDataFrame(affy, addNames=TRUE)
• Eventstxt: Path to the EventsFound.txt file generated by CDFfromGTF function.
• Filter: Boolean variable to indicate if an expression filter is applied.
• Qn: Quantile used to filter the events (Bounded between 0-1, Q1 would be 0.25).
• Statistic: Statistical test to identify differential splicing events, must be one of : “LogFC”,“Dif_LogFC” and “DRS”.
• PSI: Boolean variable to indicate if ΔΨ should be calculated for every splicing event.

If the Filter variable is TRUE, the following is performed:

The summarized expression value of all the reference paths is obtained and the threshold is set depending on the Qn value used.

An event is considered if at least one sample in all paths is expressed above the threshold.

The rest of the events are not shown unless the Filter variable is set to FALSE

   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)

## 07:18:55 AM Running EventPointer:  
##  ** Statistical Analysis: Logarithm of the fold change of both isoforms 
##  ** Delta PSI will be calculated 
##  ** No expression filter 
##             ---------------------------------------------------------------- 
##  ** Calculating PSI...Done 
##  ** Running Statistical Analysis...Done 
## 
##  07:18:55 AM  Analysis Completed

Table 1 displays the output of EventPointer function

The columns of the data.frame are:

• Gene name : Gene identifier
• Event Type: Type of alternative splicing event (Cassette, Alternative 3’,Alternative 5’, Alternative Last Exon, Alternative First Exon, Mutually Exclusive Exons or Complex Event)
• Genomic Position: Genomic Coordinates for the event
• Splicing Z Value: Corresponding Z value for the statistical test performed
• Splicing P Value: Corresponding P-value for the statistical test performed
• Delta PSI: ΔΨ parameter for each event

Overview of RNA-Seq

EventPointer is also able to identify alternative splicing events from RNA-Seq data. The only required files are the corresponding .BAM files from the experiment.

Each time an experiment is analyzed with EventPointer, the complete process needs to be executed which may cause long waiting times to get the results. To avoid this issue, we have divided every step of the algorithm in different functions so as the user can reuse previous result and thus reduce computational time.

For the examples in this part of the vignette, we will use .bam files depicted in the SGSeq vignette that correspond to paired-end RNA-seq data from four tumor and four normal colorectal samples, which are part of a data set published in Seshagiri et al. 2012. As stated in SGSeq vignette the bam files only include reads mapping to a single gene of interest (FBXO31).

Note: For sequencing data, there are two requirements for the BAM files in order to get EventPointer working correctly:

1. The BAM files should include the XS-flag, the flag can be included in the files when running the alignment. Most of the available software has parameters to incude the flag. For example, in the case of STAR the flag –outSAMattributes XS must be included when mapping

2. All files to be analyzed should have the corresponding index files (.bai) in the same directory as the BAM files. Create the index before running EventPointer.

BAM Preparation

The first step to analyze alternative splicing events in RNA-Seq data, is the creation of the splicing graphs. This step relies entirely on SGSeq R package.

The function PrepareBam_EP transforms all the information found in the bam files into splicing graph features and counts

• Samples : Name of the .bam files to be analyzed (Sample1.bam,Sample2.bam,…,etc)
• SamplePath: Path where the bam files are stored
• Ref_Transc: Reference transcriptome used to name the genes found. Options are: “Ensembl”,“UCSC” or “GTF”.
• fileTransc: Path to the GTF reference transcriptome if Ref_Transc is “GTF”.
• cores: Number of cores used for parallel processing
• Alpha: Internal SGSeq parameter to include or exclude regions (See Advanced Use)

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

The output of PrepareBam_EP function is a SGFeaturesCounts object, for more information check SGSeq Vignette. Briefly the SGFeaturesCounts contains a GRanges object with all the required elements to create the different splicing graphs found in the given samples. It also contains the number of counts associated with each element of the splicing graph.

Event Detection

After running PrepareBam_EP, we have all the elements to analyze each of the splicing graphs. The next step is to identify and classify all the events, that are present in the BAM files.

For this purpose, the function EventDetection is used.

• Input : Output of the PrepareBam_EP function
• cores: Number of cores used for parallel processing
• Path: Directory where to write the EventsFound.txt file

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

## 
##  Obtaining Events  |                                                                              |                                                                      |   0%  |                                                                              |============                                                          |  17%

## 'as(<dgCMatrix>, "dgTMatrix")' is deprecated.
## Use 'as(., "TsparseMatrix")' instead.
## See help("Deprecated") and help("Matrix-deprecated").

##   |                                                                              |=======================                                               |  33%  |                                                                              |===================================                                   |  50%  |                                                                              |===============================================                       |  67%  |                                                                              |==========================================================            |  83%  |                                                                              |======================================================================| 100%

This function retireves a list with all the events found for all the genes present in the experiment. It also generates a file called EventsFound_RNASeq.txt with the information for every detected event.

The list is organized in the following way:

Events[[i]][[j]]

The list will have as many i values as genes and j values as many events detected for the i_th gene. In other words, the command above will display the j_th event detected for the i_th gene.

Statistical Analysis

The statistical analysis of the alternative splicing events is done in exactly the same way as for junction arrays. With the Design and Contrast matrices, the algorithm gives the statistical significance and ΔΨ.

The function for the statistical analysis using EventPointer method.

• Events : Output from EventDetection function
• Design: The design matrix for the experiment.
• Contrast: The contrast matrix for the experiment.
• Statistic: Statistical test to identify differential splicing events, must be one of : “LogFC”,“Dif_LogFC” and “DRS”.
• PSI: Boolean variable to indicate if ΔΨ should be calculated for every splicing event.

The algorithm displays the different parameters that are selected to perform the analysis.

Following our example, the code chunk to obtain the results:

   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)

## 07:18:56 AM Running EventPointer:  
##  ** Statistical Analysis: Logarithm of the fold change of both isoforms 
##  ** Delta PSI will be calculated 
##             ---------------------------------------------------------------- 
##  ** Calculating PSI...Done
##  Analysis Finished
##  Done 
## 
##  07:18:57 AM  Analysis Completed

Table 2 displays the output of EventPointer function

IGV visualization

EventPointer creates one GTF file that can be loaded into IGV to visualize the alternative splicing events. Figure 6 displays an example result showed in IGV (5th ranked event in Table 2). Also, in the figure a reference transcriptome is displayed (blue track), and it can be seen that the displayed event corresponds to a novel event discovered with sequencing data and that it will not be detected using junction arrays.

Figure 6. EventPointer visualization using IGV

To create the GTF files, the algorithm uses the following code.

• Events : Data.frame generated by EventPointer_RNASeq with the events to be included in the GTF file.
• SG_RNASeq: Output from PrepareBam_EP function. Contains splicing graphs components.
• EventsTxt: Path to EventsFound.txt file generated with EventDetection function
• PathGTF: Directory where to write the GTF files.

A progress bar is displayed to the user to inform about the progress of the function.

Once the process is completed the GTF file is generated in the specified directory:

• paths_RNASeq.gtf : GTF file representing the alternative splicing events.

   # IGV Visualization
   EventsTxt<-paste(system.file("extdata",package="EventPointer"),"/EventsFound_RNASeq.txt",sep="")
   PathGTF<-tempdir()
   EventPointer_RNASeq_IGV(Events,SG_RNASeq,EventsTxt,PathGTF)
## 
##  Generating GTF Files...  |                                                                              |                                                                      |   0%  |                                                                              |====                                                                  |   6%  |                                                                              |========                                                              |  12%  |                                                                              |============                                                          |  18%  |                                                                              |================                                                      |  24%  |                                                                              |=====================                                                 |  29%  |                                                                              |=========================                                             |  35%  |                                                                              |=============================                                         |  41%  |                                                                              |=================================                                     |  47%  |                                                                              |=====================================                                 |  53%  |                                                                              |=========================================                             |  59%  |                                                                              |=============================================                         |  65%  |                                                                              |=================================================                     |  71%  |                                                                              |======================================================                |  76%  |                                                                              |==========================================================            |  82%  |                                                                              |==============================================================        |  88%  |                                                                              |==================================================================    |  94%  |                                                                              |======================================================================| 100%

Event Detection

We use EventDetection_transcriptome to identify and classify alternative splicing events of a given reference transcriptome. The required parameters are:

• inputFile = inputFile should point to the GTF file to be used.
• Transcriptome = the name of the transcriptome.
• Pathtxt = Directory to save the .txt file of the events found.
• cores = Number of cores using in the parallel processing (by default = 1).

Following code shows how to apply EventDetection_transcriptome function using the Gencode 24 transcriptome (GRCH 38) as reference annotation. In this example we do not use the full reference but only two genes (ENSG00000185252.17 and ENSG00000254709.7).

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

# Run the function
inputFile <- "/home/external_HDD/HDD_1/EventPointer/inst/extdata/gencode.v24.ann_2genes_prueba.gtf"
EventXtrans <- EventDetection_transcriptome(inputFile = inputFile,
                                             Transcriptome = Transcriptome,
                                             Pathtxt=Pathtxt,
                                             cores=1)

This function takes a couple of hours to complete (depending on the computer), and creates a file called EventsFound_.txt: a .txt file with the information of the events. This .txt file will be used as an input in further functions.

This function also returns a list containing five elements: three sparce matrices that relate which isoforms build up the paths (path1,path2 and pathRef) of each event. The fourth element contains the name of the reference annotation: only appear the name of the transcript. The final element is SG_List: a list with the information of the graph of each gene, this variable is necesary for Primers design step.

names(EventXtrans)

Get expression and PSI (Ψ)

Given (I) the output of EventDetection_transcriptome, which contains events information and (II) the result of pseudo-alignment, EventPointer is able to get the expression of the path of each event and also to compute the PSI value.

To first step is to load the data from pseudo-alignment. EventPointer works with both Kallisto and Salmon methods.

In this example we have applied kallisto to the samples ERR315326, ERR315400, ERR315407 and ERR315494 using the same reference used before. These samples have been download from EMBL-EBI dataset.

EventPointer offers the option to take advantage of bootstraps that return both Kallisto and Salmon. In case we want to leverage the Bootstrap data we have to use the function getbootstrapdata:

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

data_exp <- getbootstrapdata(PathSamples = PathSamples,type = "kallisto")
## 1 2 3 4

getbootstrapdata returns a list containing the quantification data with the bootstrap information. The list is of length equal to the number of samples. Each element of the list storage a matrix where the number of rows is equal to the number of transcripts and the number of columns is equal to the number of bootstrats + 1. (the first column corresponds to the maximum likelihood expression and the rest to the bootstrap data).

In case you don’t want to use the bootstrap returned by the pseudo-alignmnet step (or is no available), you should load a matrix with the expression of the transcripts. The following code chunk shown an example of how to do this:

# 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

Once we have the expression of the isoforms loaded, the next step is to compute the PSI value. For both cases (with or withuot bootstrap data from psudoaligment) EventPointer provides the function GetPSI_FromTranRef for this purpose.

The function GetPSI_FromTranRef returns the values of Ψ and the expression of the paths of each events. This functions requires the following inputs:

• PathsxTranscript The output of EventDetection_transcriptome.
• Samples A matrix or list containing the expression of the samples.
• Bootstrap Boolean variable to indicate if bootstrap data from pseudo-alignment is used.
• Filter Boolean variable to indicate if an expression filter is applied. Defaul TRUE.
• Qn Quantile used to filter the events (Bounded between 0-1, Q1 would be 0.25).

This function also requires that the same reference transcriptome have been used in both the pseudo-alignment and event detections steps. Besides, in the variable Samples the type of annotation used must be equal to the one used in the variable EventXtrans.

Get psi if no bootstrap from pseudo-alignmnet is used

First step is to verify if same annotation is used in RNASeq and EventXtrans:

rownames(RNASeq)[1:5]
## [1] "ENST00000400451.6|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319648.2|ZNF74-005|ZNF74|3159|protein_coding|" 
## [2] "ENST00000611540.4|ENSG00000185252.17|OTTHUMG00000150687.4|-|ZNF74-202|ZNF74|3714|protein_coding|"                    
## [3] "ENST00000403682.7|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319645.2|ZNF74-002|ZNF74|2789|protein_coding|" 
## [4] "ENST00000357502.5|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319622.3|ZNF74-001|ZNF74|2788|protein_coding|" 
## [5] "ENST00000493734.5|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319646.2|ZNF74-003|ZNF74|3758|retained_intron|"
EventXtrans$transcritnames[1:5]
## [1] "ENST00000611540.4" "ENST00000400451.6" "ENST00000403682.7"
## [4] "ENST00000357502.5" "ENST00000493734.5"

We need to change the rownames of RNASeq variable before applying GetPSI_FromTranRef function:

#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

The output of the function is a list containing two elements: a matrix with the Ψ values, and a list containing as many matrices as number of events. In each matrix is stored the expression of the different paths of an event along the samples.

If bootstrap data from pseudo-alignment is required

As done above, we need to check if same annotation is used in data_exp and EventXtrans. The difference from above is that data_exp is a list and not a matrix as RNASeq. Therefore, we are going to used the rownames of the first list of data_exp to check if same annotations is used:

rownames(data_exp[[1]])[1:5]
## [1] "ENST00000400451.6|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319648.2|ZNF74-005|ZNF74|3159|protein_coding|" 
## [2] "ENST00000611540.4|ENSG00000185252.17|OTTHUMG00000150687.4|-|ZNF74-202|ZNF74|3714|protein_coding|"                    
## [3] "ENST00000403682.7|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319645.2|ZNF74-002|ZNF74|2789|protein_coding|" 
## [4] "ENST00000357502.5|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319622.3|ZNF74-001|ZNF74|2788|protein_coding|" 
## [5] "ENST00000493734.5|ENSG00000185252.17|OTTHUMG00000150687.4|OTTHUMT00000319646.2|ZNF74-003|ZNF74|3758|retained_intron|"
EventXtrans$transcritnames[1:5]
## [1] "ENST00000611540.4" "ENST00000400451.6" "ENST00000403682.7"
## [4] "ENST00000357502.5" "ENST00000493734.5"

We need to change the rownames of the first element of the list data_exp before applying GetPSI_FromTranRef function:

#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

The output of the function is a list containing two elements: an array with the Ψ values with dimension = c(number of events, number of bootstrap+1,number of samples), and a list containing as many matrices as number of events. In each matrix is stored the expression of the different paths of an event along the samples.

Statistical Analysis

The statistical analysis of the alternative splicing events is done in exactly the same way as for junction arrays.

The function EventPointer_RNASeq_TranRef perform this statistical analysis and requires the following parameters:

• Count_Matrix The list containing the expression data taken from the ouput of GetPSI_FromTranRef
• Statistic The type of statistic to apply. Default = “LogFC” (can be “logFC,”Dif_LogFC”,“DRS”)
• Design The design matrix of the experiment.
• Contrast The Contrast matrix of the experiment.

# 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)
## Done
##  Analysis Finished
##  Done 
## 
##  07:18:59 AM  Analysis Completed

The output of this function is a data.frame with the information of the names of the event, its p.values and the corresponding z.value. If there is more than one contrast, the function returns as many data.frames as number of contrast and all these data.frame are sotred in an unique list.

Statistical Analysis (Bootstrap Test)

We have implemented a bootstrap test to evaluate the difference in Ψ among the conditions under study. This test can be done both for PSI values that contain information from the pseudo-alignment bootstrap and for cases in which we do not have this information.

The function EventPointer_Bootstraps proceed to calculate the increase in PSI and its p.values by performing bootstrap statistic test. This function can be parallelized. Users may set the number of Cores (by default 1) and the ram available (by default 1 Gb).

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)
## 
##  The Contrast and Design Matrices have been correctly designed. 
## 
## Number of chunks:  1 
## 
## Performing resampling of chunk number 1 , this may take a while...
## 
## ****Analysis 100 % Completed****
## 
## Preparing output...
##  finish chunk  1 
## 
##  The program has succesfully ended.

You can extract a table of the top-ranked events from the results with the function ResulTable

ResulTable(EP_Result = Fit,coef = 1,number = 5)
##                            deltaPSI      pvalue      lfdr   qvalues
## ENSG00000185252.17_5  -0.0001254954 0.009787789 0.4267450 0.1174535
## ENSG00000185252.17_22  0.0688425350 0.098111380 0.9357568 0.5886683
## ENSG00000254709.7_3    0.0453065788 0.312944752 1.0000000 0.9665380
## ENSG00000185252.17_23  0.0688425350 0.322179318 1.0000000 0.9665380
## ENSG00000185252.17_21 -0.0000418445 0.737901068 1.0000000 1.0000000

IGV visualization

EventPointer creates one GTF file that can be loaded into IGV to visualize the alternative splicing events. Figure XX displays an example result showed in IGV (5th ranked event in Table 2). Also, in the figure a reference transcriptome is displayed (blue track).

To create the GTF files, the algorithm uses the EventPointer_RNASeq_TranRef_IGV functions with the following parameters:

• SG_List: List with the Splicing Graph information of the events. This list is created by EventDetection_transcriptome function.
• pathtoeventstable: Complete path to the table returned by EventDetection_transcriptome that contains the information of each event, or table with specific events that the user want to load into IGV to visualize.
• PathGTF: Directory where to write the GTF files.

As EventDetection_transcriptome function returns the splicing graph information, this function does not need to create the splicing graph. (EventPointer_IGV is equivalent function for array platform).

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)
## 
## Creating .GTF ...  |                                                                              |                                                                      |   0%  |                                                                              |======                                                                |   8%  |                                                                              |============                                                          |  17%  |                                                                              |==================                                                    |  25%  |                                                                              |=======================                                               |  33%  |                                                                              |=============================                                         |  42%  |                                                                              |===================================                                   |  50%  |                                                                              |=========================================                             |  58%  |                                                                              |===============================================                       |  67%  |                                                                              |====================================================                  |  75%  |                                                                              |==========================================================            |  83%  |                                                                              |================================================================      |  92%  |                                                                              |======================================================================| 100%
## 
## .GTF created

Domain Enrichment

The main idea of protein domain analysis is to analyze whether the presence of a protein domain increases or decreases in the condition under study. To do this, we have implemented the function Protein_Domain_Enrichment. This function requires the following parameters:

• PathsxTranscript: the output of EventDetection_transcriptome.
• TxD: matrix that relates transcripts with Protein domain. Users can get it from BioMart.
• Diff_PSI: matrix with the difference of psi of the condition under study.
• method: a character string indicating which correlation coefficient is to be calculated. “spearman” (default) or “pearson” can be selected.

As mentioned, it is first necessary to know the structure of protein domains within the transcriptome or in other words: to relate the isoforms to protein domains. This can be done using ‘bioMart’ R package. The input needed by EP to perform this analysis is a matrix that relates transcripts to protein domains (TxD matrix). Therefore, you can work either with Pfam, Interpro or superfamily annotation.

The following conde chunck shows an example of how to get this TxD matrix:

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

The above code chunk build the required Transcript X Domain (TxD) sparse matrix where TxD_ij = 1 means that Domain j match to transcript i.

The next step is to perform the enrichment analysis with the function Protein_Domain_Enrichment Internally, given the TxD matrix Protein_Domain_Enrichment relates protein domains to events and performs the enrichment analysis:

NOTE: Check that EventXtrans$transcritnames annotation type match to the rownames of the TxD matrix.

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$deltaPSI)

The function returns a list containing the results of the protein domain enrichment analysis. This list contains 3 matrices in which the rows indicate the protein domains and the columns the number of contrasts. The 3 matrices are the following:

• mycor: correlation value between the deltaPSI and the DifProtDomain matrix (see more details in vignette)
• STATISTIC: the values of the test statistic
• PVAL: the pvalues of the test statistic

Primers Design

EventPointer can also be used to design primers for PCR validation. The aim of FindPrimers function is the design of PCR primers and TaqMan probes for detection and quantification of alternative splicing. Depending on the assay we want to carry out the the algorithm will design the primers for a conventional PCR or the primers and TaqMan probes if we are performing a TaqMan assay. In the case of a conventional PCR we will be able to detect the alternative splicing event. Besides, the algorithm gives as an output the length of the PCR bands that are going to appear. In the case of a TaqMan assay, we will not only detect but also quantify alternative splicing.

The Primers Design step has been developed with Primer3 software and works with versions >= 2.3.6. In order to use this option you need to install this sofware. You can download from sourceforge (for Windows, Mac OSX or Unix/Linux) or from github (for Mac OSX or Unix/Linux).

To work with Primer Design option you also need to add Primer3 to your environment PATH. For Windows you can add environment variables in:

Control Panel -> System and Security -> System -> Advanced system settings -> Environment Variables…

For Unix/Linux or Max you need to add the path to Primer3 to the variable PATH using terminal.

Then, as is explained below, one of the input variables for the main function is Primer3Path which is a string variable with the complete path where primer3_core.exe is placed. As we have added this path to the environment variable we can set down this variable with the following command as is shown later in the example:

Primer3Path <- Sys.which("primer3_core")

The main function of this step is called Find_Primers and requires the following main parameters:

• SG: Information of the graph of the gene where the selected event belongs.This information is available in the output of EventsGTFfromTranscriptomeGTF function.
• EventNum: The “EventNum” variable can be found in the returned .txt file from the EventsGTFfromTranscriptomeGTF function in the column “EventNumber” or in the output of EventPointer_RNASeq_TranRef, the number after the “_” character of the ‘Event_ID’.
• Primer3Path: Complete path where primer3_core.exe is placed.
• Dir: Complete path where primer3web_v4_0_0_default_settings.txt file and primer3_config directory are placed.
• taqman: 1 if you want to get probes and primers for TaqMan. 0 if you want to get primers for conventional PCR.

Other required parameters that have a default value but can be modified by user:

• nProbes: Number of probes for TaqMan experiments. By default 1.
• nPrimerstwo: Number of potential exon locations for primers using two primers (one forward and one reverse). By default 3.
• ncommonForward: Number of potential exon locations for primers using one primer in forward and two in reverse. By default 3.
• ncommonReverse: Number of potential exon locations for primers using two primer in forward and one in reverse. By default 3.
• nExons: Number of combinations of ways to place primers in exons to interrogate an event after sorting. By default 5.
• nPrimers: Once the exons are selected, number of primers combination sequences to search within the whole set of potential sequences.By default 5.
• shortdistpenalty: Penalty for short exons following an exponential function. By default 2000.
• maxLength: Max length of exons that are between primers and for paths once we have calculated the sequence.By default 1000.
• minsep: Distance from which it is penalized primers for being too close. By default 100.
• wminsep: Weigh of the penalization to primers for being too close. By default 200.
• valuethreePenalty: penalization for cases that need three primers instead of 2. By default 1000.
• minexonlength: Minimum length that a exon has to have to be able to contain a primer. By default 25.
• wnpaths: Penalty for each existing path. By default 200.
• qualityfilter: Results will show as maximum 3 combinations with a punctuation higher than qualityfilter.By default 5000.

The following example shows how to design primers and TaqMan probes for a specific alternative splicing event:

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)

FindPrimers return a data.frame. Either for conventional PCR and taqman option the data.frame has the following columns:

• For1Seq: Sequence of the first forward primer.
• For2Seq: Sequence of the second forward primer in case it is needed.
• Rev1Seq: Sequence of the first reverse primer.
• Rev2Seq: Sequence of the second reverse primer in case it is needed.
• For1Exon: Name of the exon of the first forward primer.
• For2Exon: Name of the exon of the second forward primer in case it is needed.
• Rev1Exon: Name of the exon of the first reverse primer.
• Rev2Exon: Name of the exon of the second reverse primer in case it is needed.
• FINALvalue: Final punctuation for that combination of exons and sequences. The lower it is this score, the better it is the combination.
• DistPath1: Distances of the bands, in base pairs, that interrogate Path1 when we perform the conventional PCR experiment.
• DistPath2: Distances of the bands, in base pairs, that interrogate Path2 `when we perform the conventional PCR experiment.
• DistNoPath: Distances of the bands, in base pairs, that they do not interrogate any of the two paths when we perform the conventional PCR experiment.

For the tacman option, the data.frame contains three additional columns:

• SeqProbeRef: Sequence of the TaqMan probe placed in the Reference.
• SeqProbeP1: Sequence of the TaqMan probe placed in the Path1.
• SeqProbeP2: Sequence of the TaqMan probe placed in the Path2.

IGV primers and probes visualization

We can visualize the result in IGV using the Find Motif tool given the output of FindPrimers function.

Forward primers are represented in blue whilst reverse primers are represented in red as long as the gene is transcribed from the forward strand. Otherwise the colors we will be swapped. TaqMan probes will always have the color of the forward primers.

Figure 8. Transcriptome, paths and primers visualization using IGV.

RBP prediction (Future Lines)

Various processes (probably many of them yet to be discovered) are involved in the regulation of alternative splicing. One of this regulators are the RNA-binding proteins (RBPs).
we have integrated several databases of CLIP-seq experiments with an algorithm that detects differential splicing events to discover RBPs that are especially enriched. Thus, EP proposes a ranking of which RBP is likely to be driving the observed changes in the splicing.

To perform this analysis it is necessary to execute the following steps:

• Create a matrix that relates AS events with RBP binding sites using the function EventPointer:::CreateExSmatrix
• Run the prediction algorithm using the function EventPointer:::SF_Prediction

A detailed explanation of how to use this pipeline is desribed in https://github.com/clobatofern/SFPointer_testPipeline.