Input data must be batch-effect correct ? Raw ?

[Simon-chevolleau]

I use NetRep with a batch-effect corrected datasets (using fastMNN), those data are corrected, centered and scaled during fastMNN process. For the adjacency, I simply use the function from WGCNA and for the network:

# SAVE NETWORK AS EDGE/VERTICES FILE
edge_list <- as.data.frame(TOM)
edge_list$gene1 <- rownames(edge_list)

# RESHAPE THE DATA FROM WIDE TO LONG FORMAT
edge_list <- reshape(edge_list,
                     varying = names(edge_list)[-ncol(edge_list)],
                     v.names = "correlation",
                     timevar = "gene2",
                     times = names(edge_list)[-ncol(edge_list)],
                     idvar = "gene1",
                     direction = "long")

edge_list <- edge_list[edge_list$gene1 != edge_list$gene2, ]

edge_list$module1 <- gene_modules[edge_list$gene1, "Module"]
edge_list$module2 <- gene_modules[edge_list$gene2, "Module"]

# only keep correlations > 0.1
selection <- edge_list$correlation > 0.1
table(selection)
edge_list <- edge_list[selection, ]

edge_list <- unique(edge_list)

Results from those data are quite weird, see my plots (calculated on reference dataset for the first one and test dataset for the second).
Do you know what I did wrong there ? should I keep higher correlation for the network ? How can I determine the best threshold for correlation? Should I use raw data with batch inside it ?

Thank for the package ! It looks very promising.

[sritchie73]

Hi @Simon-chevolleau

I'm not familiar with the fastMNN method, but based on the plots you've shown it looks like something has gone wrong in the step where you construct the adjacency matrices - it looks like nearly all possible edges are being given a value close to 1 (i.e. indicating near perfect correlation).

For context, here's an example of what a more typical network would look like:

As a general rule of thumb the edge weights in the adjacency matrix should be (1) monotonically related to the correlation between the corresponding pair of nodes, and (2) be smaller in magnitude than the correlation between the corresponding pair of nodes.

For example, a typical analysis with WGCNA on its own might look something like this:

corr_mat <- cor(dat) # columns and genes, rows are samples
adj_mat <- corr_mat^6 # shrink all correlations towards 0 by raising the power of 6

[Simon-chevolleau]

Hi @sritchie73, I may have misunterpreted how adjacency function works ... Here is my piece of code, where I use adjacency directly on expression matrix, but I set the softpower threshold instead of 6 (as you did and in the function aswell...), and so it also affect my network construction, because I use the TOM. Thank you a lot, I will try using adjcency with default power value (6). Another question, do you agree on the construction of the network ?

adjacency <- adjacency(expression_matrix, power = softpower)
TOM <- as.matrix(TOMsimilarity(adjacency))
rownames(TOM) <- colnames(expression_matrix)
colnames(TOM) <- colnames(expression_matrix)
dissTOM = 1 - TOM

writeTable(adjacency, adjacency_path)
writeTable(TOM, TOM_path)
writeTable(dissTOM, dissTOM_path)

[sritchie73]

Hi @Simon-chevolleau

By any chance are you inputting the TOM as your network into NetRep? When used with WGCNA the network argument expects the adjacency matrix, not the TOM.

[Simon-chevolleau]

Hi again @sritchie73
Yes, I was inputting the dissTOM as my network...
Just to be sure of my inputs of the modulePreservation function:

modulePreservation(
 network=network_list, data=data_list, correlation=correlation_list,
 moduleAssignments=modules_labels,
 nPerm=num_permutations, nThreads=num_threads
)

• network: list of adjacency matrices, as said in your last reply.
• correlation: list of correlation matrices obtained from cor(expression_matrix).
• data: list of batch-effect corrected data.

Thank you a lot for those feedbacks!

[Simon-chevolleau]

Sorry to repost, but is my last post correct ?

[sritchie73]

Hi Simon,

