new module: wgcna/automaticsoftthresholding

[thomgiles]

Is there an existing module for this?

• [X] I have searched for the existing module

Is there an open PR for this?

• [X] I have searched for existing PRs

Is there an open issue for this?

• [X] I have searched for existing issues

Are you going to work on this?

• [X] If I'm planning to work on this module, I added myself to the Assignees to facilitate tracking who is working on the module

[EuancRNA]

Have written an R script to get the soft-thresholding value, but unsure if it makes sense to return that and use as input in the next module (probably "network construction" or the subsequent "blockwise" WGCNA step) as the adjacency matrices could get very large, esp if methylation-level.

[thomgiles]

So I did this as a Rscript as a single value in testing as essentially only a single value is outputted. could also be done with a nested template using stdout to better fit the normal practice?

script:

    """
    softPower=\$(Rscript -e \\
    " \\
    zz <- file(\\"messages.Rout\\", open = \\"wt\\"); \\
    sink(zz, type = \\"message\\"); \\
    sink(zz); \\
    library(WGCNA); \\
    library(doParallel); \\
    registerDoParallel(makeCluster(4)); \\
    sft_threshold <- as.numeric('${params.sft_threshold}'); \\
    matrix_path <- '${matrix}'; \\
    tab_seperator <- \\\"\\\\\\\\t\\\"; \\
    matrix <- read.delim(matrix_path, sep = tab_seperator, header = TRUE, row.names = 1); \\
    matrix_transformed = t(matrix); \\
    powers <- c(c(1:10), seq(from = 12, to=30, by=1)); \\
    sft <- pickSoftThreshold(matrix_transformed,dataIsExpr = TRUE,corFnc = cor,networkType = \\"signed\\",powerVector = powers,verbose = 4); \\
    sft_df <- data.frame(sft[[\\"fitIndices\\"]]); \\
    sft_df[[\\"model_fit\\"]] <- -sign(sft_df[[\\"slope\\"]]) * sft_df[[\\"SFT.R.sq\\"]]; \\
    r2val <- sft_df[[\\"SFT.R.sq\\"]]; \\
    slope <- sft_df[[\\"slope\\"]]; \\
    softPower <- sft_df[r2val >= sft_threshold,]; \\
    softPower <- softPower[[\\"Power\\"]][1]; \\
    if(is.na(softPower)){softPower = 20}; \\
    png(filename = \\"scale_independence.png\\", width = 800, height = 600); \\
    print(plot(sft[[\\"fitIndices\\"]][,1], -sign(sft[[\\"fitIndices\\"]][,3])*sft[[\\"fitIndices\\"]][,2], \\
           xlab=\\"Soft Threshold (power)\\", ylab=\\"Scale Free Topology Model Fit,signed R^2\\", type=\\"n\\", \\
           main = \\"Scale independence\\") + \\
        text(sft[[\\"fitIndices\\"]][,1], -sign(sft[[\\"fitIndices\\"]][,3])*sft[[\\"fitIndices\\"]][,2], \\
             labels=powers, cex=sft_threshold, col=\\"red\\") + \\
        abline(h=sft_threshold, col=\\"red\\")); \\
    dev.off(); \\
    writeLines(capture.output(sessionInfo()), 'R_sessionInfo.log'); \\
    writeLines(c('R-wgcna:1.71--r43h21a89ab_4'), 'versions.yml');\\
    sink(); \\
    cat (softPower);")""" 
}