|
| 1 | +#!/usr/bin/env Rscript |
| 2 | + |
| 3 | +# Ported to nf-core/modules with template by Jonathan Manning |
| 4 | + |
| 5 | +#' Parse out options from a string without recourse to optparse |
| 6 | +#' |
| 7 | +#' @param x Long-form argument list like --opt1 val1 --opt2 val2 |
| 8 | +#' |
| 9 | +#' @return named list of options and values similar to optparse |
| 10 | + |
| 11 | +parse_args <- function(x){ |
| 12 | + args_list <- unlist(strsplit(x, ' ?--')[[1]])[-1] |
| 13 | + args_vals <- lapply(args_list, function(x) scan(text=x, what='character', quiet = TRUE)) |
| 14 | + |
| 15 | + # Ensure the option vectors are length 2 (key/ value) to catch empty ones |
| 16 | + args_vals <- lapply(args_vals, function(z){ length(z) <- 2; z}) |
| 17 | + |
| 18 | + parsed_args <- structure(lapply(args_vals, function(x) x[2]), names = lapply(args_vals, function(x) x[1])) |
| 19 | + parsed_args[! is.na(parsed_args)] |
| 20 | +} |
| 21 | + |
| 22 | +# Load |
| 23 | +library("RTN") |
| 24 | +library("snow") |
| 25 | + |
| 26 | +################################################ |
| 27 | +################################################ |
| 28 | +## Pull in module inputs ## |
| 29 | +################################################ |
| 30 | +################################################ |
| 31 | + |
| 32 | +input_expr_matrix <- '$expression_matrix' |
| 33 | +output_prefix = ifelse('$task.ext.prefix' == 'null', '$meta.id', '$task.ext.prefix') |
| 34 | +threads <- $task.cpus |
| 35 | +args_opt <- parse_args('$task.ext.args') |
| 36 | + |
| 37 | +n_perm <- ifelse('n_permutations' %in% names(args_opt), strtoi(args_opt[['n_permutations']]), 10) |
| 38 | + |
| 39 | +# Debug messages (stdout) |
| 40 | +sink(stdout(), type = "message") # sink messages to stdout |
| 41 | +message("Expression matrix file : ", input_expr_matrix) |
| 42 | +message("Nb permutations : ", n_perm) |
| 43 | +message("Nb threads : ", threads) |
| 44 | +message("Output basename : ", output_prefix) |
| 45 | +if ('tfs' %in% names(args_opt)) { |
| 46 | + message("TFs : ", args_opt[['tfs']]) |
| 47 | + tfs <- strsplit(args_opt[['tfs']], ',') |
| 48 | +} else { |
| 49 | + # Load data |
| 50 | + data(tfsData) |
| 51 | + tfs <- tfsData\$Lambert2018\$SYMBOL |
| 52 | +} |
| 53 | +sink(NULL, type="message") # close the sink |
| 54 | + |
| 55 | +# Preprocess |
| 56 | +# Input 1: 'expData', a named gene expression matrix (genes on rows, samples on cols); |
| 57 | +# Input 2: 'regulatoryElements', a vector listing genes regarded as TFs |
| 58 | +# Input 3: 'rowAnnotation', an optional data frame with gene annotation |
| 59 | +# Input 4: 'colAnnotation', an optional data frame with sample annotation |
| 60 | + |
| 61 | +exp_data <- read.csv(input_expr_matrix, sep='\t') |
| 62 | +rownames(exp_data) <- exp_data[,1] |
| 63 | +rowAnnotation <- exp_data[,1:2] |
| 64 | +colnames(rowAnnotation) <- c('PROBEID', 'SYMBOL') |
| 65 | +rowAnnotation\$SYMBOL <- toupper(rowAnnotation\$SYMBOL) |
| 66 | +exp_data[,1:2] <- NULL |
| 67 | + |
| 68 | +# Regulatory Transcriptional Network Inference |
| 69 | +tfs <- c('ENSG00000125798', 'ENSG00000125816') |
| 70 | +rtni <- tni.constructor(expData = as.matrix(exp_data), |
| 71 | + regulatoryElements = tfs, |
| 72 | + rowAnnotation = rowAnnotation) |
| 73 | + |
| 74 | +options(cluster=snow::makeCluster(spec=threads, "SOCK")) |
| 75 | + |
| 76 | +# Please set nPermutations >= 1000 |
| 77 | +rtni_permutation <- tni.permutation(rtni, nPermutations = n_perm, pValueCutoff = 1e-4) |
| 78 | + |
| 79 | +# Unstable interactions are subsequently removed by bootstrap analysis using the |
| 80 | +# tni.bootstrap() function, which creates a consensus bootstrap network, referred |
| 81 | +# here as refnet (reference network). |
| 82 | +rtni_bootstrapped <- tni.bootstrap(rtni_permutation) |
| 83 | + |
| 84 | +stopCluster(getOption("cluster")) |
| 85 | + |
| 86 | +# remove the weakest interaction in any triplet formed by two TFs and a common |
| 87 | +# target gene, preserving the dominant TF-target pair (ARACNe) |
| 88 | +rtni_filtered <- tni.dpi.filter(rtni_bootstrapped) |
| 89 | + |
| 90 | +saveRDS(rtni, file = "tni.rds") |
| 91 | +saveRDS(rtni_permutation, file = "tni_permutated.rds") |
| 92 | +saveRDS(rtni_bootstrapped, file = "tni_bootstrapped.rds") |
| 93 | +saveRDS(rtni_filtered, file = "tni_filtered.rds") |
| 94 | + |
| 95 | +# Plot |
| 96 | +#pdf(paste0(output_prefix, "_RTN.pdf")) |
| 97 | +#tni.graph(rtni_filtered, regulatoryElements = c("FOXM1", "E2F2")) |
| 98 | +#title("Regulatory Transcriptional Network") |
| 99 | +#mtext(output_prefix, side=3) |
| 100 | +#dev.off() |
| 101 | +#cat( |
| 102 | +# paste("- Threads::", threads), |
| 103 | +# fill=TRUE, labels=output_prefix, |
| 104 | +# file=paste0(output_prefix, "_intercept_slope.txt"), append=FALSE |
| 105 | +#) |
| 106 | + |
| 107 | +################################################ |
| 108 | +################################################ |
| 109 | +## R SESSION INFO ## |
| 110 | +################################################ |
| 111 | +################################################ |
| 112 | + |
| 113 | +sink(paste(output_prefix, "R_sessionInfo.log", sep = '.')) |
| 114 | +print(sessionInfo()) |
| 115 | +sink() |
| 116 | + |
| 117 | +################################################ |
| 118 | +################################################ |
| 119 | +## VERSIONS FILE ## |
| 120 | +################################################ |
| 121 | +################################################ |
| 122 | + |
| 123 | +r.version <- strsplit(version[['version.string']], ' ')[[1]][3] |
| 124 | +rtn.version <- as.character(packageVersion('RTN')) |
| 125 | + |
| 126 | +writeLines( |
| 127 | + c( |
| 128 | + '"${task.process}":', |
| 129 | + paste(' bioconductor-rtn:', rtn.version) |
| 130 | + ), |
| 131 | +'versions.yml') |
| 132 | + |
| 133 | +################################################ |
| 134 | +################################################ |
| 135 | +################################################ |
| 136 | +################################################ |
0 commit comments