1.1 Load environment

Libraries used to create and generate this report:

• R : R version 4.3.3 (2024-02-29)
• rmarkdown : 2.21
• knitr : 1.42
• rmdformats : 1.0.4
• bookdown : 0.34
• kableExtra : 1.3.4

1.2 Load libraries

Libraries used to analyse data:

library("SNFtool")
library("pheatmap")
library("igraph")

• SNFtool: 2.3.1
• pheatmap: 1.0.12
• igraph: 1.4.2

Libraries used to load data:

library("MOFAdata")
library("data.table")
library("mixOmics")
library("WallomicsData")

• MOFAdata: 1.16.1
• data.table: 1.14.8
• mixOmics: 6.24.0
• WallomicsData: 1.0

1.3 Visualization using Cytoscape

Cytoscape is used for visualization. Figures were generated using the Cytoscape v3.9.1 and several Cytoscape apps:

• yFiles Layout Algorithm: 1.1.3
• LegendCreator: 1.1.6

3.1 Metagenomic dataset from Tara Ocean project

To retrieve data: files are available in /shared/projects/tp_etbii_2024_165650/Networks/TaraOcean_mibiomics directory path in the IFB server.

• dataset:

  • TARAoceans_proNOGS.cvs
  • TARAoceans_proPhylo.csv

• metadata:

  • TARAoceans_metadata.csv

Samples come from eight oceans around the world (SPO: South Pacific Ocean, NAO: North Atlantic Ocean, IO: Indian Ocean, RS: Red Sea, MS: Mediterranean Sea, NPO: North Pacific Ocean, SO: Southern Ocean, SAO: South Atlantic Ocean).

Samples can come from different layers with different temperatures:

• SRF: Surface Water Layer (0-5 meters)
• DCM: Deep Chlorophyll Maximum (peak of chlorophyll, 0-600 meters)
• MIX: Subsurface epipelagic Mixed Layer
• MES(O): Mesopelagic zone (from 500/1000 meters)

In a previous analysis (Sunagawa et al., 2015), they identified a stratification mostly driven by the temperature rather than geography or other environmental factors.

We have two types of data:

• orthologous genes: the relative abundance of groups of orthologous genes (OGs)
• phylogenetic profil: counts of S16 rRNA

Does an integrative analysis of these two data types retrieve the stratification driver by the layers? Does it also find a geographical clustering?

Data are coming from: MiBiOmics gitlab.

3.2 Breast cancer dataset from The Cancer Genome Atlas

To retrieve data:

• dataset: using data("breast.TCGA") from the mixOmics R package

  • breast.TCGA$data.train$mirna
  • breast.TCGA$data.train$mrna
  • breast.TCGA$data.train$protein

• metadata: using data("breast.TCGA") from the mixOmics R package

  • breast.TCGA$data.train$subtype

Human breast cancer is a heterogeneous disease. Breast tumors can be classified into several subtypes (PAM50 classification), according to the mRNA expression level (Sorlie et al., 2001). In this dataset, we have three subtypes:

• Basal: considered more aggressive than LumA
• Her2: tend to grow faster than LumA and can have a worse prognosis, but are usually successfully treated
• LumA: tend to grow more slowly than other cancers, be lower grade, and have a good prognosis

We have three types of data:

• mRNA: mRNA expression level
• miRNA: microRNA expression level
• protein: protein abundance

Does an integrative analysis of these three data types retrieve the classification of the breast cancer? Or find another classification?

Data are coming from the mixOmics R package. The full data can be downloaded here.

3.3 Chronic Lymphocytic Leukaemia (CLL) dataset

To retrieve data:

• dataset: using data("CLL_data") from the MOFAdata R package

  • CLL_data_t$Drugs
  • CLL_data_t$Methylation
  • CLL_data_t$mRNA
  • CLL_data_t$Mutations

• metadata: file is available in /shared/projects/tp_etbii_2024_165650/Networks/CLL directory path in the IFB server.

  • sample_metadata.txt

The Chronic Lymphocytic Leukaemia (CLL) is type of blood and bone marrow cancer. The full data are explained in Dietrich et al., 2018 and available here.

We have four types of data:

• mRNA: transcriptom expression level
• methylation: DNA methylation assays
• drug: drug response measurements
• mutation: sommatic mutation status

3.4 Tomato plant dataset

To retrieve data: files are available in /shared/projects/tp_etbii_2024_165650/Networks/Tomato directory path in the IFB server.

• dataset:

  • mrna.tsv
  • prots.tsv

• metadata:

  • samples_metadata.tsv

In order to study the protein turnover in developing tomato fruit (Solanum lycopersicum) in Belouah et al., two omics data types were collected:

• transcript data: gene abundance
• protein data: protein abundance

Each data type was collected in nine different developmental stages: GR1, GR2, GR3, GR4, GR5, GR6, GR7, GR8 and GR9. For each developmental stages, we have three replicates.

Does an integrative analysis of these data types retrieve the different developmental stages?

Data are coming from Belouah et al., 2019.

3.5 Arabidopsis thaliana

To retrieve data:

• dataset: using data("X") from the WallomicsData R package

  • Phenomics_Rosettes
  • Transcriptomics_Rosettes_CW
  • Proteomics_Rosettes_CW
  • …

• metadata: using data("X") from the WallomicsData R package also

  • Altitude_Cluster
  • Ecotype
  • …

In order to study the cell wall plasticity of Arabidopsis thaliana plants, exposed to different temperature growth conditions, four omics data types were collected :

• phenomics data: 5 phenotype variables measured
• metabolomics: identification and quantification of the seven cell wall monosaccharides
• cell wall proteomics: counts of cell wall proteins
• transcriptomics: counts of transcripts

Each data type was collected for

• 2 organs, rosettes and stems with three replicates.
• 2 temperatures, 15°C an d 22°C
• 5 ecotypes, Roch, Grip, Hern and Hosp growing at different altitudes in the Pyrénées and Col from low altitude in Poland
• 3 genetic clusters

To see all the available data, have a look to the manual.

Data are explained in Duruflé et al., 2019, 2020 and 2021.

4.1 In summary

The preprocessing step is the most important part of the analysis. Data need to be prepared correctly in order to extract relevant information and produce a correct and pertinent interpretation of the results.

Data preprocessing could be summarized by four main steps:

1. Prepare the data (e.g. remove outliers, correct bacth effect etc…)
2. Remove and/or impute missing data
   • remove features/samples if more than 20% of missing data
   • impute missing data (e.g. using K-nearest neighbor method (KNN))
3. Normalize the data according to the data type
4. Scale (mean = 0 and standard deviation = 1)

Figure 4.1: Distribution examples expected after preprocessing

The data must conform to a specific matrix shape:

• samples (e.g. samples, organisms …) in rows
• features (e.g. genes, proteins …) in columns

4.2 Read and prepare the data

For this tutorial, we assume that the data have been already prepared: outliers are already removed and there is no batch effect.

tara_nog <- read.table(file = "../00_Data/TaraOcean_mibiomics/TARAoceans_proNOGS.csv", sep = ",", head = TRUE, row.names = 1)
tara_nog[c(1:5), c(1:5)]

##              NOG317682    NOG135470 NOG85325    NOG285859    NOG147792
## TARA_109_SRF         0 2.390962e-05        0 4.663604e-08 1.800215e-07
## TARA_149_MES         0 4.339824e-06        0 5.182915e-07 4.190123e-06
## TARA_110_MES         0 1.348252e-05        0 6.000043e-07 2.218342e-07
## TARA_102_MES         0 6.380711e-06        0 3.816016e-07 0.000000e+00
## TARA_142_SRF         0 9.484144e-06        0 6.437103e-08 1.132431e-06

4.3 Missing data

In the SNF paper, authors recommend to filter out samples with more than 20% of missing data in a certain data type. They also recommend to filter out features with more than 20% of missing data across samples. Then, they impute the remaining missing data using K nearest neighbors (KNN) imputation.

In this tutorial, we decide to remove samples with at least one missing data. To remove samples with missing data, we propose the following NARemoving() function. Input parameters are:

• data: the data type
• margin: a vector giving the subscripts which the function will be applied over (e.g. 1 indicates rows and 2 indicates columns)
• threshold: threshold above which samples/features are deleted