Yes, your last post looks correct. Apologies for missing the question - the email notification sent through to me by GitHub had completely different contents for some reason and didn't include a question!

[Simon-chevolleau]

Thx @sritchie73 ! I will give it a try and let you know if it turns out well !

[Simon-chevolleau]

It worked, however the bottom plot doesn't look at all like the ones in the tutorial.

[sritchie73]

Hi Simon,

The output you've generated looks to me like you've run things correctly - it just looks to me like the module you're plotting may not be one that replicates/is preserved in your test dataset - although its hard to say by eye - have you run the modulePreservation() function to test this statistically?

Best,

Scott

[Simon-chevolleau]

The first one is my test dataset and the second one is my ref dataset. This module is highly preserved (7 metrics with pval < 0.05) and is one with the best looking i would say. Each summary profile of all my module look empty like this one, with transparent lines sometimes.
I'm unsure the analysis is good to interpret or if it's there is still something wrong.

[sritchie73]

Hi Simon,

A few thoughts:

(1) As a general rule as modules get larger, the less visually compelling their topological replication is going to look. What you've shown above looks about right in the test data for a module of > 1000 genes.
(2) However, I would say I'm surprised that a module of that size is so tightly correlated in your discovery dataset.
(3) RE: preservation P-values, you may wish to set a more conservative P-value threshold, e.g. Bonferroni correcting for the total number of modules. In case of strong evidence of preservation, I'd normally expect to see permutation p-values coming back as the minimum possible value for the given number of permutations being run.
(4) Re: the transparent lines, these are usually technical artefacts of the limitations of the way the plots have been saved; the combination of small heatmap cell size/bar width and plotting resolution means that transparent lines can sometimes appear. You can try increasing the resolution of the plots although this will make the file sizes bigger.
(5) RE: the empty module data, I suspect this is a normalization issue - i.e. there are some outlier samples/genes in your data that are determining the range of the color gradient, meaning most of your samples will have data close to the middle and thus white on the heatmap. If you just want to fix this for plotting purposes, the plotModule function has a 'dataRange' argument that allows you to manually set the range of values you want the color gradient to use. However, its also worth considering the impact outlier samples/genes might be having on your analysis more broadly.

[Simon-chevolleau]

Hi @sritchie73 ,

1. You're right, many of our modules are large, from hundred to thousand genes.
2. We want to observe the module preservation between original cell types and new stem cells model that are supposed to mimic the original cell types.
3. Bonferroni is very conservative when dealing with large number of tests, I would maybe go with another one less conservative because we have approximately 60 modules. Does NetRep is already configured with a pvalue correction ?
4. When plotting complex heatmap, we usually use color scale where value above 95% are the same color and same for 5%, it removes the outlier impact. Also, our data are normalized, batch-effect corrected, logged, centered and scaled. Our data range from ~ -3.5 to +3.5.

Does the module data plot, is simply a heatmap representing the data passed to the data parameter of the modulePreservation function?

[sritchie73]

Hi @Simon-chevolleau

1. The module data plot is simply a heatmap representing the data passed to the the data parameter, filtered to the nodes in the module(s) (and samples) being plotted. The range on the heatmap legend is automatically determined by the range of values in that subset of the data, unless you manually override using the 'dataRange' argument.

2. There's no in-built P-value correction so you can flexibly choose the relevant test statistics of interest and apply your own correction methods. Bonferroni correction for 60 modules (P < 0.0008) would be in line with the approach we've outlined in the past: https://www.cell.com/fulltext/S2405-4712(16)30217-4 and would only require 1200 permutations (we always recommend running with at least 10k permutations).

3. That being said, given your new stem cells are supposed to mimic the original cell types, should you be expecting much stronger visual evidence of topological replication? As someone not familiar with the technology and success rates for creating new cell lines, naively I would hope that the topology would look almost identical, e.g. similar to if you were to randomly split your original data into 50/50 discovery/testing subsets and running the same network inference + module preservation analysis.