I am working with DeSeq2 and trying to convert ENSG IDs to Gene symbols after calling the results function in DeSeq2. When doing this however, I end up with duplicated (or sometimes triplicated or more) gene symbols each corresponding to a different ENSG ID. These duplicated symbols often have the same exact P-adj, logFC values, etc but sometimes they are slightly different.

I am a novice at RNA-Seq analysis so apologies if this is a silly question! I have tried looking on bioconductor, DeSeq vignette, etc and have not been able to figure out a reason or work around for this issue.

An example of one of the genes this is happening with is: AADACL2 which shows up with two ENSG IDs, ENSG00000197953 and ENSG00000261846.

My code is (hopefully) posted here:

#Creating a DeSeq Matrix from Tximport

samples <- read.table(file.path(dir,"samples_v3.txt"), header=TRUE)
files <- file.path(dir,"salmon", paste("salmon_",samples$run, ".txt", sep='') )
load("/Users/rimlandca/Desktop/R files/tx2gene_GRCh38.80.RData")
txi <- tximport(files, type="salmon", tx2gene=tx2gene)
ddsTxi <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ type)

#Collapsing Technical Replicates

ddsColl<- collapseReplicates(ddsTxi, samples$name)

#Pre-Filter

ddsColl <- ddsColl[ rowSums(counts(ddsColl)) > 1, ]

#Subset for protein coding genes (we are only interested in protein coding genes)

ddsCollinfo <- data.frame(getBM(attributes=c("ensembl_gene_id","hgnc_symbol", "transcript_biotype"), filters="ensembl_gene_id",values=rownames(ddsColl), mart=ensembl))

ddsCollinfocoding <- subset(ddsCollinfo, ddsCollinfo$transcript_biotype == "protein_coding")

ddsCollcoding <- subset(ddsColl, rownames(ddsColl) %in% ddsCollinfocoding$ensembl_gene_id)

#Run DeSeq on just a subset of the biological groups

ddsCollcoding_CBDonly<-ddsCollcoding[, ddsCollcoding$origin == "CBD"] 

ddsCollcoding_CBDonly$type<-droplevels(ddsCollcoding_CBDonly$type) 

dds_DGE_CBD <- DESeq(ddsCollcoding_CBDonly)

resultsCtCc <- results(dds_DGE_CBD, contrast=c("type", "CBD_t", "CBD_c"), alpha=0.05) 

summary(resultsCtCc)

#Add Gene Symbols to Results output and export

library( org.Hs.eg.db ) 
library(AnnotationDbi) 
resultsCtCc$symbol <- mapIds(org.Hs.eg.db, keys=row.names(resultsCtCc), column="SYMBOL", keytype="ENSEMBL", multiVals="first")

resultsCtCc.df <- as.data.frame(resultsCtCc)

write.csv(resultsGtGc.df, file="/Users/rimlandca/Desktop/R Files/resultsGtGc.csv")​

You will often find things like this, because the human genome is a complicated place and simple things like a gene symbol are not sufficient to uniquely identify a gene. Using the example you give, we can explore things.

> select(org.Hs.eg.db, c("ENSG00000197953","ENSG00000261846"), c("SYMBOL","GENENAME","ENTREZID"),"ENSEMBL")
'select()' returned 1:1 mapping between keys and columns
          ENSEMBL  SYMBOL                         GENENAME ENTREZID
1 ENSG00000197953 AADACL2 arylacetamide deacetylase like 2   344752
2 ENSG00000261846 AADACL2 arylacetamide deacetylase like 2   344752

So it appears that NCBI thinks this is one single gene, and EBI thinks it's two. This is a consistent issue - different people think different things, and NCBI and EBI are different groups that don't always agree on things. For example, if we use the Homo.sapiens package to see where this gene lives, we get this:

> select(Homo.sapiens, c("ENSG00000197953","ENSG00000261846"), c("CDSCHROM","CDSSTART","CDSEND"), "ENSEMBL")
'select()' returned 1:many mapping between keys and columns
           ENSEMBL CDSCHROM  CDSSTART    CDSEND
1  ENSG00000197953     chr3 151451824 151451961
2  ENSG00000197953     chr3 151458434 151458656
3  ENSG00000197953     chr3 151461881 151461950
4  ENSG00000197953     chr3 151463297 151463468
5  ENSG00000197953     chr3 151474780 151475382
6  ENSG00000197953     chr3 151474816 151475382
7  ENSG00000261846     chr3 151451824 151451961
8  ENSG00000261846     chr3 151458434 151458656
9  ENSG00000261846     chr3 151461881 151461950
10 ENSG00000261846     chr3 151463297 151463468
11 ENSG00000261846     chr3 151474780 151475382
12 ENSG00000261846     chr3 151474816 151475382

But if we use an EnsDb package (I just have an old one - you would do well to get an updated version), we get this:

> select(EnsDb.Hsapiens.v79, c("ENSG00000197953","ENSG00000261846"), c("GENEID","SEQNAME","GENESEQSTART","GENESEQEND"), "GENEID")
           GENEID             SEQNAME GENESEQSTART GENESEQEND
1 ENSG00000197953                   3    151733916  151761339
2 ENSG00000261846 CHR_HSCHR3_1_CTG2_1    151744454  151770036

So EBI/Ensembl think this gene lives on an unplaced scaffold or haplotype or some such thing as well. In which case they evidently want to give it a different Ensembl ID. But if the sequence is identical, you may get most or all of your reads aligning to both places, in which case you are likely to get very similar fold changes and p-values.