Hi bm.nina ,

OK first things first. For DESeq2, you want properly summarized gene-level counts, not TPMs. If using featureCounts, a safer approach is to convert the ITAG4 GFF3 to a clean GTF, confirm that exon rows contain gene_id, and then count exons grouped by gene_id. For example, after conversion: if your real goal is TPMs rather than differential expression, bypass DESeq2 and use Salmon, kallisto, or another transcript-quantification workflow. Salmon outputs TPM, effective length, and estimated read counts; kallisto also estimates transcript abundance. If the final goal is differential expression from Salmon/kallisto output, import estimated counts with tximport and then run DESeq2; do not give TPMs to DESeq2 as counts.

Beyond that there are a couple issues, I think. I would separate three steps here: STAR alignment, feature counting, and DESeq2. DESeq2 is neither a splice-site discovery tool nor a TPM-generation tool: it expects a feature-by-sample matrix of un-normalized or estimated counts for differential expression. If you are seeing “14 splice sites” per gene in something you are passing to DESeq2, that structure probably came from the annotation/counting/import step, not from DESeq2 itself. STAR’s splice-junction output and gene-count output are different things. SJ.out.tab reports splice junctions, while ReadsPerGene.out.tab from --quantMode GeneCounts reports gene-level counts with columns for unstranded and stranded interpretations. So I would debug the “413 splice sites” and the “14 splice sites per gene” issues separately.

Also,

featureCounts -t exon -g gene_id -s 0 -a ITAG4.0_gene_models.gtf -o counts.txt *.bam

Set -s according to strandedness: 0 unstranded, 1 stranded, 2 reverse-stranded. Since these are single-end reads, do not use paired-end options.

GTF files usually carry gene_id and transcript_id attributes on exon rows. GFF3 uses a hierarchy with ID and Parent, and in many GFF3 files the Parent of an exon is a transcript, not the gene. Therefore, a counting tool can easily produce transcript- or exon-level rows if the grouping attribute is wrong.³ STAR can use GFF3 for splice-junction annotation during genome generation, for example with:

--sjdbGTFfile annotation.gff3 \
--sjdbGTFfeatureExon exon \
--sjdbGTFtagExonParentTranscript Parent \
--sjdbOverhang 53

For 54 bp reads, --sjdbOverhang 53 follows STAR’s read-length-minus-one recommendation. This helps STAR interpret exon-to-transcript relationships, but it does not guarantee that downstream featureCounts grouping is gene-level.