Hello, I have a study of 70 samples. There are three study participants that have duplicate samples (two technical replicates for each of those participants). Is it possible to have these samples merged during the differentialabundance pipeline? It doesnt seem to be clear from the documentation how to merge reads/counts from these samples? Thank you, Nathan

I have run the nfcore rnaseq pipeline for a group here and subsequently ran the differential abundance pipeline on it. We ran the rnaseq pipeline using star + salmon on mouse data. Going through some of the DESeq2 results I have some odd results I would love your input on. There are genes that pass the DESeq2 results as significant but show some very large foldchange such as (Looking at WDBRN vs CDBRN):

gene_id	baseMean	log2FoldChange	lfcSE	pvalue	padj
ENSMUSG00000026180.10	0.9640223	-14.02183	4.95085698	2.061986e-05	0.00289795

The basemean is fairly low which made me think some of these might have low counts. Filtering on basemean alone does not seem to be a valid reason (https://support.bioconductor.org/p/9143776/#9143785) I took a look at the all.normalised_counts.tsv:

gene_id	CD1BRN	CD2BRN	CD3BRN	CD4BRN	PRORA1BRN	PRORA2BRN	PRORA3BRN	PRORA4BRN	WD1BRN	WD2BRN	WD3BRN	WD4BRN
ENSMUSG00000026180.10	0	3.219497278	0	0	1.93365157	0	6.415119346	0	0	0	0	0

all.vst.tsv:

gene_id	CD1BRN	CD2BRN	CD3BRN	CD4BRN	PRORA1BRN	PRORA2BRN	PRORA3BRN	PRORA4BRN	WD1BRN	WD2BRN	WD3BRN	WD4BRN
ENSMUSG00000026180.10	6.37715807	6.66063142	6.37715807	6.37715807	6.59698741	6.37715807	6.77667224	6.37715807	6.37715807	6.37715807	6.37715807	6.37715807

The raw counts from salmon:

salmon.merged.gene_counts.tsv													
gene_id	gene_name	CD1BRN	CD2BRN	CD3BRN	CD4BRN	PRORA1BRN	PRORA2BRN	PRORA3BRN	PRORA4BRN	WD1BRN	WD2BRN	WD3BRN	WD4BRN
ENSMUSG00000026180.10	Cxcr2	0	3	0	0	2	0	6	0	0	0	0	0

and finally the gene length: salmon.merged.gene_lengths.tsv

gene_id	gene_name	CD1BRN	CD2BRN	CD3BRN	CD4BRN	PRORA1BRN	PRORA2BRN	PRORA3BRN	PRORA4BRN	WD1BRN	WD2BRN	WD3BRN	WD4BRN
ENSMUSG00000026180.10	Cxcr2	2769.29671	2760.794	2769.29671	2769.29671	2773.335	2769.29671	2769.456	2769.29671	2769.29671	2769.29671	2769.29671	2769.29671

This is a fairly low expression gene which is pretty large in size. I don't understand why it comes up as significant in DESeq2 but if you plot the normalised vst counts there is no difference between samples. How should we remove such genes from our analysis and how are you dealing with it in your datasets? Thanks for your help in understanding this better!

Hello I just ran the differentialabundance_1.5.0 and I noticed in the study.html, description for the table says the following “_The

DESeq2 R

package was used for differential analysis. p-values were adjusted with the BH method to reduce the number of false positives. Genes were considered differential if, for the respective contrast, the adjusted p-value was equal to or lower than 0.1 and the absolute log2 fold change was equal to or higher than 0.”_ But looking at the plots below and the data in *deseq2.results_filtered.tsv the adjusted p-value was equal to or lower than 0.05 and the absolute log2 fold change was equal to or higher than 1. can you help me understand if the description mentioned above refers to something else ? Thank you

Hello, I am searching in how to do differential abundance analysis but in genomic information instead of transcriptomic. I have two site with three different conditions in which I would like to perform DAA to see if any of the conditions enrich a specific set of genes at the genomic level. Could this pipeline be use for that? Becuase in the documentation it says that it supports RNA-seq data but is it exclusively RNA-seq data?

Hi, is it possible to run differential abundance on samples from different nf-core/rnaseq runs or not because of the dependency on some shared count matrix?

Hi all, Is there a way to filter for minimum samples in each DESeq2 pairwise comparison? I have been running the differential abundance pipeline and within each DESeq2 pairwise comparison, I am struggling with the low number of samples (with abundance), per comparison. I have been running the pipeline using 36 samples (6 treatments) but only interested in 4 of the pairwise comparisons. The FAQs in the DESeq2 vignette suggest '_we recommend users to run samples from all groups together'_, which makes me think other people have encountered something similar? Is there a way to filter for a proportion of samples per pairwise comparison? e.g. 3 out of 12 (25%)? Many thanks for your help! If I adjust default prefilter params,

// Filtering options
 filtering_min_samples           = 5.0
 filtering_min_abundance         = 1.0

it reduces the volume,

Input was a matrix of 69389 genes for 36 samples, reduced to 41690 genes after filtering for low abundance.

but also provides some comparison DEG results

gene_id	baseMean	log2FoldChange	lfcSE	pvalue	padj
ENSSSAG00000006734	0.3893752	-11.70838	10.69647484	0.00002290416	0.02075814

with abundance in 1 sample out of 12

D1L	D1L	D1L	D1L	D1L	D1L	D1PAIR	D1PAIR	D1PAIR	D1PAIR	D1PAIR	D1PAIR
gene_id	H9_1	H9_2	H17_1	H17_2	H1_1	H1_2	H2_1	H2_2	H7_1	H7_2	H15_1	H15_2
ENSSSAG00000006734	0	0	0	0	3.548172273	0	0	0	0	0	0	0

across all 36 samples

D1H	D1H	D2PAIR	D2PAIR	D2H	D2H	D2L	D2L	D2H	D2H	D1PAIR	D1PAIR	D1H	D1H	D1L	D1L	D2L	D2L	D1L	D1L	D1PAIR	D1PAIR	D1H	D1H	D2H	D2H	D2L	D2L	D2PAIR	D2PAIR	D1PAIR	D1PAIR	D2PAIR	D2PAIR	D1L	D1L
gene_id	H10_1	H10_2	H11_1	H11_2	H12_1	H12_2	H13_1	H13_2	H14_1	H14_2	H15_1	H15_2	H16_1	H16_2	H17_1	H17_2	H18_1	H18_2	H1_1	H1_2	H2_1	H2_2	H3_1	H3_2	H4_1	H4_2	H5_1	H5_2	H6_1	H6_2	H7_1	H7_2	H8_1	H8_2	H9_1	H9_2
ENSSSAG00000006734	0	0	0	0	0	2.014868833	0	1.16021696	0	0	0	0	0	0	0	0	1.099274509	0	3.548172273	0	0	0	0	0	0	0	0	0	0	0	0	0	0	5.050406692	0	0

Dear <!channel>, We hope you’re all having a great time at the Nextflow Summit! nextflow @suzannejin and I are excited about the upcoming nf core hackathon next week. We’d really appreciate your input on your experiences with the #C045UNCS5R9 pipeline and your thoughts on future developments. Please take a moment to complete our 👉 short survey — your feedback means a lot to us!

What I used to run the Rshiny app locally

Rscript -e 'library(methods); shiny::runApp("PATH_TO_APP_FOLDER", launch.browser = TRUE)'

I am trying out this pipeline on some data I manually analyzed and am getting some pretty different results. In my manual analysis, I was using DESeq2 with interaction terms. I am hoping my results are different this time because I don't fully understand how I should utilize the blocking function. Is there a way I could chat with someone in a video call to try to reconcile the results?

Hi, all. I've still been working at reconciling the differential results from this pipeline with manual DESeq2 calls. I am using the same count matrix as input and used the same model reported in the MultiQC file for the design in my DESeq2 experiment. I am doing the same contrast on the same data with the same design. However, the results I see from manually running the analysis are completely different than what the pipeline outputs.

Long shot, but does anyone know the study, experiment, or paper where the full test data comes from? It is a knockout using a cre system but I don't know what exactly was knocked out.

Hi all, Perhaps a straightforward question, but I couldn't find an obvious answer. Is the pipeline supporting LRT as well as Wald for DESeq2? I see it's available as an option in

--deseq2_test

, but wasn't sure if this was functional? 🙂

Hi All, I noticed from my latest run of the pipeline (

dev

branch), that the p-values reported tend to fall in discrete clusters instead of being more continuously spread. For example, sorting from smallest to largest, the p-values I'm seeing are a series of the following: 0 1.00E-04 2.00E-04 3.00E-04 4.00E-04 5.00E-04 6.00E-04 7.00E-04 etc. Dozens or hundreds of genes have these same p-values, and there are no values in between these groups of numbers. In prior runs of the pipeline, p-values were more continuously spread (e.g., I might have values between 0 and 1.00E-04 or between 1.00E-04 and 2.00E-04). Is this difference due to some change in the way that the pipeline (or DESeq2) is rounding the p-values, or is this a peculiarity of my specific data set?

Hi all, I'm trying to wrap my head around how I could do a nested experimental design with Differential abundance pipeline, based on this documentation from the DESeq2 vignette I'm trying to follow a design with: • 2 groups (X and Y) • 3 patients per group (1,2,3 and 4,5,6) • 2 conditions (A and B) per group The goal is to to look at if the effects of A vs B differ between groups X and Y, while controlling for individual. (I've attached a diagram to illustrate this) I can encode the table like so:

##    grp ind cnd ind.n
## 1    X   1   A     1
## 2    X   1   B     1
## 3    X   2   A     2
## 4    X   2   B     2
## 5    X   3   A     3
## 6    X   3   B     3
## 7    Y   4   A     1
## 8    Y   4   B     1
## 9    Y   5   A     2
## 10   Y   5   B     2
## 11   Y   6   A     3
## 12   Y   6   B     3

From the documentation, this suggests the model can be:

counts ~ grp + grp:ind.n + grp:cnd

However, I'm failing to see how we can translate interaction terms into the samplesheet and contrasts? The closest that I can think of/manage is that we can instead have a column

grp_ind

that acts to combine the `grp`and

ind

columns (as they are linearly related)

##    grp ind cnd individual_and_group
## 1    X   1   A     X_1
## 2    X   1   B     X_1
## 3    X   2   A     X_2
## 4    X   2   B     X_2
## 5    X   3   A     X_3
## 6    X   3   B     X_3
## 7    Y   4   A     Y_4
## 8    Y   4   B     Y_4
## 9    Y   5   A     Y_5
## 10   Y   5   B     Y_5
## 11   Y   6   A     Y_6
## 12   Y   6   B     Y_6

As far as I can tell, with the pipeline, we can only model A vs B within each group (X and Y) accounting for patient, and then have to do I have to the group-specific condition effects outside the pipeline? My contrasts file would therefore, from what I can tell, look like:

id,variable,reference,target,blocking,exclude_samples_col,exclude_samples_values
X_cnd_B_vs_A,cnd,A,B,ind_grp,grp,Y
Y_cnd_B_vs_A,cnd,A,B,ind_grp,grp,X

where the first row is looking at A vs B only in group X, accounting for patient ID. and the second row is the same except for group Y? Hopefully what I've laid out above makes sense and I hope I haven't made an incredibly obvious blunder, thanks in advance 🙂

I plan to run the rnaseq pipeline totalRNA-sequenced results, and then run DA on that output. However, I am only interested in the protein coding results. Can I just filter out all of the gene IDs from the tsv count files that do not have the protein_coding annotation? Will there be any complications from this approach?

Hej! I'm trying to set up the minRep value. However it only accepts numbers and not "Inf"

deseq2_min_replicates_for_replace
 'minReplicatesForReplace' parameter passed to DESeq() 
Type expected for this parameter's value (null) is different from the one defined in the pipeline schema. Must be a number.
the minimum number of replicates required in order to use replaceOutliers on a sample. If there are samples with so many replicates, the model will be refit after these replacing outliers, flagged by Cook's distance. Set to Inf in order to never replace outliers. It set to Inf for fitType="glmGamPoi".

So, I'm inventing a number here to avoid it: --deseq2_min_replicates_for_replace 1000000 Is it ok to do? Why I can't pass minRep Inf? And null does not mean the same.

Hi all, I was reading through the dev documentation (soon to be 2.0) and came across this line:

This gives full control over the contrast definition but requires understanding of the model matrix. Some downstream applications (e.g. GSEA_GSEA, SHINYNGS_APP) do not support formula-based contrasts as they require a meta.variable.

I'm just wanting to confirm this will also be the case for 2.0? And that we may have to custom post-hoc GSEA for formula based designs (or design our own modules? 🙂) Thanks in advance!

Don't know if this of interest, but there's a problem in the way the digits are rounded (in R scripts, possibly in general). If you look at this line https://github.com/nf-core/modules/blob/master/modules/nf-core/deseq2/differential/templates/deseq2_differential.R#L485 the digits are rounded using round, which means that if you're using 8 digits for example, everything below 1E-8 (or 1E-9) will become 0. Since the default is 4 digits, it means everything below 1E-4. I suggest replacing it with something like

apply(data.frame(comp.results), 2, function(x) sprintf(paste0("%." opt\$round_digits, "g"), x))

This will print everything below significant digits in E format, but still keep a fixed amount of significant digits. I can open an issue, but this will need to be fixed on multiple modules. Would you like me to do that? Maybe this isn't a problem, and should stay the way it is? Should I open an issue on differentialabundane, and then have sub-issues for each module with a PR?

Hi, can the output from nf-core: rnaseq (https://github.com/nf-core/rnaseq?tab=readme-ov-file) salmon.merged.gene_counts.tsv be used directly as sample input to this nf-core:differentialabundance pipeline? Thanks!

Hey differentialabundance team, the current release 1.5.0

test_full

of differentialabundance seems to be broken due to a non existing file from the test-datasets repo: https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/rnaseq_featurecounts_merged_gene_counts.txt

ERROR ~ No such file or directory: <https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/rnaseq_featurecounts_merged_gene_counts.txt>

On dev, the gene_counts parameter points to the right file: https://raw.githubusercontent.com/nf-core/test-datasets/differentialabundance/testdata/rnaseq_featurecounts_merged_gene_counts.tsv

There is a small https://github.com/nf-core/differentialabundance/pull/737 PR open for making that 2.0 🙂 Mind reviewing? 😛

Hi differentialabundance team, I'm hitting a bug in differentialabundance 2.0.0 when I pass multiple gene set files via

gene_sets_files

(comma-separated): the GSEA_GSEA task fails with exit code 1 during post-processing

Environment: - nf-core/differentialabundance v2.0.0 - Nextflow 26.04.4 - Google Batch executor - 13 contrasts, 2 GMT files (Reactome + Hallmark)

Hi differentialabundance team, great work on version 2.0.0. I absolutely love the new report. Just a small question regarding the heatmap. Is it possible to increase the size of the plot or leave out the legend? I have a lot of metadata and they basically overlay the heatmap completely.

Hi, first of all, thank you for this great pipeline. In our centre we are usually using DIA-NN (or specifically bigBio:quantmsdiann) output for proteomics data, which requires some normalisation. I previously implemented a new 'study_type' for general 'mass_spec' with a simple count matrix in a fork and added different limma normalisation methods for it (including variance stabilised normalisation). This was version 1.5.0. I started porting this to the version 2, but I realised that you completely removed any non nf-core modules. Before I would get it ready for a PR, in case this is of general interest, would it be okay to add the limma normalisation processes simply in

modules/local

?

Hey folks, I ran into something with differentialabundance and 3′ tag/DGE RNA-seq (QuantSeq, Lexogen 3′ FWD, etc.) and wanted to bring it here first, since I'd like to help fix it rather than just file it and walk away. The DESeq2/limma paths (2.0.0 docs) expect a length-scaled or length-offset gene matrix from nf-core/rnaseq. That's right for full-length libraries, but for 3′ tag data the count per gene is length-independent, so a length offset would bias fold changes by transcript length. Standard handling is CPM via edgeR, no length correction, DESeqDataSetFromMatrix on raw counts. I didn't see this distinction anywhere in the usage/output docs — might just be something I missed. Wherever this lands, I'd like to pitch in: - If it's just a docs gap, I can write the pull request (a warning + guidance for 3′ users on which matrix to supply). - If it's worth a real parameter (e.g. --library_type 3prime|full_length switching the normalization path) — I can put together a first pass and iterate with whoever knows the DESeq2/limma module best. - Either way I have a 3′ DGE dataset I can test against, so I'm not just guessing at the fix. Let me know which direction makes sense, or if I'm missing existing support for this and it's a non-issue.

Hi all, FYI, I was testing the workflow with ATACSEQ data, I got the following error in

NFCORE_DIFFERENTIALABUNDANCE:DIFFERENTIALABUNDANCE:QUARTONOTEBOOK

(...)
   12.               │ └─base::withCallingHandlers(...)
   13.               ├─purrr:::call_with_cleanup(...)
   14.               └─.f(.x[[i]], ...)
   15.                 └─global generate_boxplots(combined, label_col)
   16.                   └─purrr::walk(...)
   17.                     └─purrr::map(.x, .f, ..., .progress = .progress)
   18.                       └─purrr:::map_("list", .x, .f, ..., .progress = .progress)
   19.                         ├─purrr:::with_indexed_errors(...)
   20.                         │ └─base::withCallingHandlers(...)
   21.                         ├─purrr:::call_with_cleanup(...)
   22.                         └─.f(.x[[i]], ...)
   23.                           ├─plotly::layout(...)
   24.                           ├─plotly::ggplotly(p_box, tooltip = c("x", "y", "color", "text"))
   25.                           └─plotly:::ggplotly.ggplot(...)
   26.                             └─plotly::gg2list(...)
   27.                               └─plotly (local) ggplotly_build(p)
   28.                                 └─base::lapply(data, scales_map_df, scales = npscales)
   29.                                   └─plotly (local) FUN(X[[i]], ...)
   30.                                     └─scales$map_df(df)
   31.                                       └─ggplot2 (local) map_df(..., self = self)
   32.                                         ├─base::unlist(...)
   33.                                         └─base::lapply(self$scales, function(scale) scale$map_df(df = df))
   34.                                           └─ggplot2 (local) FUN(X[[i]], ...)
   35.                                             └─scale$map_df(df = df)
   36.                                               └─ggplot2 (local) map_df(..., self = self)
   37.                                                 └─base::lapply(aesthetics, function(j) self$map(df[[j]]))
   38.                                                   └─ggplot2 (local) FUN(X[[i]], ...)
   39.                                                     └─self$map(df[[j]])
   40.                                                       └─ggplot2 (local) map(..., self = self)
   41.                                                         └─self$palette(n)
   42.                                                           └─ggplot2 (local) palette(...)
   43.                                                             └─ggplot2 (local) palette(...)
   44.                                                               └─ggplot2 (local) palette(...)
   45.                                                                 └─cli::cli_abort("Insufficient values in manual scale. {n} needed but only {length(values)} provided.")
   46.                                                                   └─rlang::abort(...)

nf-core/differentialabundance 2.0.0. my input was:

/nextflow run  nf-core/differentialabundance  \
                -c nfcore.cfg \
        -profile apptainer \
        -resume \
        --outdir "results" \
        --input "samplesheet.csv \
        --contrasts contrasts.csv \
        --matrix  featureCount.tsv \
        --features_metadata_cols "gene_id,gene_name" \
        --features_id_col gene_id \
        --features_name_col "gene_name"

Pierre

Hi differentialabundance team, I noticed that the report html shows blank interactive plots when using --functional_method gsea (Version 2.0). While the test profile works fine, we observed it for the full_test and with our own data. The cause might be GSEA results having a trailing tab on every line. As workaround it worked to adapt the quarto template to drop empty-named columns, but maybe there are better ways to prevent/fix it?

Hi all, I tested the differentialabundance pipeline with standard settings on one of our datasets. The results look very good! However, while checking the results I found a rare number of genes to be significantly different, but seem only found expressed in ONE individual sample (of 5x2 samples covering two conditions). I think there may be outliers among the samples influencing this outcome, but generally: is this expected behaviour? Do you have any recommendations on how to best address this, e.g., by setting filtering settings ("--filtering_min_abundance")? I´m grateful for any advise. Best wishes, Christian.

depictio Depictio dashboard for differentialabundance: your feedback wanted Hi all! 👋 I am building Depictio, an open-source platform from EMBL to explore pipeline results without writing a Shiny or Streamlit app each time. A pipeline gets a dashboard template: a short definition that turns the outputs of any run into an interactive dashboard. I am now sharing it more widely, so I would like the view of differentialabundance users and developers. ✨ What it does • MultiQC reports (aggregated across runs) next to result tables, images, phylogenetic trees, maps... with linked filters • Reusable: point the template at any new run of the same pipeline and you get the same dashboard • Fully editable: modify any dashboard from the GUI or in YAML, so the template is a starting point, not a constraint • Versioned, with comments and annotations on any component • Runs as a Python package, with Docker Compose or on Kubernetes, and can serve as a companion web app for a publication 📊 For differentialabundance: sample PCA and normalisation checks, DESeq2 volcano, MA and QQ plots, a genome view of the calls and GSEA gene set enrichment. ⚖️ Compared with shinyngs: the pipeline's own report is an R/Shiny application (shinyngs). The template does not replace it: it reads the same result tables and puts them (sample PCA and normalisation checks, DESeq2 volcano, MA and QQ plots, GSEA) in a dashboard that is not tied to R, sits behind one link like the dashboards of the other pipelines, and has versioning, comments and annotations. I would like to know what you use shinyngs for that this dashboard should offer too. 👀 What I need: I have built 25 templates and want people who really use these pipelines to review and criticise them. What is useful? What is misleading? 👉 Dashboard: differentialabundance dashboard (no account needed) 💬 Comment on the dashboard (Google sign-in) or in the GitHub discussion dedicated to differentialabundance. ⭐ A star on the repo is welcome if you find it useful. Thanks! 🙏 🎤 In you're interested in the project, I was speaking about Depictio at the Nextflow Summit (session page), and the recording will be available in the coming weeks.