Skip to content

Commit 4f2e275

Browse files
giorgiagandolfiSPPearceGiorgia Gandolfi
authored
add viber module (#8797)
* add viber module * mc parameters viber * new snapshot * add conda env specification * add viber module * mc parameters viber * new snapshot * add conda env specification * changes after reviews viber * Update modules/nf-core/viber/main.nf Indentation issue Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * Update modules/nf-core/viber/main.nf Indentation issue Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * Update modules/nf-core/viber/main.nf Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * Update modules/nf-core/viber/main.nf Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * Update modules/nf-core/viber/templates/viber_main_script.R Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * Update modules/nf-core/viber/templates/viber_main_script.R Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * Update modules/nf-core/viber/templates/viber_main_script.R Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * Update main.nf Change name of the report output and adjust the `emit`s in the output channel * Update meta.yml * Update main.nf * Update meta.yml * Update viber_main_script.R * mc main.nf * mc modules/nf-core/viber/main.nf Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * mc modules/nf-core/viber/templates/viber_main_script.R Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> * change filtering strategy in viber script * add snapshot * final changes viber module * add updated snapshot * revert to r-cnaqc 1.1.2 version --------- Co-authored-by: Simon Pearce <24893913+SPPearce@users.noreply.github.com> Co-authored-by: Giorgia Gandolfi <giorgiagandolfi@215.52.dmg.units.it>
1 parent 3fa617b commit 4f2e275

File tree

7 files changed

+825
-0
lines changed

7 files changed

+825
-0
lines changed

‎modules/nf-core/viber/environment.yml‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
channels:
2+
- conda-forge
3+
- bioconda
4+
5+
dependencies:
6+
- bioconda::r-cnaqc=1.1.2
7+
- bioconda::r-viber=1.0.0
8+
- conda-forge::r-cli=3.6.5
9+
- conda-forge::r-dplyr=1.1.4
10+
- conda-forge::r-ggplot2=3.5.2
11+
- conda-forge::r-ggpubr=0.6.1
12+
- conda-forge::r-tidyr=1.3.1

‎modules/nf-core/viber/main.nf‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
process VIBER {
2+
tag "$meta.id"
3+
label "process_single"
4+
5+
conda "${moduleDir}/environment.yml"
6+
container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
7+
'oras://community.wave.seqera.io/library/r-cnaqc_r-viber:014077a3164189d5':
8+
'community.wave.seqera.io/library/r-cnaqc_r-viber:2314592f7d2f9abe'}"
9+
10+
input:
11+
tuple val(meta), path(rds_join), val(tumour_samples)
12+
13+
output:
14+
tuple val(meta), path("*_viber_best_fit.rds") , emit: viber_rds
15+
tuple val(meta), path("*_viber_best_heuristic_fit.rds") , emit: viber_heuristic_rds
16+
tuple val(meta), path("*_viber_best_fit_plots.rds") , emit: viber_plots_rds
17+
tuple val(meta), path("*_viber_best_heuristic_fit_plots.rds"), emit: viber_heuristic_plots_rds
18+
tuple val(meta), path("*_viber_report.rds") , emit: viber_report_rds
19+
tuple val(meta), path("*_viber_report.pdf") , emit: viber_report_pdf
20+
tuple val(meta), path("*_viber_report.png") , emit: viber_report_png
21+
path "versions.yml" , emit: versions
22+
23+
when:
24+
task.ext.when == null || task.ext.when
25+
26+
script:
27+
template "viber_main_script.R"
28+
29+
stub:
30+
def prefix = task.ext.prefix ?: "${meta.id}"
31+
32+
"""
33+
touch ${prefix}_viber_best_fit.rds
34+
touch ${prefix}_viber_best_heuristic_fit.rds
35+
touch ${prefix}_viber_best_fit_plots.rds
36+
touch ${prefix}_viber_best_heuristic_fit_plots.rds
37+
touch ${prefix}_viber_report.rds
38+
touch ${prefix}_viber_report.pdf
39+
touch ${prefix}_viber_report.png
40+
41+
cat <<-END_VERSIONS > versions.yml
42+
"${task.process}":
43+
VIBER: \$(Rscript -e 'library(VIBER); sessionInfo()\$otherPkgs\$VIBER\$Version')
44+
cli: \$(Rscript -e 'library(cli); sessionInfo()\$otherPkgs\$cli\$Version')
45+
dplyr: \$(Rscript -e 'library(dplyr); sessionInfo()\$otherPkgs\$dplyr\$Version')
46+
ggplot2: \$(Rscript -e 'library(ggplot2); sessionInfo()\$otherPkgs\$ggplot2\$Version')
47+
END_VERSIONS
48+
"""
49+
}

