The GeneNetwork pipeline is a command line python3 tool to perform the decomposition and gene pathway score calculation. The main goal of the analysis is to calculate gene z-scores for every gene for every pathway in the input files which gives the contribution of a gene to the pathway. Where the input pathway file contains only a Boolean value if the gene is linked to the pathway, this analysis calculates a contribution z-score value for each gene to every pathway by using gene expression data.

The entire method and clinical application of this pipeline is described in our article “KidneyNetwork: Using kidney derived gene expression data to predict and prioritize known and novel genes involved in kidney disease” (under review). Below you find a summary of the method and a description how the pipeline can be used.

The pipelines contain two steps: 1) the decomposition and 2) the gene-phenotype score calculation. The two steps are separated in two different python scripts.

The first step of the pipeline is the calculation of the eigenvectors from the gene expression dataset. Three different decomposition methods are implemented, PCA, FastICA and a stabilized version of FastICA. The decomposition method PCA was implemented by using Singular value decomposition class of Sklearn package and the svd package of SciPy. The different algorithms van be selected by using the “pca_type” parameter. See parameter description for more information.

The input of the script in a tab separated gene expression file where the rows represented the samples and the columns the genes. The output contains one file with the gene/component’s matrix (eigenvectors) and one file with the sample/component’s matrix (projected data). If the “input_transpose” option was set, both output files are swapped to maintain the structure and layout of both files. Beside the output matrixes a log file is created.

The output is by default exported to both txt and pickle format. The .txt file is gzip compressed. If some files are not necessary, they can be disabled by using the output parameters.

The second step is the calculation of the gene-pathway scores based on the eigenvectors obtained from the gene expression data and the give(n) pathway databases. The given eigenvectors were used to fit a logistic model to identify an intercept and the coefficients which were used to calculate a gene log_odds score. This score is translated to z-scores by using a permutation strategy. The gene-components matrix is permutated to create a null-distribution. The fitted intercept and coeffects were used to calculate the gne log_odds scores corresponded to the null distribution and the mean and standard deviation (std) of these values were calculated. The mean and std are used to translate the gene log_odds score to z_scores.

The final gene z-scores were used to calculation the performance of the model by calculating the Area under the ROC-curve (AUC value) and the corresponded p-value by using the Mann-Whitney rank test. The p_values are also Bonferroni corrected.

The input files of the second step of the pipeline consist of the eigenvector matrix from created by the decomposition step. The input can be both the .txt file or the pickle file. The second input file contains the pathway matrix file which indications which gene is linked to which pathway. The file must contain a tab separated matrix where the rows correspond to the pathways and the columns to the genes. The values must be Booleans indicating is a gene is linked (value is 1) or not linked (value is 0) to a pathway.

The model is by default fitted on the complete gene list, but a gene sub selection can be given to fit the model only on a subset of genes. The subset can be given as path to a file with the gene names, where every row corresponds to one gene.

The output of the pipeline contains one file with the gene z-scores for every pathway (gene_pathway_scores) where the rows correspond to every pathway and the columns to every gene. The second output file contains the prediction AUC values and the corresponded p value and Bonferroni corrected p value. The output is by default exported to both txt and pickle format. The .txt file is gzip compressed. If some files are not necessary, they can be disabled by using the output parameters.

The permutated version of the input matrix will be saved as pandas pickle file on the same location of the input file with the suffix “_permutation_matrix”. If this matrix is already present, the script will read this file instead of creating a new permutated matrix (except when the option –force is set)

The process of calculation of z-scores for one pathway can be take a serious amount of time, especially for pathways with a lot of genes. To reduce the overall processing time, the script use by default multiple CPU cores of the machine. The maximal number of available cores are used, but this can be reduced by setting a maximal number of cores through a parameter. The script has also build-in support for multi node processing, for cluster usage. This can be easily setup with a job manager like Slurm. An example job script is included.

Both scripts use input data frames to perform the analysis. Every input data frame can be given as Pandas pickled data frame or as a tab separated txt file. The txt input file is automatically decompressed if the file paths ends withs ‘.gz’, ‘.bz2’, ‘.zip’, or ‘.xz’. An cached pandas pickle file are created on the same location as the input txt file to speed up the reading process of the input file in further analysis. These cached files automatically created on the same location and same name as the input file with the suffix “_cashed.pickle”.

Both scripts create different output files as described in the above. The files are created in the given output directory, which is also created if the directory does not exit. Output matrixes are by default created as pandas pickled dataframe and tab separated text file. The txt matrix files are by default gz compressed. The compressing can be disabled (parameter: --output_disable_gzip), even as the creation of the txt file or the pickle file (option: --output_disable_txt and --output_disable_pickle)

This program is developed in Python 3.7, performance on other versions is not guaranteed.

The program requires the following packages to be installed:

• pandas (v0.25.1; BSD 3-Clause License)
• numpy (v1.17.2; BSD License)
• sklearn (v0.22.2; BSD 3-Clause License)
• matplotlib (v3.1.0; PSF License)
• seaborn (v0.9.0; BSD 3-Clause License)

See 'Installing' on how to install these packages.

Installing and running the program can be done by executing the following steps:

Step 1: acquire the source files Download or clone the GeneNetworkAnalysis directory from the systemsgenetics repository.

Step 2: installing the required packages

1. Install the required packages. The python packager manager (pip) and the requirements file can be used to install all the nessesary packages:

   pip install -r requirements.txt
   

2. Wait for the command to finish

De repository contains different scripts to create and analyse gene expression networks. Most of the scripts contain an argument parser and help function.

The decomposition step is implemented in the “gene_network_decomposition.py” script. The help and parameter instructions can be found by using the command:

```  
python3 gene_network_decomposition.py --help
```  

The decomposition step is implemented in the “gene_pathway_score_prediction.py” script. The help and parameter instructions can be found by using the command:

```  
python3 gene_pathway_score_prediction.py --help
```  

Henry Wiersma (1)

1. Genetics Department of the University Medical Center Groningen (UMCG), The Netherlands.