diff --git a/.github/workflows/linting.yml b/.github/workflows/linting.yml index 30e66026..9c2912ae 100644 --- a/.github/workflows/linting.yml +++ b/.github/workflows/linting.yml @@ -32,6 +32,8 @@ jobs: - name: Install Nextflow uses: nf-core/setup-nextflow@v2 + with: + version: "25.10.4" - uses: actions/setup-python@e797f83bcb11b83ae66e0230d6156d7c80228e7c # v6 with: diff --git a/.nf-core.yml b/.nf-core.yml index 7f7b5a8d..72e486c1 100644 --- a/.nf-core.yml +++ b/.nf-core.yml @@ -1,6 +1,8 @@ lint: multiqc_config: - report_comment + files_unchanged: + - .github/workflows/linting.yml nf_core_version: 3.4.1 repository_type: pipeline template: diff --git a/CITATIONS.md b/CITATIONS.md index b539fef1..9b3d84ec 100644 --- a/CITATIONS.md +++ b/CITATIONS.md @@ -64,6 +64,10 @@ > Kharchenko PK, Tolstorukov MY, Park PJ "Design and analysis of ChIP-seq experiments for DNA-binding proteins" Nat. Biotech. doi:10.1038/nbt.1508 +- [CALDER2](https://github.com/CSOgroup/CALDER2) + + > Liu, Y., Nanni, L., Sungalee, S. et al. Systematic inference and comparison of multi-scale chromatin sub-compartments connects spatial organization to cell phenotypes. Nat Commun 12, 2439 (2021). doi:10.1038/s41467-021-22666-3 + ## Software packaging/containerisation tools - [Anaconda](https://anaconda.com) diff --git a/README.md b/README.md index 282189ad..3635fbad 100644 --- a/README.md +++ b/README.md @@ -38,9 +38,10 @@ Here is an outline of the analysis steps: 4. Mark duplicate reads ([`picard Markduplicates`](http://broadinstitute.github.io/picard)) 5. Filter reads and generate alignment statistics ([`samtools`](http://www.htslib.org/)) 6. Create single track profiles in bigwig format ([`deeptools bamCoverage`](https://deeptools.readthedocs.io/en/latest/)) -7. (Optional) Generate pairwise comparison tracks in bigWig format ([`spp`](https://github.com/hms-dbmi/spp)) or ([`deeptools bigwigCompare`](https://deeptools.readthedocs.io/en/develop/content/tools/bigwigCompare.html)). -8. (Optional) Identify differentially enriched solubility regions as described in [Wang et al., 2024 ](https://doi.org/10.1038/s41594-025-01622-5). -9. Generate an analysis report by collecting all generated QC and statistics ([`MultiQC`](http://multiqc.info/)) +7. (Optional) Perform chromatin compartmentalization analysis to identify A/B compartments across chromosomes ([`CALDER2`](https://github.com/CSOgroup/CALDER2)). +8. (Optional) Generate pairwise comparison tracks in bigWig format ([`spp`](https://github.com/hms-dbmi/spp)) or ([`deeptools bigwigCompare`](https://deeptools.readthedocs.io/en/develop/content/tools/bigwigCompare.html)). +9. (Optional) Identify differentially enriched solubility regions as described in [Wang et al., 2025 ](https://doi.org/10.1038/s41594-025-01622-5). +10. Generate an analysis report by collecting all generated QC and statistics ([`MultiQC`](http://multiqc.info/))

