Input object

The first object is an input object for the disease module inference methods, with class MODifieR_input. This object can be used for all disease module inference methods.

The MODifieR_input class has two subclasses, either MicroArray or RNA-seq depending on the input data type.

The input objects contain all necessary input data for the methods. Below is a table showing what methods use an expression matrix, Differentially Expressed Genes (DEGs) or both.

Table 1: Input types taken by inference methods

Module object

The second object, with class MODifieR_module is an object that is produced by the disease module inference methods. They also have a subclassses according to the method and input type (Microarray or RNA-seq) that produced them.

Like the input object, the module objects are also named lists. Module objects differ between each function that generates them (hence the subclass) but have certain features in common:

• The module object will always contain an element (character vector) called module_genes which are the genes that make up the final module

• The module object will always contain a named list called settings which holds the parameters that have been used when generating the object

Using microarray data

To create an input object with microarray data, the function create_input_microarray can be used. The function takes several inputs:

• expression_matrix A matrix with normalized expression values. The rows are probes and the columns are samples.

• annotation_table An two columns dataframe providing annotation for the probes. The first column should contain the probe ID, and the second column its corresponding annotation.

• group1_indici A vector that should have the indici in the expression matrix for the samples belonging to group 1.

• group2_indici Same, but for group 2

• group1_label Label to identify group 1, for example “Control” or “Patient”

• group2_label Same, but for group 2

• expression A Boolean value indicating if expression matrix data should be calculated.

• differential expression A Boolean value indicating if DEGs data should be calculated.

• use_adjusted use adjusted p value for differential expression analysis?

• method When collapsing the probes, the collapseRows() method from WGCNA is used. This argument determines the collapsing method to use. Please see the WGCNA documentation for more information. Options are:
  • MaxMean Highest mean value
  • MinMean Lowest mean value
  • MaxRowVariance Highest variance
  • absMaxMean Highest mean absolute value
  • absMinMean Lowest mean absolute value
  • ME Eigenrow (first principal component)

Using RNA-seq data

Creating an RNA-seq input object is quite similar to creating a microarray input object. The function is create_input_rnaseq and the arguments are:

• count matrix Matrix containing raw RNA-seq counts

• group1_indici A vector that should have the indici in the expression matrix for the samples belonging to group 1.

• group2_indici Same, but for group 2

• group1_label Label to identify group 1, for example “Control” or “Patient”

• group2_label Same, but for group 2

• expression A Boolean value indicating if expression matrix data should be calculated.

• normalize_quantiles, Normalize quantiles for WGCNA-based methods?

Example dataset microarray

Example data has been included in the package. There is an example microarray expression matrix from Gene Expression Omnibus (GEO) with accession number GSE4588. This data set has 16 columns, so 16 samples in total with 9 systemic lupus erythematosus (SLE) patients and 7 healthy controls. The matrix has 52307 rows, meaning 52307 probes. The indici for the healthy controls are 1 to 9, and 10 to 16 make up the SLE patients. The expression matrix can be accessed using expression_matrix

We can inspect it using the head function:

head(expression_matrix)
                GSM101962 GSM101963 GSM101964 GSM101965 GSM101966
AFFX-BioB-5_at      442.0     442.6     230.8     275.2     540.7
AFFX-BioB-M_at      620.7     565.5     256.9     299.8     693.7
AFFX-BioB-3_at      404.7     376.1     175.4     168.2     403.9
AFFX-BioC-5_at     1298.3    1201.6     742.4     781.2    1485.3
AFFX-BioC-3_at     1647.5    1566.8     966.2    1059.7    1743.9
AFFX-BioDn-5_at    2952.5    2819.6    2107.7    2444.5    3265.4
                GSM101967 GSM101968 GSM101969 GSM101970 GSM102710
AFFX-BioB-5_at      432.7     476.4     953.9    1397.2     197.0
AFFX-BioB-M_at      658.5     605.7    1313.7    1863.5     256.8
AFFX-BioB-3_at      440.5     379.5     889.9    1277.7     159.4
AFFX-BioC-5_at     1406.0    1301.2    2368.0    3752.2     674.3
AFFX-BioC-3_at     1720.4    1598.2    2727.0    3771.3     867.4
AFFX-BioDn-5_at    3150.1    2911.8    5480.8    7806.6    2055.5
                GSM102711 GSM102712 GSM102713 GSM102714 GSM102715