‎modules/nf-core/viber/meta.yml‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
name: viber
2+
description: Multisample subclonal deconvolution of cancer genome sequencing data.
3+
keywords:
4+
- subclonal deconvolution
5+
- genomics
6+
- cancer evolution
7+
tools:
8+
- "viber":
9+
description: |
10+
VIBER is a package that implements a variational Bayesian model to fit multi-variate Binomial mixtures.
11+
The statistical model is semi-parametric and fit by a variational mean-field approximation to the model posterior.
12+
The components are Binomial distributions which can model count data;
13+
these can be used to model sequencing counts in the context of cancer, for instance.
14+
The package implements methods to fit and visualize clustering results.
15+
homepage: "https://caravagnalab.github.io/VIBER/"
16+
documentation: "https://caravagnalab.github.io/VIBER/"
17+
tool_dev_url: "https://github.com/caravagnalab/VIBER/"
18+
doi: "10.1038/s41588-020-0675-5"
19+
licence: ["GPL v3-or-later"]
20+
identifier: ""
21+
22+
input:
23+
- - meta:
24+
type: map
25+
description: |
26+
Groovy Map containing sample information
27+
e.g. `[ id:'sample1']`
28+
- rds_join:
29+
type: file
30+
description: Either a .rds object of class mCNAqc or a .csv mutations table
31+
pattern: "*.{rds,csv}"
32+
ontologies:
33+
- edam: "http://edamontology.org/format_3752" # csv
34+
- tumour_samples:
35+
type: list
36+
description: List of tumour sample identifiers for a specific patient
37+
output:
38+
viber_rds:
39+
- - meta:
40+
type: map
41+
description: |
42+
Groovy Map containing sample information
43+
e.g. `[ id:'sample1' ]`
44+
- "*_viber_best_fit.rds":
45+
type: file
46+
description: Best viber fit as an .rds object
47+
pattern: "*_viber_best_fit.rds"
48+
ontologies: []
49+
viber_heuristic_rds:
50+
- - meta:
51+
type: map
52+
description: |
53+
Groovy Map containing sample information
54+
e.g. `[ id:'sample1' ]`
55+
- "*_viber_best_heuristic_fit.rds":
56+
type: file
57+
description: Cluster-filtering heuristics viber fit as an .rds object
58+
pattern: "*_viber_best_heuristic_fit.rds"
59+
ontologies: []
60+
viber_plots_rds:
61+
- - meta:
62+
type: map
63+
description: |
64+
Groovy Map containing sample information
65+
e.g. `[ id:'sample1' ]`
66+
- "*_viber_best_fit_plots.rds":
67+
type: file
68+
description: best fit plot
69+
pattern: "_viber_best_fit_plots.rds"
70+
ontologies: []
71+
viber_heuristic_plots_rds:
72+
- - meta:
73+
type: map
74+
description: |
75+
Groovy Map containing sample information
76+
e.g. `[ id:'sample1' ]`
77+
- "*_viber_best_heuristic_fit_plots.rds":
78+
type: file
79+
description: heuristics fit plot
80+
pattern: "*_viber_best_heuristic_fit_plots.rds"
81+
ontologies: []
82+
viber_report_rds:
83+
- - meta:
84+
type: map
85+
description: |
86+
Groovy Map containing sample information
87+
e.g. `[ id:'sample1' ]`
88+
- "*_viber_report.rds":
89+
type: file
90+
description: final report plots as a .pdf file
91+
pattern: "*_viber_report.rds"
92+
ontologies: []
93+
viber_report_pdf:
94+
- - meta:
95+
type: map
96+
description: |
97+
Groovy Map containing sample information
98+
e.g. `[ id:'sample1' ]`
99+
- "*_viber_report.pdf":
100+
type: file
101+
description: final report plots as a .pdf file
102+
pattern: "*_viber_report.pdf"
103+
ontologies: []
104+
viber_report_png:
105+
- - meta:
106+
type: map
107+
description: |
108+
Groovy Map containing sample information
109+
e.g. `[ id:'sample1' ]`
110+
- "*_viber_report.png":
111+
type: file
112+
description: final report plots as a .png file
113+
pattern: "*_viber_report.png"
114+
ontologies: []
115+
versions:
116+
- versions.yml:
117+
type: file
118+
description: File containing software versions
119+
pattern: "versions.yml"
120+
121+
ontologies:
122+
- edam: http://edamontology.org/format_3750 # YAML
123+
authors:
124+
- "@giorgiagandolfi"
125+
maintainers:
126+
- "@giorgiagandolfi"