diff --git a/bin/differential_solubility_functions.R b/bin/differential_solubility_functions.R deleted file mode 100644 index f31f9ce8..00000000 --- a/bin/differential_solubility_functions.R +++ /dev/null @@ -1,453 +0,0 @@ -#!/usr/bin/env Rscript -# differential_solubility_functions.R - -##################################################################### -## HELPER FUNCTIONS -##################################################################### - -fun1 <- function(lst, n){ - sapply(lst, `[`, n) -} - -##################################################################### -## IMPORT AND REBIN BIGWIG FUNCTION -##################################################################### - -import_and_rebin__bw <- function(files, bin_list, names, genome = NULL) { - bws <- lapply(files, function(file) { - bwR <- rtracklayer::import(file, format = "BigWig", as = "RleList") - - bin_names <- GenomeInfoDb::seqlevels(bin_list) - bw_names <- names(bwR) - chr_order <- sapply(paste0("^", bw_names, "$"), function(chr) grep(chr, bin_names)) - bins_for_bw <- bin_list - GenomeInfoDb::seqlevels(bins_for_bw) <- bin_names[as.vector(unlist(chr_order))] - - bw <- GenomicRanges::binnedAverage( - bins = bins_for_bw, - numvar = bwR[ GenomeInfoDb::seqlevels(bins_for_bw) ], - varname = "score" - ) - - # CHANGED -> set genome tag only if provided (we will provide the genome parameter later...) - if (!is.null(genome)) GenomeInfoDb::genome(bw) <- genome - - bw - }) - names(bws) <- names - bws -} - -##################################################################### -## CHECK SIGN FUNCTION -##################################################################### - -check_sign <- function(x,meann){ - if ( sign(x) == sign(meann) ){ - return("constant_solubility") - } else if (sign(x) > sign(meann)) { - return("shift_increase") - }else if (sign(x) < sign(meann)) { - return("shift_decrease") - } -} - -##################################################################### -## CONFIDENCE INTERVAL FUNCTION -##################################################################### -##standard error - -confidence_interval <- function(vector, nm="prove") { - # Standard deviation of sample - vec_serr <- sd(vector)/sqrt(length(vector)) - vec_serr2x<- vec_serr*2 - # Sample size - n <- length(vector) - # Mean of sample - vec_mean <- mean(vector) - nm_confint_low<-paste0(nm,"_serrx2_lower") - nm_confint_up<-paste0(nm,"_serrx2_upper") - name_mean <- paste0(nm,"_mean") - name_serr <- paste0(nm,"_serrX2") - result <- c(nm_confint_low = vec_mean - vec_serr2x, - nm_confint_up = vec_mean + vec_serr2x, - name_mean = vec_mean, - name_serr = vec_serr2x - ) - names(result) <- c(nm_confint_low,nm_confint_up,name_mean,name_serr) - return(result) -} - -##################################################################### -## RANGE CHECK FUNCTION -##################################################################### - -is_in_serrx2_range_and_shift <- function(vector, - xgroup_name, - y_name - ){ - #meann,lower,upper,mean_tp_name,whereis_tp_name){ - xgroup_lower_bound_confint_name<-paste0(xgroup_name,"_serrx2_lower") - xgroup_upper_bound_confint_name<-paste0(xgroup_name,"_serrx2_upper") - xgroup_mean_name<-paste0(xgroup_name,"_mean") - - xgroup_lower_bound_confint<-vector[xgroup_lower_bound_confint_name][1] - xgroup_upper_bound_confint<-vector[xgroup_upper_bound_confint_name][1] - - yname_val<-vector[y_name][1] - xgroup_mean<-vector[xgroup_mean_name][1] - - confint_check <-'' - mean_startsign<- '' - mean_sign <-'' - shift_solubility<- '' - - col_is_in_confint_name<-paste0(y_name,"_sign") # cambiato da ov_check a _sign - col_whereis_tp_name<- paste0(y_name,"_shift") # cambiato da ov_specs a _shift - shift_solubility_name<- paste0(y_name, "_sol_shift") - mean_startsign_name<- paste0(xgroup_name, "_mean_startsign") - ##########check if both or one mean is in the confint of the other mean - if ( between(yname_val, xgroup_lower_bound_confint, xgroup_upper_bound_confint) - ) { - confint_check <- sign(yname_val) - mean_sign <- "nodiff" - shift_solubility <- "constant_solubility" - mean_startsign<- sign(xgroup_mean) - ##########check if confint are not overlapping, the ygroup timepoint confint is lower than xgroup - }else if (yname_val < xgroup_lower_bound_confint ) { - - confint_check <- sign(yname_val) - - #mean_sign <- paste0(y_name ,"_lower_than_",xgroup_name) - mean_sign <- "lower" - shift_solubility <- check_sign(yname_val,xgroup_mean) - mean_startsign<- sign(xgroup_mean) - - ##########check if confint are not overlapping, the xgroup timepoint confint is lower than ygroup - }else if (xgroup_upper_bound_confint < yname_val) { - - confint_check <-sign(yname_val) - mean_sign <- "higher" - shift_solubility <- check_sign(yname_val,xgroup_mean) - mean_startsign<- sign(xgroup_mean) - } - - else{ - - confint_check <- sign(yname_val) - mean_sign <- "no_idea" - shift_solubility <- check_sign(yname_val,xgroup_mean) - mean_startsign<- sign(xgroup_mean) - } - - result <- c(confint_check,mean_sign,shift_solubility,mean_startsign) - - names(result) <- c(col_is_in_confint_name,col_whereis_tp_name,shift_solubility_name,mean_startsign_name) - return(result) -} - -##################################################################### -## GENE ANNOTATION SETUP FUNCTION -##################################################################### - -setup_gene_annotation <- function(gtf_file) { - # Create TxDb from GTF - txdb <- makeTxDbFromGFF(gtf_file) - genes <- genes(txdb) - - # Function to summarize protein coding genes - summarizeProteinCodingGenes <- function(txdb) { - stopifnot(is(txdb, "TxDb")) - protein_coding_tx <- names(cdsBy(txdb, use.names = TRUE)) - all_tx <- mcols(transcripts(txdb, columns = c("gene_id", "tx_name"))) - all_tx$gene_id <- as.character(all_tx$gene_id) - all_tx$is_coding <- all_tx$tx_name %in% protein_coding_tx - tmp <- splitAsList(all_tx$is_coding, all_tx$gene_id) - gene <- names(tmp) - n_tx <- lengths(tmp) - n_coding <- sum(tmp) - n_non_coding <- n_tx - n_coding - data.frame(gene, n_tx, n_coding, n_non_coding, stringsAsFactors = FALSE) - } - - # Get protein coding genes - geneid_codingdf <- summarizeProteinCodingGenes(txdb) - final_genes <- genes[genes$gene_id %in% geneid_codingdf[geneid_codingdf$n_coding > 0,]$gene] - - # Clean gene IDs - mcols(final_genes)$gene_id <- gsub("\\..*", "", mcols(final_genes)$gene_id) - - return(final_genes) -} - -##################################################################### -## STATISTICAL TESTING FUNCTION -##################################################################### - -func_ztest_gr_byrow <- function(gr, - x, - y, - # correction_method = "BH", - cohenthresh = 0.8) - { - ppval <- lapply(seq(nrow(as.data.frame(mcols(gr)))), function(i) { - - # Cohen's d calculation - cohend <- cohen.d( - unlist(as.vector(as.data.frame(mcols(gr))[x][i,])), - unlist(as.vector(as.data.frame(mcols(gr))[y][i,])) - ) - cohen.estimate <- cohend$estimate - cohend.magnitude <- as.character(cohend$magnitude) - - # BSDA Z-test implementation - # ztest <- z.test(x = as.data.frame(mcols(gr))[x][i,], - # y = as.data.frame(mcols(gr))[y][i,], - # sigma.x = sd(as.data.frame(mcols(gr))[x][i,]), - # sigma.y = sd(as.data.frame(mcols(gr))[y][i,]), - # alternative = 'two.sided', - # conf.level = 0.99) - - # ztest_pvalue <- ztest$p.value - - zzzz <- list(cohen.estimate, cohend.magnitude) - names(zzzz) <- c("cohen.estimate", "cohen.magnitude") - return(zzzz) - }) - - df_tomerge_mcols <- data.frame( - cohen.estimate = unlist(fun1(ppval, 1)), - cohen.magnitude = unlist(fun1(ppval, 2)) - ) - - # Apply Benjamini-Hochberg correction - # df_tomerge_mcols[[paste0("ztest_", correction_method, "_correct")]] <- p.adjust(df_tomerge_mcols$ztest, method = correction_method) - - mcols(gr) <- cbind(mcols(gr), df_tomerge_mcols) - - # Filter by Cohen's d - gr <- gr[abs(mcols(gr)$cohen.estimate) >= cohenthresh] - return(gr) -} - -##################################################################### -## MAIN BINS SELECTOR FUNCTION -##################################################################### - -Bins_selector <- function(combination, allmixeddf_grobj, fraction1 = "S2S", fraction2 = "S3", ths = 0.1) { - cat("Running Bins_selector for combination:", combination, "\n") - - # Get pr from global environment - pr <- get("pr", envir = .GlobalEnv) - - # Define groups to compare and select their samples - x <- get(pr[, combination][1]) - y <- get(pr[, combination][2]) - xgroup <- pr[, combination][1] - ygroup <- pr[, combination][2] - - cat("Groups:", paste(x, collapse = ", "), "vs", paste(y, collapse = ", "), "\n") - - new_selection <- allmixeddf_grobj - mcols(new_selection) <- mcols(new_selection)[c(x, y)] - - # Calculations - confint_mean_serr_first <- apply(as.matrix(mcols(new_selection)[c(x)]), 1, confidence_interval, nm = xgroup) - confint_mean_serr_second <- apply(as.matrix(mcols(new_selection)[c(y)]), 1, confidence_interval, nm = ygroup) - delta <- confint_mean_serr_first[1, ] - confint_mean_serr_second[1, ] - res <- cbind(t(confint_mean_serr_first), t(confint_mean_serr_second), delta) - df_toadd1 <- do.call("cbind", as.data.frame(res)) - mcols(new_selection) <- cbind(mcols(new_selection), df_toadd1) - - # Range analysis forward comparison - range_analysis <- lapply(1:length(y), function(n) { - y_name <- y[n] - z <- apply(as.matrix(mcols(new_selection)[c(paste0(xgroup, "_serrx2_lower"), paste0(xgroup, "_serrx2_upper"), paste0(xgroup, "_mean"), y_name)]), 1, - is_in_serrx2_range_and_shift, - xgroup_name = xgroup, - y_name = y_name - ) - return(as.data.frame(t(z))) - }) - - df_toadd <- do.call("cbind", range_analysis) - mcols(new_selection) <- cbind(mcols(new_selection), df_toadd) - - # Range analysis reverse comparison - range_analysis_rev <- lapply(1:length(x), function(n) { - x_name <- x[n] - z <- apply(as.matrix(mcols(new_selection)[c(paste0(ygroup, "_serrx2_lower"), paste0(ygroup, "_serrx2_upper"), paste0(ygroup, "_mean"), x_name)]), 1, - is_in_serrx2_range_and_shift, - xgroup_name = ygroup, - y_name = x_name - ) - return(as.data.frame(t(z))) - }) - - df_toadd2 <- do.call("cbind", range_analysis_rev) - mcols(new_selection) <- cbind(mcols(new_selection), df_toadd2) - - # Select bins - prvdf <- as.data.frame(new_selection) - prvdf_gr <- makeGRangesFromDataFrame(prvdf, keep.extra.columns = TRUE) - - # Calculate Cohen's d for all bins - all_bins_with_cohens <- func_ztest_gr_byrow(prvdf_gr, - x = x, - y = y, - cohenthresh = 0 # No filtering by Cohen's d here - ) - cat("Cohen's d calculated for", length(all_bins_with_cohens), "bins\n") - - prvdf_with_cohens <- as.data.frame(all_bins_with_cohens) - prvdftest <- prvdf_with_cohens[abs(prvdf_with_cohens[paste0(xgroup, "_mean")]) >= ths, ] - cat("After threshold filtering:", nrow(prvdftest), "bins remain\n") - - pprvlow <- prvdftest[paste0(y, "_shift")] == "lower" # <-- cambiato da _ov_specs a _shift - pprvhigh <- prvdftest[paste0(y, "_shift")] == "higher" # <-- cambiato da _ov_specs a _shift - - # Sum up by group 1 - prvdftest[, paste0(ygroup, "_sign_SUM")] <- rowSums(sapply(prvdftest[, paste0(y, "_sign")], as.numeric)) # cambiato da _ov_check a _sign - prvdftest$ovlow <- apply(pprvlow, 1, sum) * -1 - prvdftest$ovvhigh <- apply(pprvhigh, 1, sum) - - # Commutative group testing - pprvlow_X <- prvdftest[paste0(x, "_shift")] == "lower" # <-- cambiato da ov_specs a _shift - pprvhigh_X <- prvdftest[paste0(x, "_shift")] == "higher" # <-- cambiato da ov_specs a _shift - - # Sum up by group 2 - prvdftest[, paste0(xgroup, "_sign_SUM")] <- rowSums(sapply(prvdftest[, paste0(x, "_sign")], as.numeric)) ## <-- cambiato da _ov_check a _sign - - # Count characteristics - prvdftest$ovlow_X <- apply(pprvlow_X, 1, sum) * -1 - prvdftest$ovvhigh_X <- apply(pprvhigh_X, 1, sum) - - # Add meantosep column - prvdftest$meantosep <- prvdftest[,paste0(xgroup, "_mean_startsign")] - - prvdftest_gr <- makeGRangesFromDataFrame(prvdftest, keep.extra.columns = TRUE) - - # Apply statistical testing - up_down_to_ztest_grr <- prvdftest_gr[abs(mcols(prvdftest_gr)$cohen.estimate) >= 3] - - cat("Statistical testing completed. Regions passing threshold:", length(up_down_to_ztest_grr), "\n") - - # Separate bins according to meantosep - startmeanpos <- up_down_to_ztest_grr[up_down_to_ztest_grr$meantosep >= 0] - startmeanneg <- up_down_to_ztest_grr[up_down_to_ztest_grr$meantosep <= 0] - - # Select coherent bins based on direction of change - ovvhighconservedpos <- startmeanpos[abs(startmeanpos$ovvhigh) > abs(startmeanpos$ovlow)] - ovlowconservedpos <- startmeanpos[abs(startmeanpos$ovvhigh) < abs(startmeanpos$ovlow)] - ovvhighconservedneg <- startmeanneg[abs(startmeanneg$ovvhigh) > abs(startmeanneg$ovlow)] - ovlowconservedneg <- startmeanneg[abs(startmeanneg$ovvhigh) < abs(startmeanneg$ovlow)] - - # Add bin type labels - ovvhighconservedpos$bintype <- rep(paste0(fraction1, "_up"), length(ovvhighconservedpos)) - ovlowconservedpos$bintype <- rep(paste0(fraction1, "_down"), length(ovlowconservedpos)) - ovvhighconservedneg$bintype <- rep(paste0(fraction2, "_up"), length(ovvhighconservedneg)) - ovlowconservedneg$bintype <- rep(paste0(fraction2, "_down"), length(ovlowconservedneg)) - - # Control groups - controlstartmeanpos <- prvdftest_gr[prvdftest_gr$meantosep >= 0] - controlstartmeanneg <- prvdftest_gr[prvdftest_gr$meantosep <= 0] - - controlstartmeanpos$bintype <- rep(fraction1, length(controlstartmeanpos)) - controlstartmeanneg$bintype <- rep(fraction2, length(controlstartmeanneg)) - - controlstartmeanpos <- controlstartmeanpos[!controlstartmeanpos %in% - c(ovvhighconservedpos, ovlowconservedpos)] - controlstartmeanneg <- controlstartmeanneg[!controlstartmeanneg %in% - c(ovvhighconservedneg, ovlowconservedneg)] - - # All groups to return - all_gr_toreturn <- c(ovvhighconservedpos, - ovlowconservedpos, - ovvhighconservedneg, - ovlowconservedneg, - controlstartmeanpos, - controlstartmeanneg) - - # Make a list of bins to save and analyze - list_ofbins_to_save_and_analyse <- setNames( - list( - ovvhighconservedpos, - ovlowconservedpos, - ovvhighconservedneg, - ovlowconservedneg - ), - c( - paste0(fraction1, "_up"), - paste0(fraction1, "_down"), - paste0(fraction2, "_up"), - paste0(fraction2, "_down") - ) - ) - - # Gene analysis (if gtf provided) - list_of_vector_geneNumber <- list() - list_of_genes_vec <- list() - - final_genes_obj <- tryCatch(get("final_genes", envir = .GlobalEnv), error = function(e) NULL) - - if (!is.null(final_genes_obj)) { - for (i in names(list_ofbins_to_save_and_analyse)) { - gr_touse <- list_ofbins_to_save_and_analyse[[i]] - if (length(gr_touse) != 0) { - cat(paste0(ygroup, "_vs_", xgroup, "_", i, " has ", length(gr_touse), " regions"), "\n") - - # Calculate overlapping genes - genes_gr <- final_genes_obj[findOverlaps(gr_touse, promoters(final_genes_obj, upstream = 2500, downstream = 500))@to] - - list_of_vector_geneNumber[[i]] <- length(mcols(genes_gr)$gene_id) - list_of_genes_vec[[paste0(ygroup, "_vs_", xgroup, "_", i)]] <- unique(mcols(genes_gr)$gene_id) - } else { - list_of_vector_geneNumber[[i]] <- 0 - list_of_genes_vec[[paste0(ygroup, "_vs_", xgroup, "_", i)]] <- character(0) - } - } - } else { - for (i in names(list_ofbins_to_save_and_analyse)) { - list_of_vector_geneNumber[[i]] <- 0 - list_of_genes_vec[[paste0(ygroup, "_vs_", xgroup, "_", i)]] <- character(0) - } - } - - # Numeric coding for the groups - mcols(ovvhighconservedpos)[[paste0(ygroup, "_vs_", xgroup)]] <- rep(2, length(ovvhighconservedpos)) - mcols(ovlowconservedpos)[[paste0(ygroup, "_vs_", xgroup)]] <- rep(1, length(ovlowconservedpos)) - mcols(ovvhighconservedneg)[[paste0(ygroup, "_vs_", xgroup)]] <- rep(-1, length(ovvhighconservedneg)) - mcols(ovlowconservedneg)[[paste0(ygroup, "_vs_", xgroup)]] <- rep(-2, length(ovlowconservedneg)) - - # Prepare results list - x <- list() - x[[paste0(ygroup, "_vs_", xgroup, "_all_shifting_bins")]] <- c( - ovvhighconservedpos, - ovlowconservedpos, - ovvhighconservedneg, - ovlowconservedneg - ) - x[[names(list_ofbins_to_save_and_analyse[1])]] <- list_ofbins_to_save_and_analyse[[1]] - x[[names(list_ofbins_to_save_and_analyse[2])]] <- list_ofbins_to_save_and_analyse[[2]] - x[[names(list_ofbins_to_save_and_analyse[3])]] <- list_ofbins_to_save_and_analyse[[3]] - x[[names(list_ofbins_to_save_and_analyse[4])]] <- list_ofbins_to_save_and_analyse[[4]] - x[[paste0(ygroup, "_genes_", xgroup)]] <- list_of_vector_geneNumber - x[["genes"]] <- list_of_genes_vec - x[[paste0(ygroup, "_gr_", xgroup)]] <- prvdftest_gr - x[[paste0(ygroup, "_allgr_", xgroup)]] <- all_bins_with_cohens - - names(x) <- c( - paste0(ygroup, "_vs_", xgroup, "_all_shifting_bins"), - paste0(ygroup, "_", names(list_ofbins_to_save_and_analyse[1]), "_", xgroup), - paste0(ygroup, "_", names(list_ofbins_to_save_and_analyse[2]), "_", xgroup), - paste0(ygroup, "_", names(list_ofbins_to_save_and_analyse[3]), "_", xgroup), - paste0(ygroup, "_", names(list_ofbins_to_save_and_analyse[4]), "_", xgroup), - paste0(ygroup, "_genes_", xgroup), - "genes", - paste0(ygroup, "_gr_", xgroup), - paste0(ygroup, "_allgr_", xgroup) - ) - - cat("Results:", paste0(names(x[paste0(ygroup, "_vs_", xgroup, "_all_shifting_bins")]), "_", length(x[[paste0(ygroup, "_vs_", xgroup, "_all_shifting_bins")]]), "_bins"), "\n") - return(x) -} diff --git a/conf/modules.config b/conf/modules.config index 48519511..45399138 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -375,4 +375,38 @@ process { mode: params.publish_dir_mode, ] } + + withName: ".*CHR_SPLIT.*" { + publishDir = [ + path: { "${params.outdir}/compartments" }, + mode: params.publish_dir_mode, + enabled: false + ] + } + + withName: ".*CHR_COMPARTMENTS_CALLING.*" { + publishDir = [ + path: { "${params.outdir}/compartments" }, + mode: params.publish_dir_mode, + enabled: false + ] + } + + withName: ".*CHR_COMBINE_COMPARTMENTS.*" { + publishDir = [ + path: { "${params.outdir}/compartments" }, + mode: params.publish_dir_mode, + pattern: '*.{bed,bedgraph}', + enabled: true + ] + } + + withName: ".*BUILD_CONSENSUS.*" { + publishDir = [ + path: { "${params.outdir}/compartments/consensus" }, + mode: params.publish_dir_mode, + pattern: '*.bed', + enabled: true + ] + } } diff --git a/conf/test_groups.config b/conf/test_groups.config new file mode 100644 index 00000000..06c63eaf --- /dev/null +++ b/conf/test_groups.config @@ -0,0 +1,43 @@ +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Nextflow config file for running minimal tests +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Defines input files and everything required to run a fast and simple pipeline test. + + Use as follows: + nextflow run nf-core/sammyseq -profile test, --outdir + +---------------------------------------------------------------------------------------- +*/ + +process { + resourceLimits = [ + cpus: 4, + memory: '2.GB', + time: '1.h' + ] +} + +params { + config_profile_name = 'Groups Test profile' + config_profile_description = 'Minimal test dataset to check pipeline function with a complete analysis' + + // Input data + input = 'https://genome.isasi.cnr.it/biocomp/test-datasets/mouse_MEFACKOvsWT/test_mm_chr4chr8_ss.csv' + + // Analysis details + normalizeUsing = 'RPKM' + bw_resolution = 50 + extendReads = 250 + comparison_maker = 'spp' + comparison = 'S2SvsS3' + compartmentalization_analysis = true + differential_solubility = true + compare_groups = 'MEFACKOvsMEFwt' + binsize = 150000 + + // Genome references + fasta = 'https://genome.isasi.cnr.it/biocomp/test-datasets/refs/mm10/chr4chr8.fa.gz' + blacklist = 'https://raw.githubusercontent.com/daisymut/test-datasets/e774a4b965c23ac599b9d663f92cd1328ec83f0f/testdata/mm10-blacklist.v2.chr4chr8.bed' + gtf = 'https://raw.githubusercontent.com/daisymut/test-datasets/e774a4b965c23ac599b9d663f92cd1328ec83f0f/testdata/gencode.vM25.basic.annotation.chr4chr8.gtf' + } diff --git a/docs/images/nf-core-sammyseq_tubemap.png b/docs/images/nf-core-sammyseq_tubemap.png index cc70d35e..0555b078 100644 Binary files a/docs/images/nf-core-sammyseq_tubemap.png and b/docs/images/nf-core-sammyseq_tubemap.png differ diff --git a/docs/output.md b/docs/output.md index 3ccabb14..40cec8a5 100644 --- a/docs/output.md +++ b/docs/output.md @@ -17,9 +17,10 @@ The pipeline is built using [Nextflow](https://www.nextflow.io/) and processes d - [Samtools reads filtering](#samtools-reads-filtering) - [Signal track generation](#signal-track-generation) - [DeepTools based QC](#deeptools-based-qc) +- [Compartmentalization Analysis](#compartmentalization-analysis) - [Comparisons](#comparisons) -- [MultiQC](#multiqc) - [Differential Solubility Analysis](#differential-solubility-analysis) +- [MultiQC](#multiqc) - [Pipeline information](#pipeline-information) ### Read quality check @@ -175,6 +176,21 @@ This analysis uses computeMatrix in reference-point mode to generate coverage pr +### Compartmentalization Analysis + +When `--compartmentalization_analysis` is enabled, a module is triggered to infer A/B chromatin compartments from SAMMY-seq signal tracks. The analysis is based on fixed size genomic binning, performed per chromosome using bedtools makewindows, which divides the genome into windows of a defined size (default 50000) set by the `--binsize` parameter. To restrict the analysis to specific chromosomes, a BED file with only the chromosomes to include can be provided via the `--keep_regions_bed` parameter. The `--gtf` parameter is also required, as gene annotations are used in downstream steps. Each fraction is identified using the `experimentalID` column in the samplesheet allowing fractions of the same sample to be analyzed together. The Compartments calling is based on CALDER2 algorithm, which builds a correlation matrix across genomic bins and applies eigenvector decomposition to classify each bin as either A (open) or B (closed) compartment. After compartment calling, results from all analyzed chromosomes are combined into a single compartment BED and a single BedGraph file with eigenvalues for each sample. Consensus compartment profiles are then generated across biological replicates within each `sample_group`. A majority consensus assigns the most frequent compartment call (A or B) across replicates to each genomic bin, marking ties as NA, while a strict consensus annotates only bins where all replicates unanimously agree on the same compartment call, marking all other bins as NA. + +