AFFX-BioB-5_at      184.2     415.8     499.4    1249.3     380.3
AFFX-BioB-M_at      237.2     560.7     729.1    1796.1     557.1
AFFX-BioB-3_at      190.6     359.2     445.2    1194.3     315.8
AFFX-BioC-5_at      573.9    1368.8    1485.1    3886.1    1300.5
AFFX-BioC-3_at      841.0    1476.6    1891.9    4485.8    1706.6
AFFX-BioDn-5_at    1774.7    2927.9    3254.2    8357.6    2781.7
                GSM102716
AFFX-BioB-5_at      339.5
AFFX-BioB-M_at      576.0
AFFX-BioB-3_at      358.9
AFFX-BioC-5_at     1262.9
AFFX-BioC-3_at     1466.3
AFFX-BioDn-5_at    3058.1

The annotation for the probes is provided in a data.frame accessible through probe_annotation. The first column contains the probe ids where the second column holds the corresponding annotation, in this case ENTREZ gene identifiers.

Again using head:

head(probe_annotation)
          PROBEID ENTREZID
1  AFFX-BioB-5_at     <NA>
2  AFFX-BioB-M_at     <NA>
3  AFFX-BioB-3_at     <NA>
4  AFFX-BioC-5_at     <NA>
5  AFFX-BioC-3_at     <NA>
6 AFFX-BioDn-5_at     <NA>

Now we can put these things together to create an input object:

MODifieR_input_microarray <- create_input_microarray(expression_matrix = expression_matrix,
                               annotation_table = probe_annotation,
                               group1_indici = 1:9,
                               group2_indici = 10:16,
                               group1_label = "Control",
                               group2_label = "Patient")

When the function has finished, it might take a few minutes, we can inspect the object. Let’s start with the class:

[1] "MODifieR_input" "Expression"     "MicroArray"    

Now we can inspect the DEGs:

       gene      pvalue
1         1 0.976932866
2        10 0.347847172
3       100 0.006504456
4      1000 0.329104066
5     10000 0.115466316
6 100009676 0.433955896

The limma_probe_table is the table generated by limma with annotation for the probe added:

      logFC   AveExpr        t      P.Value    adj.P.Val          B
1 4455.8984 2263.4000 16.21819 6.973635e-11 2.143054e-06 -0.4889325
2  883.5587  710.8125 16.03447 8.194137e-11 2.143054e-06 -0.4961970
3 1995.8476 1046.2500 12.89384 1.721050e-09 2.255049e-05 -0.6651916
4  845.2206  549.2062 12.89198 1.724472e-09 2.255049e-05 -0.6653242
5  339.1746  240.1000 11.96197 4.802206e-09 5.023779e-05 -0.7382467
6 1924.1302 1276.3625 10.46105 2.902571e-08 2.280640e-04 -0.8907192
      PROBEID IDENTIFIER
1   204439_at      10964
2   213294_at       5610
3 214453_s_at      10561
4   227609_at      94240
5   219352_at      55008
6   202086_at       4599

Finally, we can take a look at the groups indici and it’s names by looking at the group_indici component of the object:

$Control
[1] 1 2 3 4 5 6 7 8 9

$Patient
[1] 10 11 12 13 14 15 16

Now we have set up an microarray input object to use in the disease inference methods.

Example dataset RNA-seq

There is an example RNA-seq expression matrix from Gene Expression Omnibus (GEO) with accession number GSE123496. This data set has 10 columns, so 10 samples taken from the frontal cortex with in total 5 Multiple Sclerosis (MS) patients and 5 healthy controls. The matrix has 16569 rows, meaning the same amount of genes. The indici for the patients are 1 to 5, and 6 to 10 make up the controls. The count matrix can be accessed using count_matrix

We can inspect it using the head function:

head(count_matrix)
      MS32CF MS38CF MS44CF MS50CF MS56CF NL74CF NL80CF NL86CF NL62CF
7105     251     90    209    336    358    191    189    172    124
8813     235    121    404    428    330    594    242    435     89
57147    224     64    280    437    387    545    259    478     86
55732    109     37    145    203    189    214    114    197     37
2268     336    237    227    345    733    262    482    130    457
3075     867   1181    951   1662   1825    769   1465    285   1235
      NL68CF
7105     335
8813     324
57147    408
55732    203
2268     423
3075     579

To create the input object:

MODifieR_input_rna <- create_input_rnaseq(count_matrix = count_matrix, 
                                          group1_indici = 1:5, 
                                          group2_indici = 6:10, 
                                          group1_label = "Patient", 
                                          group2_label = "Control",
                                          use_adjusted = F)

This dataset does not contain any significant genes after FDR correction, so unadjusted p values are used.

When the function has finished, we can inspect the object. Let’s start with the class:

[1] "MODifieR_input" "Expression"     "RNA-seq"       

