1 Libraries and environment

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")
  • MOFAdata: 1.16.1
  • data.table: 1.14.8
  • mixOmics: 6.24.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

2 Functions

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 help, we created two functions to calculate these ACC values and choose the best threshold based on the topology:

  • CCCalculation(): this function calculates the Clustering Coefficient (CC) for each node
  • ACCCalculation(); this function averages the CC for a network, in order to obtain the ACC
## CC calculation function
CCCalculation <- function(node, graph){
    #' Clustering Coefficient (CC) calculation
    #'
    #' Calculate the Clustering Coefficient (CC) for each node in a network
    #' 
    #' @param node str.
    #' @param graph igraph. Network object (e.g. the fused network object)
    #'  
    #' @return Return the corresponding CC value
    
  degNode <- degree(graph = graph, v = node, loops = FALSE)
  if(degNode > 1){
    neighborNames <- neighbors(graph = graph, v = node)
    graph_s <- subgraph(graph = graph, vids = neighborNames)
    neighborNb <- sum(degree(graph_s, loops = FALSE))
    CC <- neighborNb / (degNode * (degNode-1))
  }else{CC <- 0}
  return(CC)
}

## ACC calculation function
ACCCalculation <- function(graph){
    #' Average Clustering Coefficient (ACC) calculation
    #'
    #' It average the Clustering Coefficient (CC) of a network
    #' 
    #' @param graph igraph. Network object (e.g. the fused network object)
    #'  
    #' @return Return the corresponding ACC value
    
  nodes <- V(graph)
  ACC <- do.call(sum, lapply(nodes, CCCalculation, graph)) / length(nodes)
  return(ACC)
}

3 Choose your datasets

Choose the dataset on which you want to apply SNF!!


Different datasets are available. Note that each dataset has its specificity and some analysis steps should be adapted.