Output files + +- `compartments/` + - `_combined_compartments.bed`: BED file with genomic bins annotated as A or B compartments for each sample. Results are combined from all analyzed chromosomes. + - `_combined_compartments.bedGraph`: BedGraph file with PC1 eigenvector values for each genomic bin. +- `compartments/consensus/` + - `_consensus_majority.bed`: BED file with majority consensus compartment calls across replicates + - `_consensus_strict.bed`: BED file with strict consensus compartment calls (unanimous agreement only) + +
+ ### Comparisons Pairwise comparisons can be generated by setting the parameter `--comparison_maker spp` (default) either by `--comparison` or `--comparison_file`. The difference between each fraction read density profile, smoothed by the Gaussian kernel is calculated and saved in bigwig format, as described in Kharchenko PK, Tolstorukov MY, Park PJ "Design and analysis of ChIP-seq experiments for DNA-binding proteins" Nat. Biotech. doi:10.1038/nbt.1508 diff --git a/docs/usage.md b/docs/usage.md index 148c4bf5..76633f44 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -96,6 +96,14 @@ CTRL004_S4,/home/sammy/test_data/CTRL004_S4_chr22only.fq.gz,,CTRL004,S4,CTRL | `sample_group` | Identifier used to group samples that belong to the same biological condition condition. | | | +### Compartmentalization Analysis + +The compartmentalization analysis can be enabled by setting the `--compartmentalization_analysis` parameter and requires the `--gtf` parameter for gene annotations. + +The signal profiles from all fractions of the same sample (as defined by the `experimentalID` field in the samplesheet) are integrated to identify chromatin compartments. The genome is divided into fixed-size windows (controlled by the `--binsize` parameter, default: 50000 bp), and the signal from all fractions in each window is combined into a correlation matrix. By applying principal component analysis (PCA), the first eigenvector (PC1) of this matrix is used to classify genomic regions as belonging to the active (A) or inactive (B) compartment, following strategies analogous to Hi-C analysis. This enables the detection of large-scale chromatin reorganizations, such as regions switching between active and inactive states across different biological conditions. + +The analysis is performed chromosome-by-chromosome using the [`CALDER2`](https://github.com/CSOgroup/CALDER2) algorithm, and results are then combined into genome-wide compartment annotations for each sample. The analysis is performed chromosome-by-chromosome using the CALDER2 algorithm, and results are then combined into genome-wide compartment annotations for each sample. Consensus profiles are then generated across biological replicates within each `sample_group`. A majority consensus assigns the most frequent compartment call across replicates to each bin (ties are marked as NA), while a strict consensus annotates only bins where all replicates unanimously agree (all other bins are marked as NA). + ### Pairwise comparisons The pipeline offers two different methods for generating these comparisons, selected with the `--comparison_maker` parameter. @@ -133,6 +141,8 @@ For 4f-SAMMYseq protocols (S2S, S2L, S3, S4), valid comparisons are: S2SvsS4 - Compare S2S fraction vs S4 fraction S2LvsS4 - Compare S2L fraction vs S4 fraction S4vsS3 - Compare S4 fraction vs S3 fraction + S3vsS2S - Compare S3 fraction vs S2S fraction + S3vsS2L - Compare S3 fraction vs S2L fraction > [!NOTE] > For 3f-SAMMYseq protocols (S2, S3, S4), valid comparisons are: @@ -140,6 +150,7 @@ For 4f-SAMMYseq protocols (S2S, S2L, S3, S4), valid comparisons are: S2vsS3 - Compare S2 fraction vs S3 fraction S2vsS4 - Compare S2 fraction vs S4 fraction S4vsS3 - Compare S4 fraction vs S3 fraction + S3vsS2 - Compare S3 fraction vs S2 fraction The pipeline will automatically create comparisons only between fractions from the same `experimentalID` (biological replicate), ensuring that comparisons are made within the same experimental condition rather than across different replicates. diff --git a/modules/local/build_consensus/main.nf b/modules/local/build_consensus/main.nf new file mode 100644 index 00000000..0035b5e6 --- /dev/null +++ b/modules/local/build_consensus/main.nf @@ -0,0 +1,53 @@ +process BUILD_CONSENSUS { + tag "$sample_group consensus - ${bed_files.size()} replicates" + label 'process_low' + + input: + tuple val(sample_group), path(bed_files) + + output: + path("${sample_group}_consensus_majority.bed"), emit: consensus_majority + path("${sample_group}_consensus_strict.bed"), emit: consensus_strict + + script: + def num_samples = bed_files.size() + """ + # MAJORITY: more A or more B + echo "track name=\\"${sample_group}_majority\\"" > ${sample_group}_consensus_majority.bed + cat ${bed_files.join(' ')} | grep -v '^track' | sort -k1,1V -k2,2n | \\ + awk '{ + bin = \$1 "\\t" \$2 "\\t" \$3 + fields[bin] = \$5 "\\t" \$6 "\\t" \$7 "\\t" \$8 + if (\$4 == "A") a[bin]++ + else if (\$4 == "B") b[bin]++ + } + END { + for (bin in fields) { + if (a[bin] > b[bin]) comp = "A" + else if (b[bin] > a[bin]) comp = "B" + else comp = "NA" + color = (comp == "A") ? "90,149,143" : (comp == "B") ? "224,170,88" : "255,255,255" + print bin "\\t" comp "\\t" fields[bin] "\\t" color + } + }' | sort -k1,1V -k2,2n >> ${sample_group}_consensus_majority.bed + + # STRICT: all A or all B + echo "track name=\\"${sample_group}_strict\\"" > ${sample_group}_consensus_strict.bed + cat ${bed_files.join(' ')} | grep -v '^track' | sort -k1,1V -k2,2n | \\ + awk -v n=${num_samples} '{ + bin = \$1 "\\t" \$2 "\\t" \$3 + fields[bin] = \$5 "\\t" \$6 "\\t" \$7 "\\t" \$8 + if (\$4 == "A") a[bin]++ + else if (\$4 == "B") b[bin]++ + } + END { + for (bin in fields) { + if (a[bin] == n) comp = "A" + else if (b[bin] == n) comp = "B" + else comp = "NA" + color = (comp == "A") ? "90,149,143" : (comp == "B") ? "224,170,88" : "255,255,255" + print bin "\\t" comp "\\t" fields[bin] "\\t" color + } + }' | sort -k1,1V -k2,2n >> ${sample_group}_consensus_strict.bed + """ +} diff --git a/modules/local/build_consensus/meta.yml b/modules/local/build_consensus/meta.yml new file mode 100644 index 00000000..36fc3be5 --- /dev/null +++ b/modules/local/build_consensus/meta.yml @@ -0,0 +1,44 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "build_consensus" +description: Generate majority and strict consensus compartment calls from multiple samples +keywords: + - consensus + - compartments + - chromatin + - genomics +tools: + - awk: + description: Pattern scanning and processing language + homepage: https://www.gnu.org/software/gawk/ + documentation: https://www.gnu.org/software/gawk/manual/gawk.html + licence: ["GPL-3.0-or-later"] + - coreutils: + description: GNU core utilities + homepage: https://www.gnu.org/software/coreutils/ + documentation: https://www.gnu.org/software/coreutils/manual/ + licence: ["GPL-3.0-or-later"] + +input: + - sample_group: + type: string + description: Sample group identifier for consensus calling + pattern: "^[a-zA-Z0-9_-]+$" + - bed_files: + type: file + description: List of BED files with compartment calls (A/B) to generate consensus from + pattern: "*.bed" + +output: + - consensus_majority: + type: file + description: BED file with majority consensus (more A or more B across samples) + pattern: "*_consensus_majority.bed" + - consensus_strict: + type: file + description: BED file with strict consensus (all A or all B across samples) + pattern: "*_consensus_strict.bed" + +authors: + - "@ugoiannacchero" +maintainers: + - "@daisymut" diff --git a/modules/local/chr_combine_compartments/main.nf b/modules/local/chr_combine_compartments/main.nf new file mode 100644 index 00000000..c3a09a4f --- /dev/null +++ b/modules/local/chr_combine_compartments/main.nf @@ -0,0 +1,45 @@ +process CHR_COMBINE_COMPARTMENTS { + tag "$sample_id" + label 'process_medium' + + conda "conda-forge::coreutils=9.1" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/ubuntu:20.04' : + 'nf-core/ubuntu:20.04' }" + + input: + tuple val(sample_id), path(bed_files), path(bedgraph_files) + + output: + tuple val(sample_id), path("${sample_id}_combined_compartments.bed"), emit: combined_beds + tuple val(sample_id), path("${sample_id}_combined_compartments_eigen.bedgraph"), emit: combined_bedgraphs + path "versions.yml", emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + """ + # Combine BED: pull together all the bed files and sort them by chromosome and start position. The grep commands remove any existing track or comment lines from the input files before sorting. + { echo 'track name="${sample_id}" description="${sample_id} (Emission ordered)" visibility=1 itemRgb="On"'; cat ${bed_files.join(' ')} | grep -v '^track' | grep -v '^#' | sort -k1,1V -k2,2n; } > ${sample_id}_combined_compartments.bed + + # Combine BedGraph:pull together all the bedgraph files and sort them by chromosome and start position. The grep commands remove any existing track, browser, or comment lines from the input files before sorting. + { echo 'track type=bedGraph name="${sample_id}_eigenvector" description="${sample_id} eigenvector" visibility=full color=90,149,143 altColor=224,170,88 priority=20'; cat ${bedgraph_files.join(' ')} | grep -v '^track' | grep -v '^browser' | grep -v '^#' | sort -k1,1V -k2,2n; } > ${sample_id}_combined_compartments_eigen.bedgraph + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + coreutils: \$(sort --version | head -n1 | sed 's/^.* //g') + END_VERSIONS + """ + + stub: + """ + touch ${sample_id}_combined_compartments.bed + touch ${sample_id}_combined_compartments_eigen.bedgraph + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + coreutils: \$(sort --version | head -n1 | sed 's/^.* //g') + END_VERSIONS + """ +} diff --git a/modules/local/chr_combine_compartments/meta.yml b/modules/local/chr_combine_compartments/meta.yml new file mode 100644 index 00000000..9a6fecd5 --- /dev/null +++ b/modules/local/chr_combine_compartments/meta.yml @@ -0,0 +1,52 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/yaml-schema.json +name: "chr_combine_compartments" +description: Combine per-chromosome compartment BED and BedGraph files into genome-wide files per sample +keywords: + - compartmentalization + - merge + - bed + - bedgraph + - genomics +tools: + - coreutils: + description: GNU core utilities for basic file, shell and text manipulation + homepage: https://www.gnu.org/software/coreutils/ + documentation: https://www.gnu.org/software/coreutils/manual/ + tool_dev_url: https://git.savannah.gnu.org/cgit/coreutils.git + licence: ["GPL-3.0-or-later"] + +input: + - sample_id: + type: string + description: Sample identifier + - bed_files: + type: file + description: List of per-chromosome compartment BED files + pattern: "*_compartments.bed" + - bedgraph_files: + type: file + description: List of per-chromosome eigenvector BedGraph files + pattern: "*_comp_eigenvector.bedgraph" + +output: + - combined_beds: + type: file + description: | + Combined BED file with all chromosomes for the sample + Structure: [ sample_id, bed ] + pattern: "*_combined_compartments.bed" + - combined_bedgraphs: + type: file + description: | + Combined BedGraph file with all chromosomes for the sample + Structure: [ sample_id, bedgraph ] + pattern: "*_combined_compartments_eigen.bedgraph" + - versions: + type: file + description: File containing software versions + pattern: "versions.yml" + +authors: + - "@ugoiannacchero" +maintainers: + - "@daisymut" diff --git a/modules/local/chr_compartments_calling/environment.yml b/modules/local/chr_compartments_calling/environment.yml new file mode 100644 index 00000000..6fca0b64 --- /dev/null +++ b/modules/local/chr_compartments_calling/environment.yml @@ -0,0 +1,20 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json +channels: + - conda-forge + - bioconda + - r +dependencies: + - r-base=4.3.3 + - r-calder2=0.7 + - r-patchwork=1.3.2 + - bioconductor-gviz=1.46.1 + - bioconductor-limma=3.58.1 + - bioconductor-rtracklayer=1.62.0 + - bioconductor-genomicranges=1.54.1 + - bioconductor-genomeinfodb=1.38.1 + - bioconductor-genomicfeatures=1.54.1 + - r-data.table=1.17.8 + - r-dplyr=1.1.4 + - r-effsize=0.8.1 + - r-bsda=1.2.2 + - procps-ng diff --git a/modules/local/chr_compartments_calling/main.nf b/modules/local/chr_compartments_calling/main.nf new file mode 100644 index 00000000..c32de116 --- /dev/null +++ b/modules/local/chr_compartments_calling/main.nf @@ -0,0 +1,21 @@ +process CHR_COMPARTMENTS_CALLING { + tag "${patient}_${meta.chromosome}" + label 'process_medium' + + container 'ghcr.io/daisymut/sammyr:0.0.0.9001' + // conda "${moduleDir}/environment.yml" + + input: + tuple val(meta), path(chr_bed), val(patient), path(csv) + val binsize + path gtf + + output: + tuple val(patient), path("*.bed") , emit: bed_files + tuple val(patient), path("*.bedgraph") , emit: bedgraph_files + path "versions.yml" , emit: versions + + script: + def args = task.ext.args ?: '' + template 'chr_compartments_calling.R' +} diff --git a/modules/local/chr_compartments_calling/templates/chr_compartments_calling.R b/modules/local/chr_compartments_calling/templates/chr_compartments_calling.R new file mode 100644 index 00000000..c1935241 --- /dev/null +++ b/modules/local/chr_compartments_calling/templates/chr_compartments_calling.R @@ -0,0 +1,56 @@ +#!/usr/bin/env Rscript + +suppressMessages({ + library(parallel) + library(data.table) + library(GenomicRanges) + library(rtracklayer) + library(patchwork) + library(Gviz) + library(CALDER) + library(sammyR) +}) + +options(ucscChromosomeNames=FALSE) + +## Get parameters from Nextflow template +patient <- "${patient}" +chromosome <- "${meta.chromosome}" +binsize <- as.integer(${binsize}) +input_file <- "${csv}" +chrom_bed <- "${chr_bed}" +gene_gtf <- "${gtf}" + +## Load data +comp_df <- fread(input_file, data.table = FALSE) +comp_df_replica <- comp_df[comp_df\$Patient_name == patient, ] + +bins_gr <- import(chrom_bed, format = "BED") +genes_gr <- import(gene_gtf, format = "GTF") + +sub2_colors <- c("B" = "#4575b4", "A" = "#d73027") +subs_file <- paste0(patient, "_compartment___", chromosome, '_', binsize, ".Rdata") + +# Run SAMMY +sub_objs <- call_subcompartments_sammy( + patients = patient, + tracks_db = comp_df_replica, + bins_gr = bins_gr, + subs_file = subs_file, + binsize = binsize, + chr = chromosome, + genes_gr = genes_gr, + keeping_bins1 = "all", + sublevel = "sub.2", + sub_colors = sub2_colors +) + +generate_files(sub_objs, chromosome) + +# Write versions +pkgs <- c("sammyR", "CALDER", "GenomicRanges", "rtracklayer", "Gviz", "data.table") +writeLines(c( + '"${task.process}":', + paste0(' r-base: "', paste(R.Version()[c("major","minor")], collapse="."), '"'), + vapply(pkgs, function(p) paste0(' ', p, ': ', as.character(packageVersion(p))), character(1)) +), "versions.yml") diff --git a/modules/local/chr_split/main.nf b/modules/local/chr_split/main.nf new file mode 100644 index 00000000..c7786dbf --- /dev/null +++ b/modules/local/chr_split/main.nf @@ -0,0 +1,19 @@ +process CHR_SPLIT { + tag "$meta.id" + label 'process_low' + + input: + tuple val(meta), path(bed) + path chrom_sizes + + output: + tuple val(meta), path("*.bed"), emit: beds + + script: + """ + awk ' + NR==FNR { chroms[\$1]; next } + \$1 in chroms { print > \$1".bed" } + ' ${chrom_sizes} ${bed} + """ +} diff --git a/modules/local/chr_split/meta.yml b/modules/local/chr_split/meta.yml new file mode 100644 index 00000000..212258da --- /dev/null +++ b/modules/local/chr_split/meta.yml @@ -0,0 +1,46 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/meta-schema.json +name: "chr_split" +description: Splits a BED file into separate files per chromosome based on a chromosome sizes file +keywords: + - bed + - chromosome + - split + - filter +tools: + - "awk": + description: | + "Pattern scanning and processing language for splitting BED files by chromosome" + homepage: "https://www.gnu.org/software/gawk/" + licence: ["GPL-3.0"] +input: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. `[ id:'sample1' ]` + - bed: + type: file + description: Input BED file to split by chromosome + pattern: "*.bed" + ontologies: + - edam: "http://edamontology.org/format_3003" + - chrom_sizes: + type: file + description: File containing chromosome names and sizes (tab-separated) + pattern: "*.{txt,sizes,chrom.sizes}" +output: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. `[ id:'sample1' ]` + - beds: + type: file + description: BED files split by chromosome (one file per chromosome) + pattern: "*.bed" + ontologies: + - edam: "http://edamontology.org/format_3003" +authors: + - "@ugoiannacchero" +maintainers: + - "@daisymut" diff --git a/modules/local/differential_enrichment/environment.yml b/modules/local/differential_enrichment/environment.yml index c811ba82..6fca0b64 100644 --- a/modules/local/differential_enrichment/environment.yml +++ b/modules/local/differential_enrichment/environment.yml @@ -1,17 +1,20 @@ ---- # yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/modules/environment-schema.json channels: - conda-forge - bioconda - r dependencies: - - r-base=4.1.3 - - r-data.table=1.14.8 - - bioconductor-genomeinfodb=1.30.1 - - bioconductor-genomicranges=1.46.1 - - bioconductor-preprocesscore=1.56.0 - - bioconductor-rtracklayer=1.54.0 + - r-base=4.3.3 + - r-calder2=0.7 + - r-patchwork=1.3.2 + - bioconductor-gviz=1.46.1 + - bioconductor-limma=3.58.1 + - bioconductor-rtracklayer=1.62.0 + - bioconductor-genomicranges=1.54.1 + - bioconductor-genomeinfodb=1.38.1 + - bioconductor-genomicfeatures=1.54.1 + - r-data.table=1.17.8 + - r-dplyr=1.1.4 - r-effsize=0.8.1 - - r-dplyr=1.1.2 - - bioconductor-genomicfeatures=1.46.1 - - r-bsda=1.2.1 + - r-bsda=1.2.2 + - procps-ng diff --git a/modules/local/differential_enrichment/main.nf b/modules/local/differential_enrichment/main.nf index f9fb8e6e..200f2b9e 100644 --- a/modules/local/differential_enrichment/main.nf +++ b/modules/local/differential_enrichment/main.nf @@ -1,7 +1,7 @@ process DIFFERENTIAL_ENRICHMENT { tag "${contrast_data[0]}vs${contrast_data[1]}" - container 'docker.io/ciuki97/differential-solubility-analysis:v0.0.1' - conda "${moduleDir}/environment.yml" + container 'ghcr.io/daisymut/sammyr:0.0.0.9001' + // conda "${moduleDir}/environment.yml" label 'process_high' input: diff --git a/modules/local/differential_enrichment/templates/differential_enrichment.R b/modules/local/differential_enrichment/templates/differential_enrichment.R index b7821135..48742d32 100644 --- a/modules/local/differential_enrichment/templates/differential_enrichment.R +++ b/modules/local/differential_enrichment/templates/differential_enrichment.R @@ -3,17 +3,10 @@ suppressMessages({ library(data.table) library(rtracklayer) - library(GenomicRanges) - library(GenomeInfoDb) - library(preprocessCore) - library(effsize) - library(dplyr) - library(GenomicFeatures) - library(BSDA) + library(sammyR) + library(limma) }) -source("${projectDir}/bin/differential_solubility_functions.R") - ################################################ ################################################ ## PARSE PARAMETERS FROM NEXTFLOW ## @@ -62,7 +55,7 @@ id_col <- which(colnames(comp_db) == 'experimental_id') group_col <- which(colnames(comp_db) == 'sample_group') ################################################ -## FUNCTIONS ## +## SETUP ## ################################################ # Check if specified groups exist in the data @@ -75,49 +68,6 @@ if (length(missing_groups) > 0) { "Available groups: [", paste(available_groups, collapse = ", "), "]\\n") } -# metadata addition function -add_metadata <- function(df, comparison_name, current_ratio, fraction, direction) { - df[['comparison']] <- comparison_name - df[['ratio']] <- current_ratio - df[['fraction']] <- fraction - df[['direction']] <- direction - return(df) -} - -# save bins data to CSV in a structured way -save_bins_data <- function(data_list, current_ratio, comparison_name, file_suffix, g1 = NULL, g2 = NULL) { - if (length(data_list)) { - df <- do.call(rbind, data_list) - if (grepl("all_bins", file_suffix, ignore.case = TRUE)) { - cols_to_remove <- c('fraction', 'direction') - available_cols <- colnames(df) - cols_to_keep <- setdiff(available_cols, cols_to_remove) - df <- df[, cols_to_keep] - } else { - base_cols <- c('seqnames', 'start', 'end', 'ratio', 'comparison', 'fraction', 'direction') - essential_stats_cols <- c() - if (!is.null(g1) && !is.null(g2)) { - essential_stats_cols <- c( - paste0(g1, "_serrx2_lower"), paste0(g1, "_serrx2_upper"), - paste0(g1, "_mean"), paste0(g1, "_serrX2"), - paste0(g2, "_serrx2_lower"), paste0(g2, "_serrx2_upper"), - paste0(g2, "_mean"), paste0(g2, "_serrX2"), - "delta", - "cohen.estimate", - "cohen.magnitude" - )} - essential_cols <- c(base_cols, essential_stats_cols) - available_cols <- colnames(df) - cols_to_keep <- intersect(essential_cols, available_cols) - df <- df[, cols_to_keep] - } - output_file <- paste0(current_ratio, "_", comparison_name, "_", file_suffix, ".csv") - write.csv(df, file = output_file, quote = FALSE, row.names = FALSE) - } else { - cat("No", file_suffix, "data found for", comparison_name, "\\n") - } -} - # Setup gene annotation if (gtf_file != "" && file.exists(gtf_file)) { tryCatch({ @@ -162,7 +112,7 @@ for (current_ratio in selected_ratios) { gr1 <- GenomicRanges::makeGRangesFromDataFrame(bindf, keep.extra.columns = TRUE) gr1_names <- names(gr1@elementMetadata) - S4Vectors::mcols(gr1) <- preprocessCore::normalize.quantiles(as.matrix(S4Vectors::mcols(gr1))) + S4Vectors::mcols(gr1) <- limma::normalizeQuantiles(as.matrix(S4Vectors::mcols(gr1))) names(gr1@elementMetadata) <- gr1_names allmixeddf_grobj <- GenomicRanges::sort(gr1) unique_groups <- unique(Sample_groups) @@ -310,10 +260,10 @@ writeLines( paste(' bioconductor-rtracklayer:', as.character(packageVersion('rtracklayer'))), paste(' bioconductor-genomicranges:', as.character(packageVersion('GenomicRanges'))), paste(' bioconductor-genomeinfodb:', as.character(packageVersion('GenomeInfoDb'))), - paste(' bioconductor-preprocesscore:', as.character(packageVersion('preprocessCore'))), + paste(' bioconductor-limma:', as.character(packageVersion('limma'))), paste(' r-effsize:', as.character(packageVersion('effsize'))), - paste(' r-dplyr:', as.character(packageVersion('dplyr'))), + # paste(' r-dplyr:', as.character(packageVersion('dplyr'))), paste(' bioconductor-genomicfeatures:', as.character(packageVersion('GenomicFeatures'))), - paste(' r-bsda:', as.character(packageVersion('BSDA'))) + paste(' r-sammyr:', as.character(packageVersion('sammyR'))) ), 'versions.yml') diff --git a/nextflow.config b/nextflow.config index 95034a5a..76361bca 100644 --- a/nextflow.config +++ b/nextflow.config @@ -10,13 +10,15 @@ params { // Input options - input = null - comparison = null - comparison_file = null - comparison_maker = 'spp' - differential_solubility = false - solubility_threshold = 0.1 - compare_groups = null + input = null + comparison = null + comparison_file = null + comparison_maker = 'spp' + differential_solubility = false + solubility_threshold = 0.1 + compare_groups = null + compartmentalization_analysis = false + skip_deeptools_qc = false // References genome = null @@ -59,7 +61,7 @@ params { binsize = 50000 fragment_size = 100 // PARAM NOT USED (see module deeptools/plotfingerprint) plotfingerprint = false - bigwigcompare_pseudocount = 1e-14 + bigwigcompare_pseudocount = 1 bigwigcompare_operation = 'log2' bigwigcompare_skip_non_covered_regions = true bigwigcompare_fixed_step = false @@ -220,6 +222,7 @@ profiles { } test { includeConfig 'conf/test.config' } test_full { includeConfig 'conf/test_full.config' } + test_groups { includeConfig 'conf/test_groups.config' } test_comparison { includeConfig 'conf/test_comparison.config' } } diff --git a/nextflow_schema.json b/nextflow_schema.json index 530dd5c5..55c847bb 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -98,8 +98,8 @@ "comparison": { "type": "string", "fa_icon": "fas fa-percentage", - "description": "Fraction comparisons per experimentalID (e.g., \"S2SvsS3\" or \"S2SvsS3,S2LvsS3\")", - "pattern": "^(S2SvsS3|S2LvsS3|S2SvsS4|S2LvsS4|S2vsS3|S2vsS4|S4vsS3)(,(S2SvsS3|S2LvsS3|S2SvsS4|S2LvsS4|S2vsS3|S2vsS4|S4vsS3))*$" + "description": "Fraction comparisons per experimentalID (e.g., \"S2SvsS3\", \"S3vsS2S\")", + "pattern": "^(S2SvsS3|S2LvsS3|S2SvsS4|S2LvsS4|S2vsS3|S2vsS4|S4vsS3|S3vsS2|S3vsS2S|S3vsS2L)(,(S2SvsS3|S2LvsS3|S2SvsS4|S2LvsS4|S2vsS3|S2vsS4|S4vsS3|S3vsS2|S3vsS2S|S3vsS2L))*$" }, "comparison_file": { "type": "string", @@ -134,6 +134,10 @@ "description": "Define comparison order for sample groups in differential solubility analysis using sample_group names from input samplesheet (format: 'GroupAvsGroupB,GroupAvsGroupC,GroupBvsGroupC'). Default: all groups in samplesheet order.", "errorMessage": "Must follow format 'GroupAvsGroupB' or 'GroupAvsGroupB,GroupAvsGroupC,GroupBvsGroupC'. Group names must match sample_group column in input samplesheet." }, + "compartmentalization_analysis": { + "type": "boolean", + "description": "When set to true, this parameter activates the compartmentalization analysis workflow" + }, "stopAt": { "type": "string", "description": "Specify after which step the pipeline should stop.", @@ -326,7 +330,7 @@ "bigwigcompare_pseudocount": { "type": "number", "description": "A small number to avoid x/0. Only useful together with \u2013operation log2 or \u2013operation ratio.", - "default": 1e-14, + "default": 1, "fa_icon": "fas fa-chart-bar" }, "bigwigcompare_operation": { @@ -453,6 +457,10 @@ "help_text": "Incoming hook URL for messaging service. Currently, MS Teams and Slack are supported.", "hidden": true }, + "skip_deeptools_qc": { + "type": "boolean", + "description": "Skip deeptools correlation analysis between all tracks" + }, "skip_multiqc": { "type": "boolean", "fa_icon": "fas fa-fast-forward", diff --git a/ro-crate-metadata.json b/ro-crate-metadata.json index 169b821c..b466d22d 100644 --- a/ro-crate-metadata.json +++ b/ro-crate-metadata.json @@ -23,7 +23,7 @@ "@type": "Dataset", "creativeWorkStatus": "InProgress", "datePublished": "2025-11-06T09:32:18+00:00", - "description": "