Now we can inspect the DEGs:

       gene       pvalue
1      4241 8.162947e-05
2 105378049 1.623881e-04
3 100423062 2.969720e-04
4    155185 6.942163e-04
5       284 9.156738e-04
6     55791 9.238056e-04

The edgeR_deg_table is the table generated by edgeR:

              logFC    logCPM        F       PValue       FDR
4241       1.071991 4.0142387 36.67993 8.162947e-05 0.9579268
105378049 -2.247071 0.9811993 31.20710 1.623881e-04 0.9579268
100423062 -5.161158 0.1725241 26.93895 2.969720e-04 0.9579268
155185     1.011908 3.4099991 21.68029 6.942163e-04 0.9579268
284       -1.264745 1.8260861 20.13550 9.156738e-04 0.9579268
55791     -0.897363 3.8180308 20.08746 9.238056e-04 0.9579268

Finally, we can take a look at the groups indici and it’s names by looking at the group_indici component of the object:

$Patient
[1] 1 2 3 4 5

$Control
[1]  6  7  8  9 10

Now we have set up an RNA-seq input object to use in the disease inference methods.

Creating a custom input object

There is also the possibility to create a custom input object. This can be useful in cases where there are already DEGs, an annotated expression or both available. In order to create such an object, the functions create_custom_microarray_input_object or create_input_rna_input_object can be used.

Below are example on how to use the functions. Some of the input data is optional, for example limma_probe_table/edgeR, expression_matrix/count_matrix and annotation_table are not used in any disease module inference method. Instead, they can optionally be used to update or recalculate an input object.

generic_input_microarray_object <- create_custom_input_object(diff_genes = NULL, 
                                                              limma_probe_table = NULL,
                                                              annotated_exprs_matrix = NULL, 
                                                              expression_matrix = NULL, 
                                                              annotation_table = NULL, 
                                                              group1_indici = NULL,
                                                              group2_indici = NULL, 
                                                              group1_label = NULL, 
                                                              group2_label = NULL)

generic_input_rna_object <- create_custom_rna_input_object(diff_genes = NULL, 
                                                           edgeR_deg_table = NULL, 
                                                           annotated_exprs_matrix = NULL, 
                                                           count_matrix = NULL, 
                                                           group1_indici = NULL, 
                                                           group2_indici = NULL, 
                                                           group1_label = NULL, 
                                                           group2_label = NULL, 
                                                           settings = NULL)

Please note that the specific data in the input object that is required differs per inference method, the reader is referred to table 1 in subsection ‘Input object’ to see what inference methods needs what input data fields. For example, to use the DIAMoND method, only the diff_genes input would suffice to run the analysis.

If the inference method requires annotated_exprs_matrix, group indici and labels are also required.

Clique Sum

There are two different versions of the Clique Sum algorithm implemented, clique_sum_exact and clique_sum_permutation, that differ in the method to determine if a clique is significantly enriched with DEGs. Obtaining cliques is done in a separate function, build_clique_db, that creates an SQLite database for the cliques. This is because obtaining cliques is a computationally intensive process, may clog the memory when using a large network and is static for each network.

Both versions take in the file name of this database instead of a PPi network.

build_clique_db(ppi_network = ppi_network, 
                db_folder = "/folder/name/", 
                db_name = "example")

clique_sum_exact_module <- clique_sum_exact(MODifieR_input = MODifieR_input_microarray, 
                                            db = "folder/name/example.sqlite")

clique_sum_permutation_module <- clique_sum_permutation(MODifieR_input = MODifieR_input_microarray, 
                                                        db = "folder/name/example.sqlite")

Correlation Clique

The correlation clique is a clique-based algorithm using consensus clustering. The algorithm starts with calculating a correlation score between each interaction in the PPi network. The correlation score is obtained by subtracting the Pearson correlation p-value: \(correlation score = 1 - Pearson\: p\: value\)

Subsequently, the correlation score is multiplied by the correlation confidence. $edge score = $

When the edge scores are calculated the iterative part of the algorithm commences: All edge scores are scaled to so that the mean of all edge scores is equal to the fraction_of_iteration parameter. Next, the newly scaled scores are compared to random variables from the uniform distribution between (0,1) and only interactions where the edge score is higher than the random variable are used to construct a new PPi network. Because the scores are scaled to fraction_of_iteration, this is the average fraction of interactions used for this new network. Then, maximal cliques are inferred from this new network. The cliques are tested for significant enrichment of DEGs by Fisher’s exact test and the union of significant cliques is the disease module for this iteration. The final disease module will consist of genes that have been present in at least frequency cutoff iterations