‎modules/nf-core/viber/templates/viber_main_script.R‎

Lines changed: 167 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,167 @@
1+
#!/usr/bin/env Rscript
2+
3+
parse_args = function(x) {
4+
x = gsub("\\\\[","",x)
5+
x = gsub("\\\\]","",x)
6+
# giving errors when we have lists like c(xxx, xxx) since it will separate it
7+
# args_list = unlist(strsplit(x, ', ')[[1]])
8+
args_list = unlist(strsplit(x, ", (?=[^)]*(?:\\\\(|\$))", perl=TRUE))
9+
# args_vals = lapply(args_list, function(x) strsplit(x, split=":")[[1]])
10+
args_vals = lapply(args_list, function(x) {
11+
x_splt = strsplit(x, split=":")[[1]]
12+
c(x_splt[1], paste(x_splt[2:length(x_splt)], collapse=":"))
13+
})
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+
opt = list(
23+
prefix = ifelse('$task.ext.prefix' == 'null', '$meta.id', '$task.ext.prefix')
24+
)
25+
args_opt = parse_args('$task.ext.args')
26+
for ( ao in names(args_opt)) opt[[ao]] = args_opt[[ao]]
27+
print(opt)
28+
29+
# Script #####
30+
31+
library(VIBER)
32+
library(dplyr)
33+
library(tidyr)
34+
library(ggplot2)
35+
library(CNAqc)
36+
library(ggpubr)
37+
38+
samples = substr("$tumour_samples", 2, nchar("$tumour_samples")-1)
39+
samples = strsplit(samples, ", ")[[1]]
40+
print("$meta.patient")
41+
print("$tumour_samples")
42+
print("$rds_join")
43+
print(samples)
44+
45+
if ( grepl(".rds\$", tolower("$rds_join")) ) {
46+
input_obj = readRDS("$rds_join")
47+
if (class(input_obj) == "m_cnaqc") {
48+
shared = input_obj %>% get_sample(sample=samples, which_obj="shared")
49+
joint_table = lapply(names(shared),
50+
function(sample_name)
51+
shared[[sample_name]] %>%
52+
CNAqc::subset_by_segment_karyotype("1:1") %>%
53+
CNAqc::Mutations() %>%
54+
dplyr::mutate(sample_id=sample_name)
55+
) %>% dplyr::bind_rows()
56+
} else {
57+
cli::cli_alert_warning("Object of class {class(input_obj)} not supported.")
58+
return()
59+
}
60+
} else {
61+
joint_table = read.csv("$rds_join")
62+
}
63+
64+
print("Subset joint done")
65+
66+
## Read input joint table
67+
input_tab = joint_table %>%
68+
dplyr::mutate(VAF=replace(VAF, VAF==0, 1e-7))
69+
70+
## Convert the input table into longer format
71+
reads_data = input_tab %>%
72+
dplyr::select(chr, from, ref, alt, NV, DP, VAF, sample_id,driver_label,is_driver) %>%
73+
dplyr::rename(gene=driver_label) %>%
74+
dplyr::rename(driver=is_driver) %>%
75+
tidyr::pivot_wider(names_from="sample_id",
76+
values_from=c("NV","DP","VAF"), names_sep=".",values_fill=0)
77+
78+
## Extract DP (depth)
79+
dp = reads_data %>%
80+
dplyr::select(dplyr::starts_with("DP")) %>%
81+
dplyr::rename_all(function(x) stringr::str_remove_all(x,"DP."))
82+
83+
## Extract NV (alt_counts)
84+
nv = reads_data %>%
85+
dplyr::select(dplyr::starts_with("NV")) %>%
86+
dplyr::rename_all(function(x) stringr::str_remove_all(x,"NV."))
87+
88+
# Standard fit
89+
viber_K = as.integer(opt[["K"]])
90+
91+
message("Starting standard fit")
92+
st_fit = VIBER::variational_fit(nv, dp,
93+
K=viber_K,
94+
data=reads_data)
95+
st_fit[["description"]]="$meta.patient"
96+
message("End standard fit")
97+
best_fit = best_fit_heuristic = st_fit
98+
99+
# If all clusters are removed -> keep the origianl best fit
100+
tryCatch(expr = {
101+
# Apply the heuristic
102+
best_fit_heuristic = VIBER::choose_clusters(st_fit,
103+
binomial_cutoff=as.numeric(opt[["binomial_cutoff"]]),
104+
dimensions_cutoff=as.integer(opt[["dimensions_cutoff"]]),
105+
pi_cutoff=as.numeric(opt[["pi_cutoff"]]),
106+
re_assign=as.logical(opt[["re_assign"]])
107+
)
108+
}, error = function(e) {
109+
print(e)
110+
best_fit_heuristic <<- st_fit
111+
} )
112+
113+
# Save fits
114+
saveRDS(best_fit, file=paste0(opt[["prefix"]], "_viber_best_fit.rds"))
115+
saveRDS(best_fit_heuristic, file = paste0(opt[["prefix"]], "_viber_best_heuristic_fit.rds"))
116+
117+
# Save plots
118+
n_samples = ncol(best_fit[["x"]]) - 1
119+
if (n_samples >1) {
120+
print("multisample mode on")
121+
plot_fit = plot(best_fit)
122+
plot_fit_heuristic = plot(best_fit_heuristic)
123+
} else {
124+
plot_fit = plot_mixing_proportions(best_fit)
125+
plot_fit_heuristic = plot_mixing_proportions(best_fit_heuristic)
126+
}
127+
128+
saveRDS(plot_fit, file=paste0(opt[["prefix"]], "_viber_best_fit_plots.rds"))
129+
saveRDS(plot_fit_heuristic, file=paste0(opt[["prefix"]], "_viber_best_heuristic_fit_plots.rds"))
130+
131+
# Save report plot
132+
n_samples = ncol(best_fit[["x"]]) - 1
133+
marginals = multivariate = ggplot()
134+
135+
try(expr = {marginals <<- VIBER::plot_1D(best_fit)} )
136+
try(expr = {multivariate = plot(best_fit)})
137+
try(expr = {multivariate = ggpubr::ggarrange(plotlist = multivariate)})
138+
top_p = ggpubr::ggarrange(plotlist = list(marginals, multivariate), widths=ifelse(n_samples>2, c(1,2), c(2,1)))
139+
140+
mix_p = VIBER::plot_mixing_proportions(best_fit)
141+
binom = VIBER::plot_peaks(best_fit)
142+
elbo = VIBER::plot_ELBO(best_fit)
143+
bottom_p = ggpubr::ggarrange(plotlist = list(mix_p, binom, elbo), widths = c(1,2,1), nrow = 1)
144+
145+
report_fig = ggpubr::ggarrange(top_p, bottom_p, nrow = 2, heights = ifelse(n_samples>2, c(3,1), c(2,1)))
146+
saveRDS(report_fig, file=paste0(opt[["prefix"]], "_viber_report.rds"))
147+
ggplot2::ggsave(plot=report_fig, filename=paste0(opt[["prefix"]], "_viber_report.pdf"), height=210, width=210, units="mm", dpi = 200)
148+
ggplot2::ggsave(plot=report_fig, filename=paste0(opt[["prefix"]], "_viber_report.png"), height=210, width=210, units="mm", dpi = 200)
149+
150+
# version export
151+
f = file("versions.yml","w")
152+
cnaqc_version = sessionInfo()\$otherPkgs\$CNAqc\$Version
153+
viber_version = sessionInfo()\$otherPkgs\$VIBER\$Version
154+
cli_version = sessionInfo()\$otherPkgs\$cli\$Version
155+
dplyr_version = sessionInfo()\$otherPkgs\$dplyr\$Version
156+
tidyr_version = sessionInfo()\$otherPkgs\$tidyr\$Version
157+
ggplot2_version = sessionInfo()\$otherPkgs\$ggplot2\$Version
158+
ggpubr_version = sessionInfo()\$otherPkgs\$ggpubr\$Version
159+
writeLines(paste0('"', "$task.process", '"', ":"), f)
160+
writeLines(paste(" CNAqc:", cnaqc_version), f)
161+
writeLines(paste(" VIBER:", viber_version), f)
162+
writeLines(paste(" cli:", cli_version), f)
163+
writeLines(paste(" dplyr:", dplyr_version), f)
164+
writeLines(paste(" tidyr:", tidyr_version), f)
165+
writeLines(paste(" ggplot2:", ggplot2_version), f)
166+
writeLines(paste(" ggpubr:", ggpubr_version), f)
167+
close(f)

0 commit comments

Comments
 (0)