Hi #C0493FQQZ3P I am currently performing RNA-Seq analysis using the nf-core/`rnasplice` pipeline with the

STAR_SALMON

aligner-quantification method, and I have a question regarding the use of featureCounts output for gene-level analysis. The pipeline appears to support quantification using

featureCounts

, which is particularly useful for differential exon usage analysis with

edgeR

when the

--edger_exon

parameter is enabled. In this mode, the feature type is set to exon and cannot be changed. During my run, I also specified

--gtf_group_features

gene_id to group features by gene. I understand that the pipeline also generates a tximport folder containing various types of gene and transcript-level quantifications derived from

Salmon

. However, my collaborator has specifically requested that I perform gene-level analysis using the

featureCounts

-derived data. I’m wondering if it is appropriate to use the exon-level counts generated by the pipeline and aggregate (sum) them by gene_id to obtain gene-level counts. Would this approach be valid for downstream analysis using

edgeR

or

DESeq2

? Alternatively, would it be better to run

featureCounts

independently outside of the rnasplice pipeline to generate gene-level counts directly? Any guidance on best practices for this would be greatly appreciated.

[toufim@sumner-cba003 featurecounts]$ head -25 DD140_progenitors_B2.featureCounts.txt
# Program:featureCounts v2.0.1; Command:"featureCounts" "-F" "SAF" "-t" "exon" "-g" "gene_id" "-f" "-O" "-B" "-C" "-p" "-T" "6" "-a" "annotation.saf" "-s" "1" "-o" "DD140_progenitors_B2.featureCounts.txt" "DD140_progenitors_B2_sorted.bam" 
Geneid	Chr	Start	End	Strand	Length	DD140_progenitors_B2_sorted.bam
ENSG00000000003.16	chrX	100627108	100627108	-	1	4
ENSG00000000003.16	chrX	100627109	100629986	-	2878	1860
ENSG00000000003.16	chrX	100630759	100630866	-	108	880
ENSG00000000003.16	chrX	100632485	100632540	-	56	705
ENSG00000000003.16	chrX	100632541	100632568	-	28	643
ENSG00000000003.16	chrX	100633405	100633441	-	37	700
ENSG00000000003.16	chrX	100633442	100633539	-	98	882
ENSG00000000003.16	chrX	100633931	100634029	-	99	1164
ENSG00000000003.16	chrX	100635178	100635252	-	75	1194
ENSG00000000003.16	chrX	100635558	100635746	-	189	1417
ENSG00000000003.16	chrX	100636191	100636607	-	417	26
ENSG00000000003.16	chrX	100636608	100636689	-	82	834
ENSG00000000003.16	chrX	100636690	100636792	-	103	645
ENSG00000000003.16	chrX	100636793	100636806	-	14	96
ENSG00000000003.16	chrX	100636807	100637104	-	298	51
ENSG00000000003.16	chrX	100639945	100639991	-	47	5
ENSG00000000005.6	chrX	100584936	100585066	+	131	0
ENSG00000000005.6	chrX	100585231	100585362	+	132	0
ENSG00000000005.6	chrX	100593624	100593894	+	271	0
ENSG00000000005.6	chrX	100593895	100594035	+	141	0

featureCounts

is a fairly simple tool, and if confronted with multimapped reads, it does not know how to resolve the ambiguity and just throws them out. That is a good reason to use newer/smarter tools like Salmon or RSEM. But if you must use

featureCounts

, run it independently of the pipeline (using the BAMs) at the gene level, because if try adding up the exons, you might have alternate splice forms with overlapping ones (e.g. alternate 3'- or 5'-splice site), and a read falling into them would probably not be counted, since the exon would be ambiguous, even though the gene is clearly not. Once you have a counts table, you can use DESeq2 or edgeR. Note that this pipeline is mostly concerned with finding and quantifying alternative splice forms or transcripts, and gene expression is mostly a concern for normalization. If your primary purpose it gene expression analysis, there is #CE8SSJV3N.

sure, thanks much.