NARemoving <- function(data, margin, threshold){
    #' NA removing
    #'
    #' Calculate percentage of na
    #' Remove na from rows (margin = 1) or column (margin = 2)
    #' 
    #' @param data data.frame.
    #' @param margin int. 1 = row and 2 = column
    #' @param threshold int. Number of missing data accepted
    #'  
    #' @return Return data.frame with a specific number of na by row/column

  data_na <- apply(data, MARGIN = margin, FUN = function(v){sum(is.na(v)) / length(v) * 100})
  # print(table(data_na))
  toRemove <- split(names(data_na[data_na > threshold]), " ")[[1]]
  if(margin == 1){
    data_withoutNa <- data[!(row.names(data) %in% toRemove),]
    print(paste0("Remove ", as.character(length(toRemove)), " samples."))
  }
  if(margin == 2){
    data_withoutNa <- data[,!(colnames(data) %in% toRemove)]
    print(paste0("Remove ", as.character(length(toRemove)), " features"))
  }
  return(data_withoutNa)
}

For instance, the CLL drug data contain missing data (sample H024).

CLL_data_t$Drugs[c(1:5), c(1:5)]

##         D_001_1    D_001_2   D_001_3   D_001_4   D_001_5
## H045 0.02363938 0.04623274 0.3187471 0.8237027 0.8962777
## H109 0.07359900 0.10623002 0.2732891 0.7171379 0.8850003
## H024         NA         NA        NA        NA        NA
## H056 0.05813930 0.09022028 0.2322145 0.7225736 0.7957497
## H079 0.02042077 0.04750543 0.3638962 0.8073907 0.8794886

To remove samples with missing data in drug data, we use the following command line:

• data = CLL_data_t$Drugs: remove the samples with missing data in the drug data
• margin = 1: samples are in rows, so we want to apply the function on the rows
• threshold = 0: we remove samples with at least one missing data

CLL_drug <- NARemoving(data = CLL_data_t$Drugs, margin = 1, threshold = 0)

## [1] "Remove 16 samples."

The H024 sample is not anymore in the CLL drug data.

CLL_drug[c(1:5), c(1:5)]

##         D_001_1    D_001_2   D_001_3   D_001_4   D_001_5
## H045 0.02363938 0.04623274 0.3187471 0.8237027 0.8962777
## H109 0.07359900 0.10623002 0.2732891 0.7171379 0.8850003
## H056 0.05813930 0.09022028 0.2322145 0.7225736 0.7957497
## H079 0.02042077 0.04750543 0.3638962 0.8073907 0.8794886
## H164 0.02962725 0.08054628 0.4725991 0.8179143 0.8927961

We repeat this step for each data type. Then, we filter out samples that are not present in all data type.

sampleNames <- Reduce(intersect, list(rownames(CLL_drug), rownames(CLL_mrna)))
CLL_drug <- CLL_drug[rownames(CLL_drug) %in% sampleNames,]
CLL_mrna <- CLL_mrna[rownames(CLL_mrna) %in% sampleNames,]

These functions could help you: is.na(), lapply(), table() and view().

4.4 Normalization and scaling

4.4.2 Scaling

Each feature (column) needs to have the mean equals to zero and the standard deviation equals to one. For that, the SNF package provides the function standardNormalization().

In the Figure 4.2, you can see the data distribution for the breast cancer miRNA data before (on the left) and after (on the right) scaling. After scaling, the data distribution should be normal.

hist(as.matrix(tcga_mirna), nclass = 100, main = "Breast cancer miRNA data", xlab = "values")
hist(as.matrix(tcga_mirna_scaled), nclass = 100, main = "Breast cancer miRNA scaled data", xlab = "values")

Figure 4.2: Breast cancer miRNA data distribution before (left) and after (right) scaling.

Practice:

1. Check if your samples are in rows and your features in columns.
2. Scale your data with standardNormalization().
3. Plot the data distribution before and after scaling for each data type.
4. What can you say about these distribution plots?

These functions could help you: hist(), as.matrix().

5.1 In summary

In this section, we create a sample network for each type of data, based on the similarity between samples. The main steps are (Figure 5.1):

1. Compute the distance between each pair of samples (distance matrix D).
2. Transform distances (from distance matrix D) into weights (into similarity matrix W) using distance with nearest neighbors.
3. Create the corresponding similarity network (G).

