Alignment preparation

In order to avoid version incompatibilities with hisat2 aligner and the installation of the required tools, pre-compiled versions will be used to run the R functions. Consequently, all the tools were compressed and stored in an external server to be downloaded whenever it is required (http://iwbbio.ugr.es/utils/unixUtils.tar.gz). If the tools are directly downloaded from the link, the compressed files must be decompressed in the current project folder in R or RStudio. The name of the resultant folder must be “utils”. Nevertheless, this file can be downloaded automatically by just calling the function rawAlignment, in case the folder utils is not detected in the project folder. This is all needed to run hisat2 through the function rawAlignment. It is not possible to run the alignment without the utils folder. It must be mentioned too that the different files included in the compressed .tar.gz are not only the aligner but also functions needed in the raw alignment process. The tools included are the following:

• Bowtie2 [7].
• Hisat2 [8].
• Htseq-count [9].
• Samtools [12].
• Sratoolkit [13].
• GDC-client.

Launching Raw Alignment step

The rawAlignment function allows running hisat2 aligner. The function takes as single input a CSV from GEO or ArrayExpress loaded in R. There is the possibility to process data from GDC-portal, but a previous authorization (token file) from this platform is required. Furthermore, there is a set of logical parameters to edit the default pipeline followed for the function. With the parameters the user can select if the BAM/SAM/Count files are created. The user can choose if wants to download the reference genome, the GTF, and which version. Even if the user has custom FASTA and GTF files, this can be specified by setting the parameter referenceGenome to “custom” and using the parameters customFA and customGTF to indicates the paths to the custom files. Other functionality is the possibility to process BAM files from the GDC Portal database by setting to TRUE the parameter fromGDC. Then the function will download the specific genome reference of GDC and process the BAM files to Count files. Furthermore, if the user has access to the controlled data, with the token and the manifest acquired from GDC Portal web platform, the samples can be downloaded automatically. An example to run the function with hisat2 aligner is showed below:

# Downloading one series from NCBI/GEO and one series from ArrayExpress
downloadPublicSeries(c("GSE74251"))

# Using read.csv for NCBI/GEO files (read.csv2 for ArrayExpress files)
GSE74251csv <- read.csv("ReferenceFiles/GSE74251.csv")

# Performing the alignment of the samples by using hisat2 aligner 
rawAlignment(GSE74251csv,downloadRef=TRUE,downloadSamples=TRUE,BAMfiles = TRUE,
SAMfiles = TRUE,countFiles = TRUE,referenceGenome = 38, fromGDC = FALSE, customFA = "", 
customGTF = "", tokenPath = "", manifest = "")

RawAlignment function creates a folder structure in the current project folder which will store all the downloaded and created files. The main folder of this structure is the folder ReferenceFiles but inside of it there are more folders that allows storing the different files used by the process in an organized way.

Another important requirement to take into account is the format of the csv file used to launch the function. It could be from three repositories, two publics (NCBI/GEO and ArrayExpress) and one controlled (GDC Portal). Each of these repositories has its own format in the csv file that contains the information to download and process the desired samples. The necessary format for each repository is explained below.

# GDC portal controlled data processing with automatic raw data download

rawAlignment(x, downloadRef=TRUE, downloadSamples=TRUE, fromGDC = TRUE,
tokenPath = "~/pathToToken", manifestPath = "~/pathToManifest")

Preparing the count files

From now on, the data that will be used for the documentation are real count files, but with a limited number of genes (around 1000). Furthermore, to reduce the computational cost of this example, only 5 samples from each of the two selected series will be taken into account. Showed in the code snippet below, two RNA-seq series from NCBI/GEO are downloaded automatically and the existing count files prepared to be merged in one matrix with the purpose of preparing the data for further steps:

suppressMessages(library(KnowSeq))

dir <- system.file("extdata", package="KnowSeq")

# Using read.csv for NCBI/GEO files and read.csv2 for ArrayExpress files

GSE74251 <- read.csv(paste(dir,"GSE74251.csv",sep = "/"))
GSE81593 <- read.csv(paste(dir,"GSE81593.csv",sep = "/"))

# Creating the csv file with the information about the counts files location and the labels
Run <- GSE74251$Run
Path <- paste(dir,"/countFiles/",GSE74251$Run,sep = "")
Class <- rep("Tumor", length(GSE74251$Run))
GSE74251CountsInfo <-  data.frame(Run = Run, Path = Path, Class = Class)

Run <- GSE81593$Run
Path <- paste(dir,"/countFiles/",GSE81593$Run,sep = "")
Class <- rep("Control", length(GSE81593$Run))
GSE81593CountsInfo <-  data.frame(Run = Run, Path = Path, Class = Class)

mergedCountsInfo <- rbind(GSE74251CountsInfo, GSE81593CountsInfo)

write.csv(mergedCountsInfo, file = "mergedCountsInfo.csv")

However, the user can run a complete example by executing the following code:

dir <- system.file("script", package="KnowSeq")

# Code to execute the example script
source(paste(dir,"/KnowSeqExample.R",sep=""))


# Code to edit the example script
file.edit(paste(dir,"/KnowSeqExample.R",sep=""))

Processing count files

After the raw alignment step, a list of count files of the samples is available at ReferenceFiles/Samples/RNAseq/CountFiles. The next step in the pipeline implemented in this package is the processing of those count files in order to obtain a gene expression matrix by merging all of them.

# Merging in one matrix all the count files indicated inside the CSV file

countsInformation <- countsToMatrix("mergedCountsInfo.csv", extension = "count")

## 
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR2753159/SRR2753159.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR2753160/SRR2753160.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR2753161/SRR2753161.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR2753162/SRR2753162.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR2753163/SRR2753163.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR3541296/SRR3541296.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR3541297/SRR3541297.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR3541298/SRR3541298.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR3541299/SRR3541299.count
## /tmp/RtmpHkmplp/Rinst1de01bde14c5/KnowSeq/extdata/countFiles/SRR3541300/SRR3541300.count
## Merging 10 counts files...

# Exporting to independent variables the counts matrix and the labels

countsMatrix <- countsInformation$countsMatrix

labels <- countsInformation$labels

# Downloading human annotation
myAnnotation <- getGenesAnnotation(rownames(countsMatrix))

## Getting annotation of the Homo Sapiens...
## Using reference genome 38.

# Downloading mus musculus annotation
myAnnotationMusMusculus <- getGenesAnnotation(rownames(countsMatrix), 
notHSapiens = TRUE,notHumandataset = "mm129s1svimj_gene_ensembl")

## Downloading annotation mm129s1svimj_gene_ensembl...

## 
## Connection error, trying again...
## 
## Connection error, trying again...

# Calculating gene expression values matrix using the counts matrix

expressionMatrix <- calculateGeneExpressionValues(countsMatrix,myAnnotation,
genesNames = TRUE)

## Calculating gene expression values...
## RQ fit ..........
## SQN

## .

# Plotting the boxplot of the expression of each samples for all the genes

dataPlot(expressionMatrix,labels,mode = "boxplot", toPNG = TRUE,
toPDF = TRUE)

## Creating PNG...

## Creating PDF...

Batch effect removal

A crucial step in this pipeline is the batch effect treatment. It is widely known that this is a crucial step in the omics data processing due to the intrinsic deviations that the data can present due to its origin, sequencing design, etc… Besides, when working with public data it is very difficult to know if exists a real batch effect among the selected datasets. This package allows removing batch effect if the batch groups are known by calling the function batchEffectRemoval, that makes use of sva package [18], with the parameter mode equal to “combat” [19]. This step allows obtaining an expression matrix with the batch effect treated by combat method. An example to do this is below:

# Removing batch effect by using combat and known batch groups 

batchGroups <- c(1,1,1,1,2,2,1,2,1,2)

expressionMatrixCorrected <- batchEffectRemoval(expressionMatrix, labels, 
batchGroups = batchGroups, method = "combat")

## Correcting batch effect by using combat method...

## Using the 'mean only' version of ComBat

## Found2batches

## Adjusting for0covariate(s) or covariate level(s)

## Standardizing Data across genes

## Fitting L/S model and finding priors

## Finding parametric adjustments

## Adjusting the Data

There is another method in the function that removes the batch effect that uses surrogate variable analysis or sva. The only requirement to use it is to set the parameter method equal to “sva”. This method returns a matrix with the batch effect corrected which has to be used as input of the function DEGsExtraction.

# Calculating the surrogate variable analysis to remove batch effect

expressionMatrixCorrected <- batchEffectRemoval(expressionMatrix, labels, method = "sva")

## Calculating sva model to batch effect correction...
## Number of significant surrogate variables is:  1 
## Iteration (out of 5 ):1  2  3  4  5

Differential Expressed Genes extraction and visualization

There is a long way between the raw data and the DEGs extraction, for that in this step the samples have to have had a strong pre-processing step applied. At this point of the pipeline the DEGs existing among two or more classes will be extracted by using the novel parameter coverage (cov) along with limma R-bioc package [20]. The parameter cov represents the number of different pathologies that a certain gen is able to discern. By default, the parameter is set to 1, so all genes that has the capability to discern among the comparison of two classes would be selected as DEGs. To understand better this parameter, our multiclass study applied to different leukemia sub-types introduces it, and it’s publicly available [3].

The function DEGsExtraction receives an expression matrix, the labels of the samples and the restriction imposed for considering a gene as differential expressed gene. The function returns a list containing the table with statistical values of each DEGs and the expression matrix of the DEGs instead all of the genes. The call to the function is listed below:

# Extracting DEGs that pass the imposed restrictions

DEGsInformation <- DEGsExtraction(expressionMatrixCorrected, labels,
lfc = 1.0, pvalue = 0.01, number = 100, cov = 1)

## Two classes detected, applying limma biclass

topTable <- DEGsInformation$DEG_Results$DEGs_Table

DEGsMatrix <- DEGsInformation$DEG_Results$DEGs_Matrix

DEGs are genes that have a truly different expression among the studied classes, for that it is important to try to see graphically if those DEGs comply with this requirement. In order to provide a tool to perform this task, the function dataPlot encapsulate a set of graphs that allows plotting in different ways the expression of the DEGs.

dataPlot function also allows representing an ordered boxplot that internally orders the samples by class and plots a boxplot for each samples and for the first top 12 DEGs in this example. With this plot, the difference at gene expression level between the classes can be seen graphically. The code to reproduce this plot is the following:

# Plotting the expression of the first 12 DEGs for each of the samples in an ordered way

dataPlot(DEGsMatrix[1:12,],labels,mode = "orderedBoxplot",toPNG = FALSE,toPDF = FALSE)

In the previous boxplot the expression of a set of DEGs for each sample its showed, however it is interesting to see the differentiation at gene expression level for each of the top 12 genes used before separately. It is recommended to use this function with a low number of genes, because with a larger number the plot it is difficult to distinguish the information provided and R would not have enough memory to calculate the plot. For that, the function dataPlot with the mode genesBoxplot allows to do that by executing the next code:

# Plotting the expression of the first 12 DEGs separatelly for all the samples

dataPlot(DEGsMatrix[1:12,],labels,mode = "genesBoxplot",toPNG = FALSE,toPDF = FALSE)

Finally, it is possible to plot one of the most widespread visualization methods in the literature, the heatmap. By setting the parameter method to heatmap, the function calculates the heatmap for the given samples and classes. The code to do this is the same than for the previous boxplot but changing the method parameter:

# Plotting the heatmap of the first 12 DEGs separatelly for all the samples

dataPlot(DEGsMatrix[1:12,],labels,mode = "heatmap",toPNG = FALSE,toPDF = FALSE)

Performing the machine learning processing: classifier design and assessment and gene selection

Normally, in the literature, the last step in the pipeline for differential gene expression analysis is the DEGs extraction step. However, in this package a novel machine learning step is implemented with the purpose of giving to the user an automatic tool to assess the DEGs, and evaluate their robustness in the discernment among the studied pathologies. This library has three possible classification methodologies to take into account. These options are k-NN [21], SVM [22] and Random Forest [23], three of the most popular classifiers in the literature. Furthermore, it includes two different working procedures for each of them. The first one implements a cross-validation process, in order to assess the expected accuracy with different models and samples the DEGs with a specific number of folds. These functions return a list with 4 objects that contain the confusion matrices, the accuracy, the sensitivity and the specificity.

The second one is to assess a specific test dataset by using a classifier trained using the training dataset separately. Moreover, the function featureSelection allows performing a feature selection process by using, mRMR [24], Random Forest (as feature selector instead of classifier) or Disease Association based algorithms with the purpose of finding the best DEGs order to assess the data. Da-FS is a new novel method designed in KnowSeq with the purpose of giving to the expert a biological based feature selection method. This method makes use of targetValidation webplatform to acquire an association score for each DEGs with the required genetic disease, breast cancer in the example. This score takes values between 1 and 0, meaning 1 a total association and 0 no association. Therefore, the DEGs are sorted by this score, achieving a ranking in which in the first positions those DEGs with more biological relation to breast cancer are placed.

Moreover, targetValidation webplatform allows acquiring evidences that tie a gene with a certain disease. DA-RED-FS is a novel iterative method based on DA-FS that use these evidences to calculate the redundance between genes based on biological information. This redundances take values between 1 and 0. For example, if gene A has a redundance of 1 with gene B, means that all found evidences for gene A and a certain disease are also found in gene B and the disease. Likewise, if this redundance is 0, means that there are not any gene A evidences in gene B evidences. DA-RED-FS starts with an empty set of selected genes, S_G, and a set of possible genes G. In the first step, the gene with the highest DA score is selected. In the following steps, genes that verify the following equation are added to selected genes set.

$$ max_{g \in G - S_G} DA(g) - \frac{\sum_{g_i \in S_g} RED(g,g_i)}{|S_G|} \times DA(g) $$

Where RED(g, g_i) is the redundance between gene g and gen g_i and DA(g) is the DA score of gene g. The algorithms ends when a certain number of genes are selected, this number can be fixed by the maxGenes parameter.

To invoke these functions, it is necessary an expression matrix with the samples in the rows and the genes in the columns and the labels of the samples, the genes that will be assessed and the number of fold in the case of the cross-validation function. In the case of the test functions, it is necessary the matrix and the labels for both the training and the test datasets:

DEGsMatrixML <- t(DEGsMatrix)

# Feature selection process with mRMR and RF 
mrmrRanking <- featureSelection(DEGsMatrixML,labels,colnames(DEGsMatrixML), mode = "mrmr")

## Calculating the ranking of the most relevant genes by using mRMR algorithm...
## mRMR ranking: BCAS1 ATP2B4 RPS17P5 CXorf56 GALC VIM SCIN ROS1 SH2D2A ETV1 AGPS SNAI2 LAMC2 GNA15 MAGEC2 EHD3 GYG2 PRSS3 FSTL4 BIRC3 ADAM22 COL17A1 CALCR CCN5 ARHGAP44 DNASE1L1 TNC DLX3 FUT8 CCDC85A CDH1 EHD2 HOXC8 PRSS21 MATK ABCB4 DKK3 YBX2 COL9A2 RIPOR3 VCAN NFIX LY75 SLC7A2 CCT8L1P ERBB3 ARSD BARX2 APPBP2 CNTN1 DCN NLRP2 APBA2 TYMP MPPED2 PRKCQ TENM1 FUZ CHDH NEXMIF PLAUR CYP26B1 DLEC1 E2F2 FAS RAI14 EPHA3 AGPAT4 RCN1 ATP2C2 GAS7 CPS1 LAMA3 ARHGAP31 TIMP2 PRSS8 DGKA PLEKHG6 TRAF1 BTN3A1 EPN3 CP ABCC2 PRSS22 ME1 ZIC2 SYT13 SEZ6 PLEKHB1 CYFIP2 ABCC8 ZMYND12 SPATA20 CEACAM7 DAPK2 GRAMD1B IYD PHF21B SARM1 NDC1

rfRanking <- featureSelection(DEGsMatrixML,labels,colnames(DEGsMatrixML), mode = "rf")

## Calculating the ranking of the most relevant genes by using Random Forest algorithm...
## Random Forest ranking: HOXC8 DKK3 RIPOR3 CEACAM7 ATP2B4 GYG2 ABCC8 LAMA3 CHDH EPHA3 ATP2C2 FAS FSTL4 GALC DAPK2 PHF21B SPATA20 IYD YBX2 BIRC3 E2F2 GNA15 AGPAT4 PRSS3 ETV1 MAGEC2 ARHGAP44 SCIN CCDC85A EHD2 ZIC2 MPPED2 FUZ PLEKHG6 COL17A1 ARHGAP31 PLEKHB1 NFIX LY75 GRAMD1B ADAM22 VIM SYT13 PLAUR TIMP2 MATK CYP26B1 BCAS1 ROS1 DLEC1 RPS17P5 EHD3 PRSS21 VCAN ARSD DGKA CPS1 ABCC2 CXorf56 EPN3 ERBB3 ZMYND12 AGPS NLRP2 CP PRSS22 LAMC2 FUT8 SLC7A2 CNTN1 NEXMIF DNASE1L1 ABCB4 CCT8L1P SEZ6 CCN5 TNC DLX3 RCN1 CDH1 CYFIP2 BTN3A1 SNAI2 NDC1 PRSS8 ME1 COL9A2 APPBP2 DCN APBA2 PRKCQ TENM1 RAI14 BARX2 TYMP GAS7 SH2D2A CALCR SARM1 TRAF1

daRanking <- featureSelection(DEGsMatrixML,labels,colnames(DEGsMatrixML), mode = "da", disease = "breast")

## Calculating ranking of biological relevant genes by using DA implementation...
## Getting annotation of the Homo Sapiens...
## Using reference genome 38.
## Obtaining scores for breast...
## Disease scores acquired successfully!
## Disease Association ranking:    BCAS1  RPS17P5     GALC      VIM     SCIN     ROS1   SH2D2A     ETV1 
##        0        0        0        0        0        0        0        0 
##     AGPS    SNAI2    LAMC2    GNA15   MAGEC2     EHD3     GYG2    PRSS3 
##        0        0        0        0        0        0        0        0 
##    FSTL4    BIRC3   ADAM22  COL17A1    CALCR     CCN5 ARHGAP44 DNASE1L1 
##        0        0        0        0        0        0        0        0 
##      TNC     DLX3   ATP2B4  PLEKHB1     FUT8  CCDC85A     CDH1     EHD2 
##        0        0        0        0        0        0        0        0 
##    HOXC8   PRSS21     MATK    ABCB4     DKK3     YBX2   COL9A2   RIPOR3 
##        0        0        0        0        0        0        0        0 
##     VCAN     NFIX    ABCC8     LY75   SLC7A2  CCT8L1P    ERBB3     ARSD 
##        0        0        0        0        0        0        0        0 
##    BARX2   APPBP2    CNTN1      DCN    NLRP2  CXorf56    APBA2    SARM1 
##        0        0        0        0        0        0        0        0 
##     TYMP   MPPED2    PRKCQ    TENM1      FUZ     CHDH   NEXMIF    PLAUR 
##        0        0        0        0        0        0        0        0 
##  CYP26B1    DLEC1     E2F2      FAS    RAI14    EPHA3   AGPAT4     RCN1 
##        0        0        0        0        0        0        0        0 
##   ATP2C2     GAS7     CPS1    LAMA3     NDC1  CEACAM7    DAPK2  SPATA20 
##        0        0        0        0        0        0        0        0 
## ARHGAP31  ZMYND12    TIMP2    PRSS8    SYT13     DGKA  PLEKHG6    TRAF1 
##        0        0        0        0        0        0        0        0 
##   CYFIP2   BTN3A1     SEZ6     EPN3      IYD       CP    ABCC2   PRSS22 
##        0        0        0        0        0        0        0        0 
##   PHF21B  GRAMD1B      ME1     ZIC2 
##        0        0        0        0

# CV functions with k-NN, SVM and RF
results_cv_knn <- knn_trn(DEGsMatrixML,labels,names(mrmrRanking)[1:10],5)

## Tuning the optimal K...

## Loading required package: ggplot2

## Loading required package: lattice

## Optimal K: 7 
## Running K-Fold Cross-Validation...
## Training fold 1...
## Running K-Fold Cross-Validation...
## Training fold 2...
## Running K-Fold Cross-Validation...
## Training fold 3...
## Running K-Fold Cross-Validation...
## Training fold 4...
## Running K-Fold Cross-Validation...
## Training fold 5...
## Classification done successfully!

results_cv_svm <- svm_trn(DEGsMatrixML,labels,rfRanking[1:10],5)

## Tuning the optimal C and G...

## Optimal cost: 0.25 
## Optimal gamma: 0.5 
## Training fold 1...
## Training fold 2...
## Training fold 3...
## Training fold 4...
## Training fold 5...
## Classification done successfully!

results_cv_rf <- rf_trn(DEGsMatrixML,labels,names(daRanking)[1:10],5)

## Tuning the optimal mtry...

## Training fold 1...

## Training fold 2...

## Training fold 3...

## Training fold 4...

## Training fold 5...

## Classification done successfully!

It is important to show graphically the results of the classifiers and for that purpose, the function dataPlot implements some methods. Concretely, to plot the accuracy, the sensitivity or the specificity reached by the classifiers, the function dataPlot has to be run with the parameter method equal to classResults. This method generated as many random colors as folds or simulations in the rows of the matrix passed to the function but, through the parameter colors a vector of desired colors can be specified. For the legend, the function uses the rownames of the input matrix but these names can be changed with the parameter legend. An example of this method is showed below:

# Plotting the accuracy of all the folds evaluated in the CV process

dataPlot(rbind(results_cv_knn$accuracy$meanAccuracy,results_cv_knn$accuracy$standardDeviation),mode = "classResults", legend = c("Mean Acc","Standard Deviation"),
main = "Mean Accuracy with k-NN", xlab = "Genes", ylab = "Accuracy")

# Plotting the sensitivity of all the folds evaluated in the CV process

dataPlot(rbind(results_cv_knn$sensitivity$meanSensitivity,results_cv_knn$sensitivity$standardDeviation),mode = "classResults", legend = c("Mean Sens","Standard Deviation"),
main = "Mean Sensitivity with k-NN", xlab = "Genes", ylab = "Sensitivity")

# Plotting the specificity of all the folds evaluated in the CV process

dataPlot(rbind(results_cv_knn$specificity$meanSpecificity,results_cv_knn$specificity$standardDeviation),mode = "classResults", legend = c("Mean Spec","Standard Deviation"),
main = "Mean Specificity with k-NN", xlab = "Genes", ylab = "Specificity")

# Plotting all the metrics depending on the number of used DEGs in the CV process

dataPlot(results_cv_knn,labels,mode = "heatmapResults")

Furthermore, the function dataPlot counts with another similar mode to the previous but this time to represents confusion matrices. This mode is called confusionMatrix and allows creating graphically a confusion matrix with the most important statistical measures. The following code allows doing this:

# Plotting the confusion matrix with the sum of the confusion matrices 
# of each folds evaluated in the CV process

allCfMats <- results_cv_knn$cfMats[[1]]$table + results_cv_knn$cfMats[[2]]$table + 
results_cv_knn$cfMats[[3]]$table + results_cv_knn$cfMats[[4]]$table + 
results_cv_knn$cfMats[[5]]$table

dataPlot(allCfMats,labels,mode = "confusionMatrix")

# Test functions with k-NN, SVM and RF
trainingMatrix <- DEGsMatrixML[c(1:4,6:9),]
trainingLabels <- labels[c(1:4,6:9)]
testMatrix <- DEGsMatrixML[c(5,10),]
testLabels <- labels[c(5,10)]

results_test_knn <- knn_test(trainingMatrix, trainingLabels, testMatrix,
testLabels, names(mrmrRanking)[1:10], bestK = results_cv_knn$bestK)

## Testing with 1 variables...
## Testing with 2 variables...
## Testing with 3 variables...
## Testing with 4 variables...
## Testing with 5 variables...
## Testing with 6 variables...
## Testing with 7 variables...
## Testing with 8 variables...
## Testing with 9 variables...
## Testing with 10 variables...
## Classification done successfully!

results_test_svm <- svm_test(trainingMatrix, trainingLabels, testMatrix,
testLabels, rfRanking[1:10], bestParameters = results_cv_svm$bestParameters)

## Testing with 1 variables...
## Testing with 2 variables...
## Testing with 3 variables...
## Testing with 4 variables...
## Testing with 5 variables...
## Testing with 6 variables...
## Testing with 7 variables...
## Testing with 8 variables...
## Testing with 9 variables...
## Testing with 10 variables...
## Classification done successfully!

results_test_rf <- rf_test(trainingMatrix, trainingLabels, testMatrix,
testLabels, colnames(DEGsMatrixML)[1:10], results_cv_rf$bestParameters)

## Testing with 1 variables...

## Testing with 2 variables...

## Testing with 3 variables...

## Testing with 4 variables...

## Testing with 5 variables...

## Testing with 6 variables...

## Testing with 7 variables...

## Testing with 8 variables...

## Testing with 9 variables...

## Testing with 10 variables...

## Classification done successfully!

# Plotting the accuracy achieved in the test process

dataPlot(results_test_knn$accVector,mode = "classResults",
main = "Accuracy with k-NN", xlab = "Genes", ylab = "Accuracy")

dataPlot(results_test_svm$accVector,mode = "classResults",
main = "Accuracy with SVM", xlab = "Genes", ylab = "Accuracy")

dataPlot(results_test_rf$accVector,mode = "classResults",
main = "Accuracy with RF", xlab = "Genes", ylab = "Accuracy")

Gene Ontology

Gene ontology (GO) provides information about the biological functions of the genes. In order to complete this pipeline, it is important to know if the DEGs have functions related with the studied pathologies. In this sense, this package brings the possibility to know the GOs from the three different ontologies (BP, MF and CC) by using the function geneOntologyEnrichment that internally uses information from DAVID web-platform [20]. The function returns a list that contains a matrix for each ontology and a matrix with the GOs of the three ontologies together. Moreover, the matrices have different statistical measures and the description of the functionality of each GO.

# Retrieving the GO information from the three different ontologies

GOsList <- geneOntologyEnrichment(names(mrmrRanking)[1:10], geneType='GENE_SYMBOL', pvalCutOff=0.1)

## Getting gene symbols...Getting annotation of the Homo Sapiens...
## Using reference genome 38.
## Retrieving Gene Ontology terms related to the list of DEGs...

## Error in geneOntologyEnrichment(names(mrmrRanking)[1:10], geneType = "GENE_SYMBOL", : Error in request. Please check parameters or use geneType as ENTREZ_GENE_ID or GENE_SYMBOL

For example, in this example, the top 10 GOs from the BP ontology for the extracted DEGs are shown in the following image.

Pathway Information

Another important step in the enrichment methodology in this pipeline is the related pathway information. The function uses the DEGs to retrieve pathways with any relation with DEGs. For that, DEGsToPathways makes use of the information gathered from KEGG database in order to supply as more details as possible .

# Retrieving the pathways related to the DEGs
paths <- DEGsToPathways(names(mrmrRanking)[1:10])

## Retrieving information about KEGG pathways...

print(paths)

##    KEGG_Path         Name  Description        Class  Genes
## 1   map04020 Calcium .... Ca2+ tha.... Environm.... ATP2B4
## 2   map04022 cGMP-PKG.... Cyclic G.... Environm.... ATP2B4
## 3   map04024 cAMP sig.... cAMP is .... Environm.... ATP2B4
## 4   map04261 Adrenerg.... Cardiac .... Organism.... ATP2B4
## 5   map04925 Aldoster.... Aldoster.... Organism.... ATP2B4
## 6   map04961 Endocrin.... Calcium .... Organism.... ATP2B4
## 7   map04970 Salivary.... Saliva h.... Organism.... ATP2B4
## 8   map04972 Pancreat.... The panc.... Organism.... ATP2B4
## 9   map04978 Mineral .... Minerals.... Organism.... ATP2B4
## 10  map00600      nothing      nothing Metaboli....   GALC
## 11  map01100      nothing      nothing      nothing   GALC
## 12  map04142     Lysosome Lysosome.... Cellular....   GALC
## 13  map05169 Epstein-.... Epstein-.... Human Di....    VIM
## 14  map05206 MicroRNA.... MicroRNA.... Human Di....    VIM
## 15  map04666 Fc gamma.... Phagocyt.... Organism....   SCIN
## 16  map04810      nothing      nothing Cellular....   SCIN
## 17  map05203 Viral ca.... There is.... Human Di....   SCIN
## 18  map04370 VEGF sig.... There is.... Environm.... SH2D2A
## 19  map05202 Transcri.... In tumor.... Human Di....   ETV1

Related diseases

Finally, the last enrichment method implemented is the related diseases enrichment. In this step, the function DEGsToDisease searches the diseases related to a list of genes or DEGs indicated as parameter. The function returns a list of diseases only for genes and also for groups of genes with several statistical values to know the relation between the diseases and the gene or group of genes. This information is retrieved from the targetValidation web platform [27]. Also, it is possible to obtain the evidences that tie diseases with each gene from targetValidation platform, setting getEvidences parameter to TRUE. There are seven types of evidences: known drug, literature, affected pathway, animal model, RNA expression, genetic association or somatic mutation. Each of them contains different fields with the most important information for each one.

# Downloading the information about the DEGs related diseases

diseasesTargetValidation <- DEGsToDiseases(rownames(DEGsMatrix)[1:5], getEvidences = FALSE)

## Getting annotation of the Homo Sapiens...
## Using reference genome 38.
## Obtaining related diseases with the DEGs from targetValidation platform...
## Diseases acquired successfully!

diseasesTargetValidation

## $BCAS1
## $BCAS1$summary
##       Disease                       Overall Score       Literature RNA Expr.
##  [1,] "Red cell distribution width" "0.469297159128405" "0"        "0"      
##  [2,] "creatinine measurement"      "0.450835042093441" "0"        "0"      
##  [3,] "erythrocyte count"           "0.447028983031038" "0"        "0"      
##  [4,] "urolithiasis"                "0.441888406496364" "0"        "0"      
##  [5,] "hair color"                  "0.437874508196795" "0"        "0"      
##  [6,] "mean corpuscular volume"     "0.425908260615748" "0"        "0"      
##  [7,] "mean platelet volume"        "0.405252455375351" "0"        "0"      
##  [8,] "hair colour measurement"     "0.403632508913504" "0"        "0"      
##  [9,] "uric acid measurement"       "0.399312013686817" "0"        "0"      
## [10,] "vitamin D measurement"       "0.390952149137878" "0"        "0"      
##       Genetic Assoc.      Somatic Mut. Known Drug Animal Model
##  [1,] "0.771958191577288" "0"          "0"        "0"         
##  [2,] "0.741589410940584" "0"          "0"        "0"         
##  [3,] "0.7353287327888"   "0"          "0"        "0"         
##  [4,] "0.726872874729183" "0"          "0"        "0"         
##  [5,] "0.72027031681414"  "0"          "0"        "0"         
##  [6,] "0.700586748177613" "0"          "0"        "0"         
##  [7,] "0.666609517016521" "0"          "0"        "0"         
##  [8,] "0.663944828094342" "0"          "0"        "0"         
##  [9,] "0.656837941515046" "0"          "0"        "0"         
## [10,] "0.643086599122487" "0"          "0"        "0"         
##       Affected Pathways
##  [1,] "0"              
##  [2,] "0"              
##  [3,] "0"              
##  [4,] "0"              
##  [5,] "0"              
##  [6,] "0"              
##  [7,] "0"              
##  [8,] "0"              
##  [9,] "0"              
## [10,] "0"              
## 
## 
## $GALC
## $GALC$summary
##       Disease                                  Overall Score      
##  [1,] "Krabbe disease"                         "0.860359048265709"
##  [2,] "genetic disorder"                       "0.508468930000455"
##  [3,] "body mass index"                        "0.505583386717068"
##  [4,] "basophil count"                         "0.423572575869943"
##  [5,] "inflammatory bowel disease"             "0.420334782589922"
##  [6,] "spastic ataxia"                         "0.406537840153172"
##  [7,] "Crohn's disease"                        "0.380355006681393"
##  [8,] "adult Krabbe disease"                   "0.378987509885701"
##  [9,] "infantile Krabbe disease"               "0.37336804819522" 
## [10,] "late-infantile/juvenile Krabbe disease" "0.369579854684702"
##       Literature           RNA Expr. Genetic Assoc.      Somatic Mut.
##  [1,] "0.569144483451211"  "0"       "0.947373469444716" "0"         
##  [2,] "0.0121586159522324" "0"       "0.835784849430134" "0"         
##  [3,] "0"                  "0"       "0.831646280634826" "0"         
##  [4,] "0"                  "0"       "0.69674472412655"  "0"         
##  [5,] "0"                  "0"       "0.691418800036603" "0"         
##  [6,] "0"                  "0"       "0.668723877372783" "0"         
##  [7,] "0"                  "0"       "0.625655104455465" "0"         
##  [8,] "0"                  "0"       "0.607930797611621" "0"         
##  [9,] "0.124625813510382"  "0"       "0.607930797611621" "0"         
## [10,] "0"                  "0"       "0.607930797611621" "0"         
##       Known Drug Animal Model        Affected Pathways
##  [1,] "0"        "0"                 "0"              
##  [2,] "0"        "0"                 "0"              
##  [3,] "0"        "0"                 "0"              
##  [4,] "0"        "0"                 "0"              
##  [5,] "0"        "0"                 "0"              
##  [6,] "0"        "0"                 "0"              
##  [7,] "0"        "0"                 "0"              
##  [8,] "0"        "0.309497569064076" "0"              
##  [9,] "0"        "0"                 "0"              
## [10,] "0"        "0"                 "0"              
## 
## 
## $VIM
## $VIM$summary
##       Disease                                            Overall Score      
##  [1,] "Partial congenital cataract"                      "0.456330257326305"
##  [2,] "Developmental cataract"                           "0.434644159966231"
##  [3,] "pulverulent cataract"                             "0.369579854684702"
##  [4,] "aspartate aminotransferase measurement"           "0.365287737275524"
##  [5,] "total cholesterol measurement"                    "0.353072224839048"
##  [6,] "high density lipoprotein cholesterol measurement" "0.324873530627104"
##  [7,] "low density lipoprotein cholesterol measurement"  "0.297570315868196"
##  [8,] "glioblastoma multiforme"                          "0.287266840421116"
##  [9,] "hypertension"                                     "0.274481885258201"
## [10,] "head and neck squamous cell carcinoma"            "0.272128616119786"
##       Literature           RNA Expr.           Genetic Assoc.      Somatic Mut.
##  [1,] "0"                  "0"                 "0.750628622729908" "0"         
##  [2,] "0.0287078432205488" "0"                 "0.571454949754924" "0"         
##  [3,] "0"                  "0"                 "0.607930797611621" "0"         
##  [4,] "0"                  "0"                 "0.600870590387312" "0"         
##  [5,] "0"                  "0"                 "0.580776999991057" "0"         
##  [6,] "0"                  "0"                 "0.534392289226727" "0"         
##  [7,] "0"                  "0"                 "0.489480574166107" "0"         
##  [8,] "0.745817510119163"  "0.210149679447545" "0"                 "0"         
##  [9,] "0.0303965398805811" "0"                 "0.449982031994787" "0"         
## [10,] "0.736957329136789"  "0"                 "0"                 "0"         
##       Known Drug Animal Model Affected Pathways  
##  [1,] "0"        "0"          "0"                
##  [2,] "0"        "0"          "0"                
##  [3,] "0"        "0"          "0"                
##  [4,] "0"        "0"          "0"                
##  [5,] "0"        "0"          "0"                
##  [6,] "0"        "0"          "0"                
##  [7,] "0"        "0"          "0"                
##  [8,] "0"        "0"          "0.861142542095881"
##  [9,] "0"        "0"          "0"                
## [10,] "0"        "0"          "0.821566086980618"
## 
## 
## $SCIN
## $SCIN$summary
##       Disease                           Overall Score       
##  [1,] "body height"                     "0.356486429466707" 
##  [2,] "sitting height measurement"      "0.205806121295621" 
##  [3,] "metabolite measurement"          "0.167857497464434" 
##  [4,] "bone density"                    "0.113325829545007" 
##  [5,] "urate measurement"               "0.113325829545007" 
##  [6,] "adolescent idiopathic scoliosis" "0.0563883824115441"
##  [7,] "neoplasm"                        "0.0491782008359755"
##  [8,] "lymphocyte count"                "0.0429726846885094"
##  [9,] "hemoglobin measurement"          "0.0425605722167021"
## [10,] "hematocrit"                      "0.040446582063385" 
##       Literature          RNA Expr. Genetic Assoc.       Somatic Mut.
##  [1,] "0"                 "0"       "0.586393107352409"  "0"         
##  [2,] "0"                 "0"       "0.338535442034145"  "0"         
##  [3,] "0"                 "0"       "0.276112837388558"  "0"         
##  [4,] "0"                 "0"       "0.18641238441979"   "0"         
##  [5,] "0"                 "0"       "0.18641238441979"   "0"         
##  [6,] "0"                 "0"       "0.0927546073222105" "0"         
##  [7,] "0.404472030609257" "0"       "0"                  "0"         
##  [8,] "0"                 "0"       "0.0706868032633587" "0"         
##  [9,] "0"                 "0"       "0.070008909540214"  "0"         
## [10,] "0"                 "0"       "0.0665315562598368" "0"         
##       Known Drug Animal Model Affected Pathways
##  [1,] "0"        "0"          "0"              
##  [2,] "0"        "0"          "0"              
##  [3,] "0"        "0"          "0"              
##  [4,] "0"        "0"          "0"              
##  [5,] "0"        "0"          "0"              
##  [6,] "0"        "0"          "0"              
##  [7,] "0"        "0"          "0"              
##  [8,] "0"        "0"          "0"              
##  [9,] "0"        "0"          "0"              
## [10,] "0"        "0"          "0"