Status_case_vs_ctr_results.txt outputs GeneIDs instead of gene names

[smoudour]

Description of the bug

I'm running a simple analysis on 4 samples (2 case replicates and 2 ctr replicates), just the DESeq2 part. When checking the results, the produced .txt file along with the plots show the geneIDS instead of Gene names. On the contrary, when I run the test pipeline, I get the gene names. I couldn't find any argument specifying this detail, but I'm not completely sure.

After checking a bit more, I find that software version files are pretty much identical. I'm suspecting that the root of this is the structure of tx2gene.tsv file output by the nf-core/rnaseq pipeline: while in the test_dataset folder (https://github.com/lconde-ucl/DGE2/blob/1c8d3da5ee82c04d121c6a8d78c4367c4d9e76de/assets/test_datasets/results_rnaseq/star_salmon/tx2gene.tsv) it looks like this:

rna0	DDX11L1	DDX11L1
rna1	WASH7P	WASH7P
...

while mine looks different:

transcript_id   gene_id gene_name
ENST00000511072 ENSG00000142611 PRDM16
ENST00000607632 ENSG00000142611 PRDM16
...

I'm guessing that when DESeq2 is run in R, it picks the second collumn instead of the third one, that's why I get the gene_ids instead of gene_names as output in the final results text file.

/

Command used and terminal output

nextflow run -bg lconde-ucl/DGE2 -profile docker -params-file params.yaml

Relevant files

params.zip

emtab9014_metadata_1.txt

System information

CREATE_PARAM_FILE:
bash: 5.0.17 3
CUSTOM_DUMPSOFTWAREVERSIONS:
python: 3.11.7
yaml: 5.4.1
DESEQ2:
deseq2 version: 1.42.0
r_version: R version 4.3.2 (2023-10-31)
REPORT_DESIGN:
ComplexHeatmap version: 2.18.0
DESeq2 version: 1.42.0
ReportingTools version: 2.42.2
ashr version: 2.2.63
ggplot2 version: 3.4.4
hwriter version: 1.3.2.1
r_version: R version 4.3.2 (2023-10-31)
Workflow:
Nextflow: 24.10.5
lconde-ucl/DGE2: '1.0'

[lconde-ucl]

Hi @smoudour,

This is not a bug, but the expected behaviour. The pipeline uses txtimport to import the salmon data, and tximport uses the second column of the txt2gene file as the gene ID (if there is a third column, that one is ignored)

If you run the nf-core/rnaseq pipeline, it will output txt2gene file with 3 columns where the gene_id and the gene_name is picked from assembly (specifically, the GTF file) that you used. For example, these are the GTF files for the Ensembl and NCBI assemblies in iGenomes:

==> Ensembl/GRCh37/Annotation/Genes/genes.gtf <==
1	processed_transcript	exon	11869	12227	.	+	.	exon_id "ENSE00002234944"; exon_number "1"; gene_biotype "pseudogene"; gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; transcript_id "ENST00000456328"; transcript_name "DDX11L1-002"; transcript_source "havana"; tss_id "TSS15133";
1	processed_transcript	transcript	11869	14409	.	+	.	gene_biotype "pseudogene"; gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; transcript_id "ENST00000456328"; transcript_name "DDX11L1-002"; transcript_source "havana"; tss_id "TSS15133";
1	transcribed_unprocessed_pseudogene	exon	11872	12227	.	+	.	exon_id "ENSE00002234632"; exon_number "1"; gene_biotype "pseudogene"; gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; transcript_id "ENST00000515242"; transcript_name "DDX11L1-201"; transcript_source "ensembl"; tss_id "TSS178604";
1	transcribed_unprocessed_pseudogene	transcript	11872	14412	.	+	.	gene_biotype "pseudogene"; gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_source "ensembl_havana"; transcript_id "ENST00000515242"; transcript_name "DDX11L1-201"; transcript_source "ensembl"; tss_id "TSS178604";

==> NCBI/GRCh38/Annotation/Genes/genes.gtf <==
chr1	BestRefSeq	exon	11874	12227	.	+	.	gene_id "DDX11L1"; gene_name "DDX11L1"; transcript_id "rna0"; tss_id "TSS31672";
chr1	BestRefSeq	exon	12613	12721	.	+	.	gene_id "DDX11L1"; gene_name "DDX11L1"; transcript_id "rna0"; tss_id "TSS31672";
chr1	BestRefSeq	exon	13221	14409	.	+	.	gene_id "DDX11L1"; gene_name "DDX11L1"; transcript_id "rna0"; tss_id "TSS31672";
chr1	BestRefSeq	exon	14362	14829	.	-	.	gene_id "WASH7P"; gene_name "WASH7P"; transcript_id "rna1"; tss_id "TSS15911";

As you can see, Ensembl uses Ensembl IDs as gene_id, while NCBI uses the standard HGNC gene name, so these are the ones that will be output as the second column in the txt2gene file. If you want the DE results to contain gene names rather than ensembl IDs, you’ll need a tx2gene.tsv file where the second column contains gene names. You can either rerun the nf-core/rnaseq pipeline using an assembly with gene_id = gene names (e.g. NCBI instead of Ensembl), or modify your existing tx2gene.tsv file so that the gene names are in the second column .

Best
Lucia

[smoudour]

Hi again Lucia,
Thanks for the response, I get it, I will need to check some more with rnaseq because currently I'm not able to locate GTF files from NCBI, only GFF3. I think rnaseq does perform a transformation from GFF3 -> GTF.

A kind of irrelevant question, can I create a metadata file that will point to multiple desired comparisons? For example, if I have a metadata.txt:

SampleID        Status
Control_1       ctr
Control_2       ctr
PLX4720_1       case1
PLX4720_2       case1
PLX7904_1       case2
PLX7904_2       case2
PLX8394_1       case3
PLX8394_2       case3

and I want to perform 3 binary comparisons (ctr vs each case), is that possible by using --treatment argument (like passing a list of comma separated treatments)?
I want to specifically define --design --treatment etc. because if not, the comparisons are done in reverse (case vs ctr / ctr vs case) and I get reversed signs on the logfold changes.

[lconde-ucl]

@smoudour yes, you can specify the comparison of interest by using the --design, --condition, --treatment, and --control parameters (please see the usage doc for an example). If you want 3 different comparisons, you will need to submit the pipeline 3 times, one per comparison.

By default, if you don't specify a design, the pipeline will run all possible comparisons, and the levels in each condition (i.e., the treatment and control groups) are picked alphabetically. If these are not in the order you want, you could potentially rename the groups, for example "a_ctr" instead of "ctr". But in your case (and in general) I suggest you run the pipeline specifying the design of interest as above.

[smoudour]

@lconde-ucl What I was thinking is performing all desired binary comparisons (in my case ctr vs case1, ctr vs case2 & ctr vs case3) in the same pipeline run (one command), but from your answer I am deducting that the only way to do this is submitting the pipeline 3 times, each time with different --treatment and --control arguments. Thanks a lot for the help!

[lconde-ucl]

@smoudour To clarify, you can run all comparisons in a single pipeline execution if you don’t specify a design. However, in that case, the treatment/control groups will be chosen alphabetically (which is why you’re seeing comparisons "in reverse"), and you’ll also get other comparisons that are likely not of interest (e.g., case1 vs case2).