Four datasets are available: **Metagenomic** dataset from Tara Ocean (image from Sunagawa et al., 2015), **Breast cancer** dataset from TCGA (image from TCGA [website](https://portal.gdc.cancer.gov/)), **CLL** dataset (Dietrich et al., 2018) and **tomato plant** dataset (figure from google image).

Figure 3.1: Four datasets are available: Metagenomic dataset from Tara Ocean (image from Sunagawa et al., 2015), Breast cancer dataset from TCGA (image from TCGA website), CLL dataset (Dietrich et al., 2018) and tomato plant dataset (figure from google image).

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.

4 Tara Ocean dataset

4.1 Input data

The metagenomic dataset from the Tara Ocean project contains 2 data types:

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

Data files are available in /shared/projects/tp_etbii_2024_165650/Networks/TaraOcean_mibiomics directory path in the IFB server.

4.1.1 Load dataset

First, we load the data type from TARAoceans_proNOGS.cvs and TARAoceans_proPhylo.csv files.

4.1.1.1 Orthologous genes (nog) data

The data file contains header (head = TRUE) and the first column contains row names (row.names = 1). Below, the first rows and columns are displayed.

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

Dimensions of the data.

dim(tara_nog)
## [1] 139 638

The nog data contain 139 samples (rows) and 638 features (columns). Data are in the right shape: samples in rows and features in columns.

4.1.1.2 Phylogentic profil (phy) data

The data file contains header (head = TRUE) and the first column contains row names (row.names = 1). Below, the first rows and columns are displayed.

tara_phy <- read.table(file = "../00_Data/TaraOcean_mibiomics/TARAoceans_proPhylo.csv", sep = ",", head = TRUE, row.names = 1)
tara_phy[c(1:5), c(1:3)]
##              EU638706.1.1353 JN537192.1.1500 CAFJ01000195.23.1515
## TARA_109_SRF               0               0                    0
## TARA_149_MES               0               0                    0
## TARA_110_MES               0               4                    0
## TARA_102_MES               0               1                    0
## TARA_142_SRF               0               3                    0

Number of rows in the data.

nrow(tara_phy)
## [1] 139

Number of columns in the data.

ncol(tara_phy)
## [1] 356

The phy data contain 139 samples (rows) and 356 features (columns). Data are in the right shape: samples in rows and features in columns.

4.1.2 Load metadata

The TARAoceans_metadata.csv file contains metadata. It contains header (head = TRUE) and row names (row.names = 1).

tara_metadata <- read.table(file = "../00_Data/TaraOcean_mibiomics/TARAoceans_metadata.csv", sep = ",", head = TRUE, row.names = 1)
head(tara_metadata)
##              ocean depth
## TARA_109_SRF   SPO   SRF
## TARA_149_MES   NAO   MES
## TARA_110_MES   SPO   MES
## TARA_102_MES   SPO   MES
## TARA_142_SRF   NAO   SRF
## TARA_109_DCM   SPO   DCM

The metadata contains two types of information: ocean and depth.

names(tara_metadata)
## [1] "ocean" "depth"

Samples come from 8 oceans.

unique(tara_metadata$ocean)
## [1] "SPO" "NAO" "IO"  "SO"  "SAO" "NPO" "RS"  "MS"

Samples were collected in 4 different depths.

unique(tara_metadata$depth)
## [1] "SRF" "MES" "DCM" "MIX"

4.1.3 Missing data

The data don’t contain missing values. We can go to the following steps.

table(is.na(tara_nog))
## 
## FALSE 
## 88682
table(is.na(tara_phy))
## 
## FALSE 
## 49484

4.1.4 Scaling

We assume that data have been already prepared and normalized.

4.1.4.1 nog data

Nog data are scaled: each column will scaled to have the mean equals to zero and the standard deviation equals to one.

tara_nog_scaled <- standardNormalization(tara_nog)

The following figures are the distribution of the data, before (left) and after (right) scaling. We expected a normal distribution of the data after scaling.

hist(as.matrix(tara_nog), nclass = 100, main = "Orthologous genes - Prepared data", xlab = "values")
hist(tara_nog_scaled, nclass = 100, main = "Orthologous genes - Scaled data", xlab = "values")

Before scaling, data values are almost all closed to zero. After scaling, data values seem to follow something close to a normal distribution. Data values are centered to zero.

4.1.4.2 phy data

Phy data are scaled: each column will scaled to have the mean equals to zero and the standard deviation equals to one.

tara_phy_scaled <- standardNormalization(tara_phy)

The following figures are the distribution of the data, before (left) and after (right) scaling. We expected a normal distribution of the data after scaling.

hist(as.matrix(tara_phy), nclass = 100, main = "Phylogenetic profile - Prepared data", xlab = "values")
hist(tara_phy_scaled, nclass = 100, main = "Phylogenetic profile - Scaled data", xlab = "values")

This kind of data are sparse: there are a lot of zero values. Majority of the values are in the first range on the histogram. After scaling, data values seem to follow something close to a normal distribution. Data values are centered to zero.

4.2 Similarity network

In this part, we create the similarity network for each data type.

4.2.1 Distance calculation

We calculate the Euclidean distance between each pair of samples for each type of data.

tara_nog_dist <- dist2(tara_nog_scaled, tara_nog_scaled)
tara_phy_dist <- dist2(tara_phy_scaled, tara_phy_scaled)

The distance matrix dimensions are 139 rows and 139 columns. Indeed, we calculated pairwise distance, so the matrix contains samples in rows and in columns.

dim(tara_nog_dist)
## [1] 139 139

The diagonal is composed of zero values (or values very closed). Indeed, there is no distance between the same sample.

tara_nog_dist[c(1:5), c(1:5)]
##              TARA_109_SRF TARA_149_MES TARA_110_MES TARA_102_MES TARA_142_SRF
## TARA_109_SRF 1.705303e-13     536.8095     329.8441 5.991295e+02     646.4804
## TARA_149_MES 5.368095e+02       0.0000     261.2395 3.325547e+02     738.4168
## TARA_110_MES 3.298441e+02     261.2395       0.0000 2.078134e+02     638.5818
## TARA_102_MES 5.991295e+02     332.5547     207.8134 2.273737e-13     877.6855
## TARA_142_SRF 6.464804e+02     738.4168     638.5818 8.776855e+02       0.0000

High distance values mean that samples are not similar. And small distance values mean that samples are similar.

4.2.2 Similarity calculation

The distance matrix is then transformed into similarity matrix for each data type. We set two parameters:

  • K = 20: number of nearest neighbors
  • sigma = 0.5: hyperparameter
K <- 20
signma <- 0.5

The affinityMatrix() function transforms the distance into similarity according the distance with the nearest neighbors.

tara_nog_W <- affinityMatrix(tara_nog_dist, K, signma)
tara_phy_W <- affinityMatrix(tara_phy_dist, K, signma)

The following figures are the heatmap of the similarity matrix (W) of each data type. The left heatmap are the nog data and the right heatmap are the phy data. Samples are clustered using hierarchical clustering. For a better visualization, we log-transform similarities.

pheatmap(log(tara_nog_W, 10), show_rownames = FALSE, show_colnames = FALSE, annotation = tara_metadata, main = "Orthologous genes - log10-transformed similarity values")
pheatmap(log(tara_phy_W, 10), show_rownames = FALSE, show_colnames = FALSE, annotation = tara_metadata, main = "Phylogenetic profil - log10-transformed similarity values")

Red color means a high similarity value between two samples whereas blue color means a small similarity value between two samples.

The two heatmaps are different. Data seem to probably carry different type of information about the samples:

  • for the nog data (left), we can see two main groups of samples. Each group seems to be composed of the same ocean depth. It’s not clear for the ocean regions.
  • for the phy data (right), we do not identify clear groups. But, samples seem to be group by the depth too.

4.3 Fusion

We created a similarity matrix for each data type. We saw that each network carries common information and its own information. Now, we will integrate all this information into only one fused similarity matrix.

4.3.1 Create the fused similarity matrix

We create the fused similarity matrix using these three parameters:

  • list that contains the nog and phy similarity matrices
  • K = 20: number of the nearest neighbors
  • T = 10: number of iterations
K <- 20
T <- 10
tara_W <- SNF(list(tara_nog_W, tara_phy_W), K, T)
tara_W[c(1:5), c(1:5)]
##              TARA_109_SRF TARA_149_MES TARA_110_MES TARA_102_MES TARA_142_SRF
## TARA_109_SRF 5.000000e-01 9.778842e-05 0.0002535809 0.0004410605 0.0019200171
## TARA_149_MES 9.778842e-05 5.000000e-01 0.0129882372 0.0126485110 0.0004977051
## TARA_110_MES 2.535809e-04 1.298824e-02 0.5000000000 0.0266185584 0.0029147232
## TARA_102_MES 4.410605e-04 1.264851e-02 0.0266185584 0.5000000000 0.0009179439
## TARA_142_SRF 1.920017e-03 4.977051e-04 0.0029147232 0.0009179439 0.5000000000

The dimension of the fused similarity matrix are 139 rows and 139 columns, such as the previous similarity matrices. The fused similarity matrix contains similarities between samples, we can also called them weights.

The fused similarity matrix contains 19321 weights.

length(tara_W)
## [1] 19321

The fused similarity matrix doesn’t contain zero:

length(tara_W[length(tara_W) == 0])
## [1] 0

The following figure is the heatmap of the fused similarity matrix. Samples are automatically clustered with a hierarchical clustering. Weights are log-transformed for a better visualization.

pheatmap(log(tara_W, 10), show_rownames = FALSE, show_colnames = FALSE, annotation = tara_metadata, main = "Tara Ocean data - log10-transformed fused similarity matrix")

Groups are more clearly defined in this fused similarity matrix. One corresponds to the depth MESO and other to a mix of DCM and SRF. The fused similarity matrix seems to be a mix of each similarity data type matrix.

4.3.2 Visualized the fused similarity network

Now, we create a fused similarity network from the fused similarity matrix. Self loops are remove (diag = FALSE) and only the upper values of the matrix are taken (mode = "upper", avoid duplicate information).

tara_W_net <- graph_from_adjacency_matrix(tara_W, weighted = TRUE, mode = "upper", diag = FALSE)

Then, the fused similarity network is saved into a the TaraOcean_W_edgeList.txt file:

write.table(as_data_frame(tara_W_net), "../02_Results/01_TaraOcean/TaraOcean_W_edgeList.txt", quote = FALSE, col.names = TRUE, row.names = FALSE, sep = "\t")
This files is loaded into Cytoscape. The Figure 4.1 shows the fused similarity network of the Tara Ocean dataset.
The fused similarity network of the Tara Ocean dataset.

Figure 4.1: The fused similarity network of the Tara Ocean dataset.

According to Cytoscape, the network contains 139 samples (nodes) and 9591 connections (edges). The number of edges is smaller than in the similarity matrix because in the similarity matrix the weights are duplicates. The similarity matrix contains also the weights for each sample compare to itself (self loops).

For now, the network is fully connected: each sample is connected to every sample. Connections between samples are weights: some connections are strong (samples are similar) some other are weak (samples are not similar).

4.4 Threshold selection

So in this section, we will choose a threshold to keep the strongest connections.

4.4.1 Fused similarity network

4.4.1.1 Arbitrary threshold

We extract the weight to display the corresponding distribution to try to find a threshold. The distribution in the left is created using the raw weights. The distribution in the right is created using the log-transformed weights.

tara_weights <- edge.attributes(tara_W_net)$weight
hist(tara_weights, nclass = 100, main = "Fused similarity network weight distribution", xlab = "weights")
hist(log(tara_weights, 10), nclass = 100, main = "Fused similarity network weight distribution", xlab = "weights")
abline(v = log(0.001, 10), col = "cyan", lwd = 3)

The log-transformed weight distribution shows a binomial distribution. We would probably like to cut between the two peaks and choose the corresponding weight: 0.001 With this threshold, we select 6440 connections.

4.4.1.2 Mean and third quantile

We calculate the median of the weights.

tara_W_median <- median(x = tara_weights)
tara_W_median
## [1] 0.002143026

With the mean (0.002143) as threshold, we select 4796 connections.

length(tara_weights[tara_weights>=tara_W_median])
## [1] 4796

Calculate the third quantile of the weights:

tara_W_q75 <- quantile(x = tara_weights, 0.75)
tara_W_q75
##         75% 
## 0.004119391

With the third quantile (0.0041194) as threshold, we select 2398 connections.

length(tara_weights[tara_weights>=tara_W_q75])
## [1] 2398

In the following figure, we display the log-transformed weight distribution with the two previous calculated metrics as markers.

hist(log(tara_weights, 10), nclass = 100, main = "Fused similarity network weight distribution", xlab = "log10(weights)")
abline(v = log(tara_W_median, 10), col = "blue", lwd = 3)
text(log(tara_W_median, 10), 400, pos = 2, "Median", col = "blue", cex = 1)
abline(v = log(tara_W_q75, 10), col = "purple", lwd = 3)
text(log(tara_W_q75, 10), 400, pos = 4, "quantile 75%", col = "purple", cex = 1)

4.4.1.3 Topology network

To determine a range of thresholds to try, we check the weights.

summary(c(tara_W))
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000570 0.0004456 0.0021613 0.0071942 0.0042193 0.5000000

We set the range between 0 and 0.0005.

thresholds <- seq(0, 0.1, 0.0005)
length(thresholds)
## [1] 201

Then, we calculate the Average Clustering Coefficient for each threshold.

tara_W_ACC <- do.call(rbind, lapply(thresholds, function(t, net){
  net_sub <- subgraph.edges(net, E(net)[weight >= t])
  df <- data.frame("ACC" = ACCCalculation(net_sub), "thresholds" = t, "EN" = length(E(net_sub)))
  return(df)
}, tara_W_net))

Calculated values are displayed in the following figures.

plot(x = tara_W_ACC$thresholds, y = tara_W_ACC$ACC, xlab = "thresholds", ylab = "ACC", main = "ACC calculation of the Fused network W", type = "o")
points(x = tara_W_ACC$thresholds[1], y = tara_W_ACC$ACC[1], col = "red", pch = 16, cex = 1.2)
points(x = tara_W_ACC$thresholds[21], y = tara_W_ACC$ACC[21], col = "pink", pch = 16, cex = 1.2)
points(x = tara_W_ACC$thresholds[29], y = tara_W_ACC$ACC[29], col = "purple", pch = 16, cex = 1.2)
abline(v = tara_W_ACC$thresholds[29], col = "purple")
text(tara_W_ACC$thresholds[29], 0.7, pos = 4, paste0("Threshold = ",  tara_W_ACC$thresholds[29]), col = "purple")
text(tara_W_ACC$thresholds[29], 0.6, pos = 4, paste0("ACCmax = ",  round(tara_W_ACC$ACC[29], 2)), col = "purple")
plot(x = tara_W_ACC$thresholds, y = tara_W_ACC$EN, xlab = "thresholds", ylab = "number of edges", main = "EN of the Fused network W", type = "o")
abline(v = tara_W_ACC$thresholds[29], col = "purple")
text(tara_W_ACC$thresholds[29], 2800, pos = 4, paste0("Threshold = ",  tara_W_ACC$thresholds[29]), col = "purple")
text(tara_W_ACC$thresholds[29], 2000, pos = 4, paste0("ACCmax = ",  round(tara_W_ACC$ACC[29], 2)), col = "purple")
text(tara_W_ACC$thresholds[29], 1200, pos = 4, paste0("EN = ",  tara_W_ACC[29, "EN"]), col = "purple")
**Left**: Average Clustering Coeeficient (ACC) values for each threshold - **Right**: Number of edges on network for each threshold**Left**: Average Clustering Coeeficient (ACC) values for each threshold - **Right**: Number of edges on network for each threshold

Figure 4.2: Left: Average Clustering Coeeficient (ACC) values for each threshold - Right: Number of edges on network for each threshold

  • The red dot value corresponds to the fully connected network.
  • The pink dot value is the smallest value before the local maxima.
  • The purple dot value is the local maxima. The one we are interested in.

If we selected the purple local maxima, we will have 404 edges. It could be not enough edges. Let’s see during the visualization.

tara_W_ACC$thresholds[29]
## [1] 0.014

4.4.1.4 Visualization using Cytoscape

The network visualization on the left was created with the third quantile (0.0041194). The network visualization on the right was created with the ACC method (0.014).