The function can be called using correlation_clique()

In the example code, the number of iterations has been reduced to 10 (from 50) to speed up the analysis

correlation_clique_module <- correlation_clique(MODifieR_input = MODifieR_input_microarray, 
                                                ppi_network = ppi_network, 
                                                iteration = 10)

correlation_clique_module <- correlation_adjust_cutoff(frequency_cutoff = 0.2, 
                                                       correlation_module = correlation_clique_module)

correlation_clique_module <- correlation_set_module_size(size = 10, 
                                                         correlation_module = correlation_clique_module)

DIAMoND

A slightly modified version of the original DIAMOnD python script (Ghiassian, Menche, and Barabási 2015) is called from within R. The only change to the original algorithm is the option to include the seed genes to the module.

In MODifieR, the function can be called using diamond()

diamond_module <- diamond(MODifieR_input = MODifieR_input_microarray, 
                          ppi_network = ppi_network)

diamond_module <- diamond_add_seed_genes(diamond_module = diamond_module)

diamond_module <- diamond_remove_seed_genes(diamond_module = diamond_module)

DiffCoEx

DiffCoEx is a method for identifying correlation pattern changes, which builds on the commonly used Weighted Gene Co expression Network Analysis (WGCNA) framework for co expression analysis. (Tesson, Breitling, and Jansen 2010)

In MODifieR, the function can be called using diffcoex()

diffcoex_module <- diffcoex(MODifieR_input = MODifieR_input_microarray)

diffcoex_module_list <- diffcoex_split_module_by_color(diffcoex_module = diffcoex_module)

MCODE

A clique based algorithm to identify disease modules from differentially expressed genes originally by Bader et al. (Bader and Hogue 2003)

Much of the code and documentation has been taken from the now defunct package “ProNet”

In MODifieR, the function can be called using mod_mcode()

The MCODE module is the union of all submodules with a score above module_score.

mcode_module <- mod_mcode(MODifieR_input = MODifieR_input_microarray, 
                               ppi_network = ppi_network)

mcode_module <- mcode_adjust_score(module_cutoff = 12, mcode_module = mcode_module)

mcode_module_list <- mcode_split_modules(module_cutoff = 3, mcode_module)

MODA

This implementation follows a workflow as described in the MODA vignette. (Li et al. 2016) First, two separate networks are constructed, a background network containing expression data from all samples and a condition specific network consisting of all samples minus the condition specific samples. Then, hierarchical clustering is performed and cutting height estimated from either maximal average density or modularity

Condition specific co-expression modules are then extracted using the Jaccard index and specificTheta.

The final module will consist of the co-expression module that has the minimal Jaccard index complemented by co-expression modules that have a Jaccard index below this minimal + specificTheta. The group of interest signifies which group has the trait of interest. In this case that is group 2, the SLE patients.

In MODifieR, the function can be called using moda()

moda_module <- moda(MODifieR_input = MODifieR_input_microarray, 
                    group_of_interest = 2)

moda_module <- moda_change_specific_threshold(moda_module = moda_module, 
                                              specificTheta = 0.05)

Module Discoverer

This is an implementation of the single seed Module Discoverer algorithm. The code has been adapted from the original code by Vlaic et al. (Vlaic et al. 2018) For details, please see the paper referenced.

In MODifieR, the function can be called using modulediscoverer()

module_discoverer_module <- modulediscoverer(MODifieR_input = MODifieR_input_microarray, 
                                             ppi_network = ppi_network, 
                                             permutations = 1000, 
                                             repeats = 3)

WGCNA trait-based

Wgcna is an implementation of WGCNA (Langfelder and Horvath 2008), (Langfelder and Horvath 2012) that associates co-expression modules (denoted by color) to a trait. Co-expression modules with an adjusted p-value < pval_cutoff will make up the final disease module.

The algorithm infers co-expression modules from combined expression data set from both group1 and group2. Co-expression modules are then correlated to trait (group 1 ~ group 2). The group of interest signifies which group has the trait of interest. In this case that is group 2, the SLE patients.

In MODifieR, the function can be called using wgcna()

wgnca_module <- wgcna(MODifieR_input = MODifieR_input_microarray,
                      group_of_interest = 2)

wgnca_module <- wgcna_adjust_significance(pval_cutoff = 0.1, 
                                          wgcna_module = wgnca_module,
                                          use_unadjusted = FALSE)

wgnca_module <- wgcna_get_module_genes_by_sign(wgcna_module = wgnca_module, 
                                               mode = "p")

wgnca_module <- wgcna_split_module_by_color(wgcna_module = wgcna_module)