\n \n \n \"nf-core/sammyseq\"\n \n

\n\n[![Open in GitHub Codespaces](https://github.com/codespaces/badge.svg)](https://github.com/codespaces/new/nf-core/sammyseq)\n[![GitHub Actions CI Status](https://github.com/nf-core/sammyseq/actions/workflows/nf-test.yml/badge.svg)](https://github.com/nf-core/sammyseq/actions/workflows/nf-test.yml)\n[![GitHub Actions Linting Status](https://github.com/nf-core/sammyseq/actions/workflows/linting.yml/badge.svg)](https://github.com/nf-core/sammyseq/actions/workflows/linting.yml)[![AWS CI](https://img.shields.io/badge/CI%20tests-full%20size-FF9900?labelColor=000000&logo=Amazon%20AWS)](https://nf-co.re/sammyseq/results)[![Cite with Zenodo](http://img.shields.io/badge/DOI-10.5281/zenodo.XXXXXXX-1073c8?labelColor=000000)](https://doi.org/10.5281/zenodo.XXXXXXX)\n[![nf-test](https://img.shields.io/badge/unit_tests-nf--test-337ab7.svg)](https://www.nf-test.com)\n\n[![Nextflow](https://img.shields.io/badge/version-%E2%89%A525.04.0-green?style=flat&logo=nextflow&logoColor=white&color=%230DC09D&link=https%3A%2F%2Fnextflow.io)](https://www.nextflow.io/)\n[![nf-core template version](https://img.shields.io/badge/nf--core_template-3.4.1-green?style=flat&logo=nfcore&logoColor=white&color=%2324B064&link=https%3A%2F%2Fnf-co.re)](https://github.com/nf-core/tools/releases/tag/3.4.1)\n[![run with conda](http://img.shields.io/badge/run%20with-conda-3EB049?labelColor=000000&logo=anaconda)](https://docs.conda.io/en/latest/)\n[![run with docker](https://img.shields.io/badge/run%20with-docker-0db7ed?labelColor=000000&logo=docker)](https://www.docker.com/)\n[![run with singularity](https://img.shields.io/badge/run%20with-singularity-1d355c.svg?labelColor=000000)](https://sylabs.io/docs/)\n[![Launch on Seqera Platform](https://img.shields.io/badge/Launch%20%F0%9F%9A%80-Seqera%20Platform-%234256e7)](https://cloud.seqera.io/launch?pipeline=https://github.com/nf-core/sammyseq)\n\n[![Get help on Slack](http://img.shields.io/badge/slack-nf--core%20%23sammyseq-4A154B?labelColor=000000&logo=slack)](https://nfcore.slack.com/channels/sammyseq)[![Follow on Bluesky](https://img.shields.io/badge/bluesky-%40nf__core-1185fe?labelColor=000000&logo=bluesky)](https://bsky.app/profile/nf-co.re)[![Follow on Mastodon](https://img.shields.io/badge/mastodon-nf__core-6364ff?labelColor=FFFFFF&logo=mastodon)](https://mstdn.science/@nf_core)[![Watch on YouTube](http://img.shields.io/badge/youtube-nf--core-FF0000?labelColor=000000&logo=youtube)](https://www.youtube.com/c/nf-core)\n\n## Introduction\n\n**nf-core/sammyseq** is a bioinformatics pipeline for the analysis of Sequential Analysis of MacroMolecules accessibilitY sequencing (SAMMY-seq) data, a cheap and effective methodology to analyze chromatin state as described in:\n\n> Lucini, F., Petrini, C., Salviato, E., Pal, K., Rosti, V., Gorini, F., Santarelli, P., Quadri, R., Lembo, G., Graziano, G., Di\u00a0Patrizio\u00a0Soldateschi, E., Tagliaferri, I., Pinatel, E., Sebesty\u00e9n, E., Rotta, L., Gentile, F., Vaira, V., Lanzuolo, C., Ferrari, F., 2024. Biochemical properties of chromatin domains define genome compartmentalization. Nucleic Acids Research 52, e54\u2013e54. [doi](https://doi.org/10.1093/nar/gkae454) [pubmed](https://pubmed.ncbi.nlm.nih.gov/38808669/)\n\n> Sebesty\u00e9n, E., Marullo, F., Lucini, F., Petrini, C., Bianchi, A., Valsoni, S., Olivieri, I., Antonelli, L., Gregoretti, F., Oliva, G., Ferrari, F., Lanzuolo, C., 2020. SAMMY-seq reveals early alteration of heterochromatin and deregulation of bivalent genes in Hutchinson-Gilford Progeria Syndrome. Nature Communications 11, 6274. [doi](https://doi.org/10.1038/s41467-020-20048-9) [pubmed](https://pubmed.ncbi.nlm.nih.gov/33293552/)\n\n> [!WARNING]\n> Please note that this pipeline is under active development and has not been released yet.\n\nHere is an outline of the analysis steps:\n\n1. Read QC ([`FastQC`](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/))\n2. Trim reads to remove adapter sequences and low quality ends ([`Trim Galore!`](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore) or [`Trimmomatic`](http://www.usadellab.org/cms/?page=trimmomatic))\n3. Align on a reference genome ([`BWA`](https://bio-bwa.sourceforge.net/) or [`Bowtie 2`](https://bowtie-bio.sourceforge.net/bowtie2))\n4. Mark duplicate reads ([`picard Markduplicates`](http://broadinstitute.github.io/picard))\n5. Filter reads and generate alignment statistics ([`samtools`](http://www.htslib.org/))\n6. Create single track profiles in bigwig format ([`deeptools bamCoverage`](https://deeptools.readthedocs.io/en/latest/))\n7. (Optional) Generate pairwise comparison tracks in bigWig format ([`spp`](https://github.com/hms-dbmi/spp)) or ([`deeptools bigwigCompare`](https://deeptools.readthedocs.io/en/develop/content/tools/bigwigCompare.html)).\n8. (Optional) Identify differentially enriched solubility regions as described in [Wang et al., 2024 ](https://doi.org/10.1038/s41594-025-01622-5).\n9. Generate an analysis report by collecting all generated QC and statistics ([`MultiQC`](http://multiqc.info/))\n\n

\n \n

\n\n## Usage\n\n> [!NOTE]\n> If you are new to Nextflow and nf-core, please refer to [this page](https://nf-co.re/docs/usage/installation) on how to set-up Nextflow. Make sure to [test your setup](https://nf-co.re/docs/usage/introduction#how-to-run-a-pipeline) with `-profile test` before running the workflow on actual data.\n\nFirst, prepare a samplesheet with your input data that looks as follows:\n\n`samplesheet.csv`:\n\n```csv\nsample,fastq_1,fastq_2,experimentalID,fraction,sample_group\nCTRL004_S2,/home/sammy/test_data/CTRL004_S2_chr22only.fq.gz,,CTRL004,S2,CTRL\nCTRL004_S3,/home/sammy/test_data/CTRL004_S3_chr22only.fq.gz,,CTRL004,S3,CTRL\nCTRL004_S4,/home/sammy/test_data/CTRL004_S4_chr22only.fq.gz,,CTRL004,S4,CTRL\n```\n\nEach row represents a fastq file (single-end) or a pair of fastq files (paired end), `experimentalID` represents the biological specimen of interest and `sample` the library produced for each fraction, it usually is a unique combination of `experimentalID` and `fraction`. The `sample_group` field is used to group samples that belong to the same biological condition.\n\nNow, you can run the pipeline using:\n\n```bash\nnextflow run nf-core/sammyseq \\\n -profile \\\n --fasta reference_genome.fa \\\n --input samplesheet.csv \\\n --outdir \n```\n\nor\n\n```bash\nnextflow run nf-core/sammyseq \\\n -profile \\\n --fasta reference_genome.fa \\\n --input samplesheet.csv \\\n --outdir \\\n --comparison S2SvsS3\n```\n\n> [!WARNING]\n> Please provide pipeline parameters via the CLI or Nextflow `-params-file` option. Custom config files including those provided by the `-c` Nextflow option can be used to provide any configuration _**except for parameters**_; see [docs](https://nf-co.re/docs/usage/getting_started/configuration#custom-configuration-files).\n\nFor more details and further functionality, please refer to the [usage documentation](https://nf-co.re/sammyseq/usage) and the [parameter documentation](https://nf-co.re/sammyseq/parameters).\n\n## Pipeline output\n\n\n\nFor more details about the output files and reports, please refer to the\n[output documentation](https://nf-co.re/sammyseq/output).\n\n## Credits\n\nThe SAMMY-seq data analysis procedure was originally developed by the laboratory of Francesco Ferrari (IFOM-ETS, Milan; IGM-CNR, Pavia) in collaboration with the laboratory of Chiara Lanzuolo (INGM, Milan; ITB-CNR, Segrate).\nThe orginal pipeline backbone was mainly the result of work by Cristiano Petrini (IFOM) and Endre Sebesty\u00e9n (IFOM), with significant contributions by Ilario Tagliaferri (IFOM), Giovanni Lembo (IFOM) and Emanuele Di Patrizio Soldateschi (INGM). The project also benefited from the collaboration and input by Eva Maria Pinatel (ITB-CNR). The product of this effort resulted in a first pipeline implemented in bash and adapted to work on Sun Grid Engine (SGE) scheduler.\n\nThe nf-core pipeline (nf-core/sammyseq) is being implemented by [Lucio Di Filippo](https://github.com/lucidif) (ISASI-CNR, Pozzuoli; IBBTEC, Santander), [Ugo Maria Iannacchero](https://github.com/ugoiannacchero) (ITB-CNR), [Nadia Sanseverino](https://github.com/nadiaxaidan) (ISASI-CNR) and [Margherita Mutarelli](https://github.com/daisymut) (ISASI-CNR).\n\n\n\nMany thanks to others who have helped out and contributed along the way too, including (but not limited to): [Phil Ewels](https://github.com/ewels), [Maxime Ulysse Garcia](https://github.com/maxulysse), [Friederike Hanssen](https://github.com/FriederikeHanssen), [Matthias H\u00f6rtenhuber](https://github.com/mashehu), [Marinicla Pascale](https://github.com/Marinicla), [J\u00falia Mir-Pedrol](https://github.com/mirpedrol) and [Marcel Ribeiro-Dantas](https://github.com/mribeirodantas).\n\n## Acknowledgements\n\nThe development of this pipeline was made possible thanks to the projects Progetti@CNR Myo-CoV-2 B93C20046330005, AFM T\u00e9l\u00e9thon EDMD-GenomeSCAN B53C22009260007 and PIR01_00011 I.Bi.S.Co. Infrastruttura per Big data e Scientific COmputing (PON 2014-2020).\n\n## Contributions and Support\n\nIf you would like to contribute to this pipeline, please see the [contributing guidelines](.github/CONTRIBUTING.md).\n\nFor further information or help, don't hesitate to get in touch on the [Slack `#sammyseq` channel](https://nfcore.slack.com/channels/sammyseq) (you can join with [this invite](https://nf-co.re/join/slack)).\n\n## Citations\n\n\n\n\n\nAn extensive list of references for the tools used by the pipeline can be found in the [`CITATIONS.md`](CITATIONS.md) file.\n\nYou can cite the `nf-core` publication as follows:\n\n> **The nf-core framework for community-curated bioinformatics pipelines.**\n>\n> Philip Ewels, Alexander Peltzer, Sven Fillinger, Harshil Patel, Johannes Alneberg, Andreas Wilm, Maxime Ulysse Garcia, Paolo Di Tommaso & Sven Nahnsen.\n>\n> _Nat Biotechnol._ 2020 Feb 13. doi: [10.1038/s41587-020-0439-x](https://dx.doi.org/10.1038/s41587-020-0439-x).\n", + "description": "

\n \n \n \"nf-core/sammyseq\"\n \n

\n\n[![Open in GitHub Codespaces](https://github.com/codespaces/badge.svg)](https://github.com/codespaces/new/nf-core/sammyseq)\n[![GitHub Actions CI Status](https://github.com/nf-core/sammyseq/actions/workflows/nf-test.yml/badge.svg)](https://github.com/nf-core/sammyseq/actions/workflows/nf-test.yml)\n[![GitHub Actions Linting Status](https://github.com/nf-core/sammyseq/actions/workflows/linting.yml/badge.svg)](https://github.com/nf-core/sammyseq/actions/workflows/linting.yml)[![AWS CI](https://img.shields.io/badge/CI%20tests-full%20size-FF9900?labelColor=000000&logo=Amazon%20AWS)](https://nf-co.re/sammyseq/results)[![Cite with Zenodo](http://img.shields.io/badge/DOI-10.5281/zenodo.XXXXXXX-1073c8?labelColor=000000)](https://doi.org/10.5281/zenodo.XXXXXXX)\n[![nf-test](https://img.shields.io/badge/unit_tests-nf--test-337ab7.svg)](https://www.nf-test.com)\n\n[![Nextflow](https://img.shields.io/badge/version-%E2%89%A525.04.0-green?style=flat&logo=nextflow&logoColor=white&color=%230DC09D&link=https%3A%2F%2Fnextflow.io)](https://www.nextflow.io/)\n[![nf-core template version](https://img.shields.io/badge/nf--core_template-3.4.1-green?style=flat&logo=nfcore&logoColor=white&color=%2324B064&link=https%3A%2F%2Fnf-co.re)](https://github.com/nf-core/tools/releases/tag/3.4.1)\n[![run with conda](http://img.shields.io/badge/run%20with-conda-3EB049?labelColor=000000&logo=anaconda)](https://docs.conda.io/en/latest/)\n[![run with docker](https://img.shields.io/badge/run%20with-docker-0db7ed?labelColor=000000&logo=docker)](https://www.docker.com/)\n[![run with singularity](https://img.shields.io/badge/run%20with-singularity-1d355c.svg?labelColor=000000)](https://sylabs.io/docs/)\n[![Launch on Seqera Platform](https://img.shields.io/badge/Launch%20%F0%9F%9A%80-Seqera%20Platform-%234256e7)](https://cloud.seqera.io/launch?pipeline=https://github.com/nf-core/sammyseq)\n\n[![Get help on Slack](http://img.shields.io/badge/slack-nf--core%20%23sammyseq-4A154B?labelColor=000000&logo=slack)](https://nfcore.slack.com/channels/sammyseq)[![Follow on Bluesky](https://img.shields.io/badge/bluesky-%40nf__core-1185fe?labelColor=000000&logo=bluesky)](https://bsky.app/profile/nf-co.re)[![Follow on Mastodon](https://img.shields.io/badge/mastodon-nf__core-6364ff?labelColor=FFFFFF&logo=mastodon)](https://mstdn.science/@nf_core)[![Watch on YouTube](http://img.shields.io/badge/youtube-nf--core-FF0000?labelColor=000000&logo=youtube)](https://www.youtube.com/c/nf-core)\n\n## Introduction\n\n**nf-core/sammyseq** is a bioinformatics pipeline for the analysis of Sequential Analysis of MacroMolecules accessibilitY sequencing (SAMMY-seq) data, a cheap and effective methodology to analyze chromatin state as described in:\n\n> Lucini, F., Petrini, C., Salviato, E., Pal, K., Rosti, V., Gorini, F., Santarelli, P., Quadri, R., Lembo, G., Graziano, G., Di\u00a0Patrizio\u00a0Soldateschi, E., Tagliaferri, I., Pinatel, E., Sebesty\u00e9n, E., Rotta, L., Gentile, F., Vaira, V., Lanzuolo, C., Ferrari, F., 2024. Biochemical properties of chromatin domains define genome compartmentalization. Nucleic Acids Research 52, e54\u2013e54. [doi](https://doi.org/10.1093/nar/gkae454) [pubmed](https://pubmed.ncbi.nlm.nih.gov/38808669/)\n\n> Sebesty\u00e9n, E., Marullo, F., Lucini, F., Petrini, C., Bianchi, A., Valsoni, S., Olivieri, I., Antonelli, L., Gregoretti, F., Oliva, G., Ferrari, F., Lanzuolo, C., 2020. SAMMY-seq reveals early alteration of heterochromatin and deregulation of bivalent genes in Hutchinson-Gilford Progeria Syndrome. Nature Communications 11, 6274. [doi](https://doi.org/10.1038/s41467-020-20048-9) [pubmed](https://pubmed.ncbi.nlm.nih.gov/33293552/)\n\n> [!WARNING]\n> Please note that this pipeline is under active development and has not been released yet.\n\nHere is an outline of the analysis steps:\n\n1. Read QC ([`FastQC`](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/))\n2. Trim reads to remove adapter sequences and low quality ends ([`Trim Galore!`](https://www.bioinformatics.babraham.ac.uk/projects/trim_galore) or [`Trimmomatic`](http://www.usadellab.org/cms/?page=trimmomatic))\n3. Align on a reference genome ([`BWA`](https://bio-bwa.sourceforge.net/) or [`Bowtie 2`](https://bowtie-bio.sourceforge.net/bowtie2))\n4. Mark duplicate reads ([`picard Markduplicates`](http://broadinstitute.github.io/picard))\n5. Filter reads and generate alignment statistics ([`samtools`](http://www.htslib.org/))\n6. Create single track profiles in bigwig format ([`deeptools bamCoverage`](https://deeptools.readthedocs.io/en/latest/))\n7. (Optional) Perform chromatin compartmentalization analysis to identify A/B compartments across chromosomes ([`CALDER2`](https://github.com/CSOgroup/CALDER2)).\n8. (Optional) Generate pairwise comparison tracks in bigWig format ([`spp`](https://github.com/hms-dbmi/spp)) or ([`deeptools bigwigCompare`](https://deeptools.readthedocs.io/en/develop/content/tools/bigwigCompare.html)).\n9. (Optional) Identify differentially enriched solubility regions as described in [Wang et al., 2025 ](https://doi.org/10.1038/s41594-025-01622-5).\n10. Generate an analysis report by collecting all generated QC and statistics ([`MultiQC`](http://multiqc.info/))\n\n

\n \n

\n\n## Usage\n\n> [!NOTE]\n> If you are new to Nextflow and nf-core, please refer to [this page](https://nf-co.re/docs/usage/installation) on how to set-up Nextflow. Make sure to [test your setup](https://nf-co.re/docs/usage/introduction#how-to-run-a-pipeline) with `-profile test` before running the workflow on actual data.\n\nFirst, prepare a samplesheet with your input data that looks as follows:\n\n`samplesheet.csv`:\n\n```csv\nsample,fastq_1,fastq_2,experimentalID,fraction,sample_group\nCTRL004_S2,/home/sammy/test_data/CTRL004_S2_chr22only.fq.gz,,CTRL004,S2,CTRL\nCTRL004_S3,/home/sammy/test_data/CTRL004_S3_chr22only.fq.gz,,CTRL004,S3,CTRL\nCTRL004_S4,/home/sammy/test_data/CTRL004_S4_chr22only.fq.gz,,CTRL004,S4,CTRL\n```\n\nEach row represents a fastq file (single-end) or a pair of fastq files (paired end), `experimentalID` represents the biological specimen of interest and `sample` the library produced for each fraction, it usually is a unique combination of `experimentalID` and `fraction`. The `sample_group` field is used to group samples that belong to the same biological condition.\n\nNow, you can run the pipeline using:\n\n```bash\nnextflow run nf-core/sammyseq \\\n -profile \\\n --fasta reference_genome.fa \\\n --input samplesheet.csv \\\n --outdir \n```\n\nor\n\n```bash\nnextflow run nf-core/sammyseq \\\n -profile \\\n --fasta reference_genome.fa \\\n --input samplesheet.csv \\\n --outdir \\\n --comparison S2SvsS3\n```\n\n> [!WARNING]\n> Please provide pipeline parameters via the CLI or Nextflow `-params-file` option. Custom config files including those provided by the `-c` Nextflow option can be used to provide any configuration _**except for parameters**_; see [docs](https://nf-co.re/docs/usage/getting_started/configuration#custom-configuration-files).\n\nFor more details and further functionality, please refer to the [usage documentation](https://nf-co.re/sammyseq/usage) and the [parameter documentation](https://nf-co.re/sammyseq/parameters).\n\n## Pipeline output\n\n\n\nFor more details about the output files and reports, please refer to the\n[output documentation](https://nf-co.re/sammyseq/output).\n\n## Credits\n\nThe SAMMY-seq data analysis procedure was originally developed by the laboratory of Francesco Ferrari (IFOM-ETS, Milan; IGM-CNR, Pavia) in collaboration with the laboratory of Chiara Lanzuolo (INGM, Milan; ITB-CNR, Segrate).\nThe orginal pipeline backbone was mainly the result of work by Cristiano Petrini (IFOM) and Endre Sebesty\u00e9n (IFOM), with significant contributions by Ilario Tagliaferri (IFOM), Giovanni Lembo (IFOM) and Emanuele Di Patrizio Soldateschi (INGM). The project also benefited from the collaboration and input by Eva Maria Pinatel (ITB-CNR). The product of this effort resulted in a first pipeline implemented in bash and adapted to work on Sun Grid Engine (SGE) scheduler.\n\nThe nf-core pipeline (nf-core/sammyseq) is being implemented by [Lucio Di Filippo](https://github.com/lucidif) (ISASI-CNR, Pozzuoli; IBBTEC, Santander), [Ugo Maria Iannacchero](https://github.com/ugoiannacchero) (ITB-CNR), [Nadia Sanseverino](https://github.com/nadiaxaidan) (ISASI-CNR) and [Margherita Mutarelli](https://github.com/daisymut) (ISASI-CNR).\n\n\n\nMany thanks to others who have helped out and contributed along the way too, including (but not limited to): [Phil Ewels](https://github.com/ewels), [Maxime Ulysse Garcia](https://github.com/maxulysse), [Friederike Hanssen](https://github.com/FriederikeHanssen), [Matthias H\u00f6rtenhuber](https://github.com/mashehu), [Marinicla Pascale](https://github.com/Marinicla), [J\u00falia Mir-Pedrol](https://github.com/mirpedrol) and [Marcel Ribeiro-Dantas](https://github.com/mribeirodantas).\n\n## Acknowledgements\n\nThe development of this pipeline was made possible thanks to the projects Progetti@CNR Myo-CoV-2 B93C20046330005, AFM T\u00e9l\u00e9thon EDMD-GenomeSCAN B53C22009260007 and PIR01_00011 I.Bi.S.Co. Infrastruttura per Big data e Scientific COmputing (PON 2014-2020).\n\n## Contributions and Support\n\nIf you would like to contribute to this pipeline, please see the [contributing guidelines](.github/CONTRIBUTING.md).\n\nFor further information or help, don't hesitate to get in touch on the [Slack `#sammyseq` channel](https://nfcore.slack.com/channels/sammyseq) (you can join with [this invite](https://nf-co.re/join/slack)).\n\n## Citations\n\n\n\n\n\nAn extensive list of references for the tools used by the pipeline can be found in the [`CITATIONS.md`](CITATIONS.md) file.\n\nYou can cite the `nf-core` publication as follows:\n\n> **The nf-core framework for community-curated bioinformatics pipelines.**\n>\n> Philip Ewels, Alexander Peltzer, Sven Fillinger, Harshil Patel, Johannes Alneberg, Andreas Wilm, Maxime Ulysse Garcia, Paolo Di Tommaso & Sven Nahnsen.\n>\n> _Nat Biotechnol._ 2020 Feb 13. doi: [10.1038/s41587-020-0439-x](https://dx.doi.org/10.1038/s41587-020-0439-x).\n", "hasPart": [ { "@id": "main.nf" diff --git a/subworkflows/local/compartmentalization_analysis/main.nf b/subworkflows/local/compartmentalization_analysis/main.nf new file mode 100644 index 00000000..2a062001 --- /dev/null +++ b/subworkflows/local/compartmentalization_analysis/main.nf @@ -0,0 +1,170 @@ +// +// Run genome compartmentalization analysis +// + +include { CHR_SPLIT } from '../../../modules/local/chr_split/main' +include { CHR_COMPARTMENTS_CALLING } from '../../../modules/local/chr_compartments_calling/main' +include { CHR_COMBINE_COMPARTMENTS } from '../../../modules/local/chr_combine_compartments/main' +include { BUILD_CONSENSUS } from '../../../modules/local/build_consensus/main' + +workflow COMPARTMENTALIZATION_ANALYSIS { + + take: + bigwig_tracks // channel: [ meta, bigwig ] + binned_genome // channel: path BED file + chrom_sizes // path: chromosome sizes file + outdir // val: output directory + binsize // val: bin size + gtf // path: GTF file + + main: + + ch_versions = Channel.empty() + + // + // Extract track information from BigWig files + // + ch_compartmentTracks = bigwig_tracks + .map { meta, bigwig -> + [meta.experimentalID, meta.fraction, meta.sample_group, bigwig] + } + + // + // Create CSV file with all tracks + // + ch_compartments_csv = ch_compartmentTracks + .map { experimentalID, fraction, sample_group, bigwig -> + "${experimentalID},${fraction},${sample_group},${bigwig}\n" + } + .collectFile( + name: 'compartments_tracks.csv', + seed: "Patient_name,Fraction,Status,File\n", + storeDir: "${outdir}/csv", + sort: true + ) + + // + // Get unique samples (experimentalID) + // + ch_uniqueSamples = ch_compartmentTracks + .map { experimentalID, fraction, sample_group, bigwig -> experimentalID } + .unique() + + // + // Create meta for BED file + // + ch_bed_with_meta = binned_genome + .map { bed -> + [ [id: bed.baseName], bed ] + } + + // + // Split binned genome by chromosome + // + CHR_SPLIT( + ch_bed_with_meta, + chrom_sizes + ) + + // + // Transpose to get one item per chromosome file + // + ch_chromBeds = CHR_SPLIT.out.beds + .transpose() + .map { meta, chr_bed -> + def chr_name = chr_bed.baseName + def new_meta = meta + [chromosome: chr_name] + [ new_meta, chr_bed ] + } + + // + // Combine each chromosome with each unique sample and CSV + // + ch_chromSampleTuples = ch_chromBeds + .combine(ch_uniqueSamples) + .combine(ch_compartments_csv) + .map { meta, chr_bed, patient, csv -> + [ meta, chr_bed, patient, csv ] + } + + // + // Call compartments for each chromosome/sample combination + // + CHR_COMPARTMENTS_CALLING( + ch_chromSampleTuples, + binsize, + gtf + ) + ch_versions = ch_versions.mix(CHR_COMPARTMENTS_CALLING.out.versions.first()) + + // + // Group BED files by sample + // + ch_beds_by_sample = CHR_COMPARTMENTS_CALLING.out.bed_files + .groupTuple(by: 0) + + // + // Group BedGraph files by sample + // + + ch_bedgraphs_by_sample = CHR_COMPARTMENTS_CALLING.out.bedgraph_files + .groupTuple(by: 0) + + // + // Combine BED and BedGraph channels for same sample + // + ch_combine_input = ch_beds_by_sample + .join(ch_bedgraphs_by_sample) + + // + // Merge all chromosomes per sample + // + CHR_COMBINE_COMPARTMENTS( + ch_combine_input + ) + ch_versions = ch_versions.mix(CHR_COMBINE_COMPARTMENTS.out.versions.first()) + + // + // Extract sample_group from original tracks + // + ch_sample_groups = ch_compartmentTracks + .map { experimentalID, fraction, sample_group, bigwig -> + [experimentalID, sample_group] + } + .unique() + + // + // Add sample_group to combined beds + // + ch_beds_with_group = CHR_COMBINE_COMPARTMENTS.out.combined_beds + .combine(ch_sample_groups) + .filter { patient_bed, bed, patient_group, group -> + patient_bed == patient_group + } + .map { patient_bed, bed, patient_group, group -> + [group, patient_bed, bed] + } + + // + // Group combined BEDs by sample_group + // + ch_consensus_input = ch_beds_with_group + .map { group, patient, bed -> [group, bed] } + .groupTuple() + + // + // Generate majority and strict consensus + // + BUILD_CONSENSUS( + ch_consensus_input + ) + + emit: + bed_files = CHR_COMPARTMENTS_CALLING.out.bed_files // channel: [ patient, bed ] + bedgraph_files = CHR_COMPARTMENTS_CALLING.out.bedgraph_files // channel: [ patient, bedgraph ] + combined_beds = CHR_COMBINE_COMPARTMENTS.out.combined_beds // channel: [ patient, bed ] + combined_bedgraphs = CHR_COMBINE_COMPARTMENTS.out.combined_bedgraphs // channel: [ patient, bedgraph ] + consensus_majority = BUILD_CONSENSUS.out.consensus_majority // channel: [ path(bed) ] + consensus_strict = BUILD_CONSENSUS.out.consensus_strict // channel: [ path(bed) ] + versions = ch_versions // channel: [ versions.yml ] +} diff --git a/subworkflows/local/compartmentalization_analysis/meta.yml b/subworkflows/local/compartmentalization_analysis/meta.yml new file mode 100644 index 00000000..e8790e65 --- /dev/null +++ b/subworkflows/local/compartmentalization_analysis/meta.yml @@ -0,0 +1,78 @@ +# yaml-language-server: $schema=https://raw.githubusercontent.com/nf-core/modules/master/subworkflows/yaml-schema.json +name: "compartmentalization_analysis" +description: Perform genome compartmentalization analysis by splitting genome into chromosomes and calling compartments per sample +keywords: + - compartmentalization + - chromatin + - epigenetics + - chromosome + - bigwig +components: + - chromosome_split + - compartments_calling + - combine_compartments +input: + - bigwig_tracks: + type: file + description: | + Channel containing BigWig coverage tracks with metadata + Structure: [ meta, bigwig ] + pattern: "*.{bigWig,bw}" + - binned_genome: + type: file + description: | + BED file containing binned genome intervals (all chromosomes) + pattern: "*.bed" + - chrom_sizes: + type: file + description: | + File containing chromosome names and sizes (tab-separated) + pattern: "*.{txt,sizes,chrom.sizes}" + - outdir: + type: string + description: Output directory path for results + - binsize: + type: integer + description: Bin size used for genome binning (in base pairs) + - gtf: + type: file + description: GTF annotation file for gene features + pattern: "*.gtf" +output: + - bed_files: + type: file + description: | + BED files containing compartment assignments per patient and chromosome + Structure: [ patient, bed ] + pattern: "*_compartments.bed" + - bedgraph_files: + type: file + description: | + BedGraph files containing eigenvector values per patient and chromosome + Structure: [ patient, bedgraph ] + pattern: "*_comp_eigenvector.bedgraph" + - combined_beds: + type: file + description: | + Combined BED files with all chromosomes per patient + Structure: [ patient, bed ] + pattern: "*_combined_compartments.bed" + - combined_bedgraphs: + type: file + description: | + Combined BedGraph files with all chromosomes per patient + Structure: [ patient, bedgraph ] + pattern: "*_combined_compartments_eigen.bedgraph" + - versions: + type: file + description: File containing software versions + pattern: "versions.yml" +notes: + - | + The CALDER algorithm may fail with "Non-positive eigenvalues" error for samples + with insufficient coverage or when the correlation matrix is singular. + Consider using errorStrategy 'ignore' or increasing binsize if this occurs. +authors: + - "@ugoiannacchero" +maintainers: + - "@daisymut" diff --git a/tests/default.nf.test.snap b/tests/default.nf.test.snap index ff904971..2a1d66da 100644 --- a/tests/default.nf.test.snap +++ b/tests/default.nf.test.snap @@ -327,9 +327,9 @@ ] ], "meta": { - "nf-test": "0.9.3", + "nf-test": "0.9.2", "nextflow": "25.10.0" }, - "timestamp": "2025-11-06T12:23:44.709651414" + "timestamp": "2026-07-25T13:02:27.28887" } } \ No newline at end of file diff --git a/workflows/sammyseq.nf b/workflows/sammyseq.nf index d8a21e40..ccc06fc4 100644 --- a/workflows/sammyseq.nf +++ b/workflows/sammyseq.nf @@ -39,6 +39,7 @@ include { BIGWIG_PLOT_DEEPTOOLS } from '../subworkflows/local/bigw include { DEEPTOOLS_QC } from '../subworkflows/local/deeptools_qc' include { GENERATE_COMPARISONS } from '../subworkflows/local/generate_comparisons' include { DIFFERENTIAL_SOLUBILITY_ANALYSIS } from '../subworkflows/local/differential_solubility_analysis' +include { COMPARTMENTALIZATION_ANALYSIS } from '../subworkflows/local/compartmentalization_analysis' /* @@ -272,24 +273,37 @@ if (params.stopAt == 'ALIGNMENT') { return } + +// Initialise so downstream (MultiQC) is safe when QC is skipped +ch_dt_corrmatrix = Channel.empty() +ch_dt_pcadata = Channel.empty() +ch_dt_fpmatrix_global = Channel.empty() +ch_dt_fpmetrics_global = Channel.empty() +ch_dt_fpmatrix_region = Channel.empty() +ch_dt_fpmetrics_region = Channel.empty() + +if (!params.skip_deeptools_qc) { DEEPTOOLS_QC ( - FILTER_BAM_SAMTOOLS.out.bam, - FILTER_BAM_SAMTOOLS.out.bai, - DEEPTOOLS_BAMCOVERAGE.out.bigwig, - params.corr_method, - params.blacklist ? PREPARE_GENOME.out.blacklist : Channel.value(tuple([ id:'no_blacklist' ], [])) + FILTER_BAM_SAMTOOLS.out.bam, + FILTER_BAM_SAMTOOLS.out.bai, + DEEPTOOLS_BAMCOVERAGE.out.bigwig, + params.corr_method, + params.blacklist ? PREPARE_GENOME.out.blacklist : Channel.value(tuple([ id:'no_blacklist' ], [])) ) - ch_dt_corrmatrix = DEEPTOOLS_QC.out.correlation_matrix - ch_dt_pcadata = DEEPTOOLS_QC.out.pca_data + ch_dt_corrmatrix = DEEPTOOLS_QC.out.correlation_matrix + ch_dt_pcadata = DEEPTOOLS_QC.out.pca_data + if (params.plotfingerprint) { - ch_dt_fpmatrix_global = DEEPTOOLS_QC.out.fingerprint_matrix_global + ch_dt_fpmatrix_global = DEEPTOOLS_QC.out.fingerprint_matrix_global ch_dt_fpmetrics_global = DEEPTOOLS_QC.out.fingerprint_metrics_global if (params.region) { - ch_dt_fpmatrix_region = DEEPTOOLS_QC.out.fingerprint_matrix_region + ch_dt_fpmatrix_region = DEEPTOOLS_QC.out.fingerprint_matrix_region ch_dt_fpmetrics_region = DEEPTOOLS_QC.out.fingerprint_metrics_region } } ch_versions = ch_versions.mix(DEEPTOOLS_QC.out.versions) +} + if (params.tss_bed) { ch_bw_by_sample = DEEPTOOLS_BAMCOVERAGE.out.bigwig @@ -311,14 +325,14 @@ if (params.stopAt == 'ALIGNMENT') { } // - // GENOME BINNING: Run only if comparison_file or comparison is provided + // GENOME BINNING: Run only if comparison_file or comparison or compartmentalization_anaylsis is provided // if (params.comparison_file && params.comparison) { error "Cannot specify both --comparison_file and --comparison parameters. Please use only one method." } - if (params.comparison_file || params.comparison) { + if (params.comparison_file || params.comparison || params.compartmentalization_analysis) { GENOME_BINNING( PREPARE_GENOME.out.filtered_bed, params.keep_regions_bed, @@ -377,6 +391,27 @@ if (params.stopAt == 'ALIGNMENT') { } } + // + // COMPARTMENTALIZATION ANALYSIS + // + + if (params.compartmentalization_analysis) { + + if (!params.gtf) { + exit 1, "ERROR: The --gtf parameter must be provided when --compartmentalization_analysis is enabled." + } + + COMPARTMENTALIZATION_ANALYSIS( + DEEPTOOLS_BAMCOVERAGE.out.bigwig, + ch_genome_bins, + PREPARE_GENOME.out.chrom_sizes, + params.outdir, + params.binsize, + PREPARE_GENOME.out.gtf + ) + + } + // // Collate and save software versions // @@ -429,12 +464,12 @@ if (params.stopAt == 'ALIGNMENT') { ch_multiqc_files = ch_multiqc_files.mix(FILTER_BAM_SAMTOOLS.out.flagstat.collect{it[1]}.ifEmpty([])) ch_multiqc_files = ch_multiqc_files.mix(FILTER_BAM_SAMTOOLS.out.idxstats.collect{it[1]}.ifEmpty([])) - ch_multiqc_files = ch_multiqc_files.mix(DEEPTOOLS_QC.out.correlation_matrix.collect{it[1]}.ifEmpty([])) - ch_multiqc_files = ch_multiqc_files.mix(DEEPTOOLS_QC.out.pca_data.collect{it[1]}.ifEmpty([])) + ch_multiqc_files = ch_multiqc_files.mix(ch_dt_corrmatrix.collect{it[1]}.ifEmpty([])) + ch_multiqc_files = ch_multiqc_files.mix(ch_dt_pcadata.collect{it[1]}.ifEmpty([])) if (params.plotfingerprint) { - ch_multiqc_files = ch_multiqc_files.mix(DEEPTOOLS_QC.out.fingerprint_matrix_global.collect{it[1]}.ifEmpty([])) - ch_multiqc_files = ch_multiqc_files.mix(DEEPTOOLS_QC.out.fingerprint_metrics_global.collect{it[1]}.ifEmpty([])) + ch_multiqc_files = ch_multiqc_files.mix(ch_dt_fpmatrix_global.collect{it[1]}.ifEmpty([])) + ch_multiqc_files = ch_multiqc_files.mix(ch_dt_fpmetrics_global.collect{it[1]}.ifEmpty([])) } if (params.tss_bed) {