Figure 5.1: Similarity network creation overview. Data type 1 is in the first row, represented with the green matrices. Data type 2 is in the last row, represented with the purple matrices.

The similarity network is a sample network. It is defined as a network G with nodes (or vertices) called V and connections (or edges) called E. The connections between samples are weighted. Weights come from the similarity matrix.

We can define the similarity network G like this:

• V: Nodes are samples
• E: Edges are connections between samples
• W: Edge weights are the similarity (or weight) between samples

5.2 Distance calculation

First, we calculate the distance between each pair of samples for each data type using the preprocessed data. The distance method used needs to be adapted to the feature type (e.g. continuous, discrete).

In the SNF paper (Wang et al.), authors suggest to use:

• distance (e.g. euclidean) or correlation (e.g. pearson) for continuous features
• chi-squared distance for discrete features
• agreement based measure for binary features

In the SNF R package, the dist2() function performs a squared Euclidean distances between samples.

These functions could help you: as.matrix(), dim(), nrow() and ncol().

5.3 Similarity calculation

The distance matrix D is then transformed into the similarity matrix W. Distances are converted to weights using the scaled exponential similarity kernel (µ) and average distances between samples and their nearest neighbors (ε):

\[ W = exp(-\frac{D^2}{µ\varepsilon})\] In other words, distances are converted to weights according the distance with the nearest neighbors of each sample pair.

The SNF R package proposes the affinityMatrix() to calculate this similarity matrix. In this case, affinity and similarity are equivalent. This function needs three parameters:

• diff: distance matrix D
• K: number of nearest neighbors (between 10 and 30)
• sigma: hyperparameter or variance (between 0.3 and 0.8)

The K neighbors is used to set the similarities outside of the neighborhood to zero. The sigma parameter allows scaling exponential similarity kernel, which is used to calculated the similarity.

Then, you can visualized the similarity matrix using the pheatmap() function. This function creates an heatmap of samples. By default samples are clustered using hierarchical clustering. For a better visualization, we recommend to:

• remove labels with show_rownames = FALSE and show_colnames = FALSE
• give metadata to the parameter annotation
• apply a log 10 transformation to the similarity matrix (log(x, 10))

6.1 In summary

Previously, we created a similarity matrix W and its corresponding similarity network G for each data type (Figure 6.1).

Figure 6.1: In the previous step, we create a similarity matrix that contains weights for each data type. We also created the corresponding similarity network.

Now, we integrate these similarity matrices (in the Figure 6.1 there are two data types). For that, we use an iterative fusion method. The number of iteration T, needs to be defined.

First, two matrices are created from the similarity matrix of each data type:

• the P matrix: this matrix, also called status matrix, contains normalized weights (come from the similarity matrix W).
• the S matrix: this matrix is the kernel matrix. It contains information about the nearest neighbors.

The P matrix carries the full information about the similarity of each samples to all others.

The S matrix carries the similarity to the K most similar samples for each sample (i.e. topology). The similarities between non-neighboring nodes are set to zero because authors assume that local similarities (high weights) are more reliable than the remote ones.

Then, the P matrix of each data type is iteratively updated with information from P matrices of the other data type, making them more similar at each step. In the Figure 6.2, you have an example with two data type. It’s a bit more complex with more data type (see the paper if you are interested in).

Figure 6.2: Example of the fusion method applied to two data types. Each data type is represented using a color: data type 1 in green and data type 2 in purple.

Finally after T iterations, the P matrices of each data type are merged together to create the final fused similarity matrix, and the corresponding fused similarity network.

Then, you can visualize the fused similarity network using Cytoscape.

6.2 Create the fused similarity matrix

First, we have to define the number of iteration called T. It should be between 10 and 20 (recommended by the authors).

Then, to perform the fusion, SNF R package proposes the SNF() function. You have to provide:

• the list of similarity matrices of each data type
• the number of nearest neighbors K (same as previously)
• the number of iteration T

These functions could help you: list(), length().

6.3 Visualize the fused similarity network

write.table(as_data_frame(W_net), "CLL_W_edgeList.txt", quote = FALSE, col.names = TRUE, row.names = FALSE, sep = "\t")