wgnca_module <- wgcna_set_module_size(size = 200,
                                      wgcna_module = wgnca_module)

Count based

A count based consensus module can be derived with the function create_module_set(). It will return a module object of class MODifieR_module, same as the individual methods. The subclass is Module_set.

The method derives a consensus module from a set of modules is by tabulating the presence of each gene in the individual methods. The next step is to set a minimum for this presence; for example, a gene has to be present in at least 3 of the individual modules.

To start deriving a consensus module, it is convenient to wrap all the previously acquired modules into a list:

module_list <- list(clique_sum_permutation_module, 
                    correlation_clique_module, 
                    diamond_module, 
                    diffcoex_module, 
                    mcode_module, 
                    moda_module, 
                    module_discoverer_module, 
                    wgcna_module)

Now the maximal number of modules that share at least one gene can be retrieved:

max_frequency <- get_max_frequency(module_list)

The maximum in this particular example turns out to be 5. This value can be used to retrieve all consensus modules, from a minimum presence of 1 (union of modules) up to 5:

consensus_modules <- lapply(X = 1:max_frequency, FUN = create_module_set, 
                            module_list = module_list)

A list of length 5 with Module_set objects is returned. In this tutorial the fourth element in the list, so the consensus module where genes are present in at least 4 individual methods will be inspected.

consensus_module4 <- consensus_modules[[4]]

Like other module objects, it has a vector containing the module genes and a list containing the settings. In addition, it has the following elements:

• module_gene_list A named list containing the module genes from the original modules

• gene_frequencyTable containing all the genes present in the modules and their frequency

• method_by_gene A named list where the elements are the modules the genes have been found in and the names are the gene names

• gene_by_method A table containing the gene frequencies by combination of methods

The resulting module object can now be used for further analysis.

In this example the module genes from consensus_module4 are written to a file:

write.table(x = consensus_module4$module_genes, 
            file = "consensus4_genes.txt", 
            quote = F, 
            row.names = F, 
            col.names = F)

And uploaded into String. This resulted in Figure 1 and Table 2.

Figure 1: Plot from https://string-db.org using the module genes from consensus_module4 as input

Table 2: Excerpt from Molecular Function GO term table on the result page of String

Specific betweenness (S2B)

An S2B consensus module can be derived with the function S2B(). It will return a module object of class MODifieR_module, same as the individual methods. The subclass is S2B.

The method derives a consensus module between two modules. Originally conceived by (Garcia-Vaquero et al. 2018) to find an overlap between two gene sets from related diseases, it can also be applied to find an overlap between gene sets tht have been inferred with different methods but are from the same disease.

The method works by prioritizing genes frequently and specifically present in shortest paths linking two disease modules. For a detailed description of the algorithm, see (Garcia-Vaquero et al. 2018).

The function is called using S2B()

In the example below, the clique_sum_permutation and MODA modules are used, and the number of iterations (nrep and nrep2) is set to 50 instead of 100 to shorten computation time.

 S2B_module <- S2B(MODifieR_module1 = clique_sum_permutation_module, 
                   MODifieR_module2 = moda_module, 
                   ppi_network = ppi_network, 
                   nrep = 50, 
                   nrep2 = 50)

As with all module objects, the S2B object has a vector called module_genes for the genes that make up the module and also a list containing the settings that have been used to generate the object.

In addition it has two related fields, input_class1 and input_class2, that show the class of the MODifieR_input objects that have been used as input. In the case of the tutorial they are Clique_Sum_permutation and MODA

Plotting

A very crude plotting function has been implemented to plot MODifieR module objects. To do this, the regular function plot() can be used. It takes an additional argument, a PPi network to overlay the module genes on. The size of the nodes on the plot is dependent on degree, nodes with a high degree are larger. In general, the plots are not great but can be used in an exploratory fashion when inferring modules.

Plotting consensus_module4 again with the example PPi included in MODifieRDev:

plot(consensus_module4, ppi_network = ppi_network)

Figure 2: Using the plot function in MODifieRDev with consensus_module4 and the included PPi network as inputs

Figure 3: Using the plot function in MODifieRDev with the diamond module and the included PPi network as inputs

Exporting

MODifieR objects can be exported to both csv and an excel workbook.

export(MODifieR_object = consensus_modules[[4]], folder = "modules/consensus_module_4")

export(MODifieR_object = MODifieR_input_microarray, folder = "inputs/microarray")

export_xlsx(MODifieR_object = diamond_module, folder = "modules/diamond_module")

export_xlsx(MODifieR_object = MODifieR_input_microarray, folder = "inputs/microarray")