This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Connect gene names and GO terms after STAR for DESeq2

EDIT: I have created a new post to address the remaining unanswered part to this question: Connect gene names to GO terms/function for loading in DESeq2

Hi, I want to ask how to go about doing differential expression analysis after getting gene counts with STAR. I would also like to connect the gene names to GO terms before loading them into DESeq2. The counts are from mouse reads mapped to the genome and annotation files from GENCODE.

I seen in another post someone's counts looking like this:

Sox17   9368    81      9287
Mrpl15  4473    30      4443
Lypla1  3122    20      3102
Tcea1   397     6       391

Mine look like this, without gene names:

ENSMUSG00000095475.1    5       0       5
ENSMUSG00000094855.1    0       0       0
ENSMUSG00000095019.1    0       0       0
ENSMUSG00000095041.7    919     16      903

I intend on making graphs similar to this (taken from the paper that performed this study already): http://cancerres.aacrjournals.org/content/canres/78/5/1334/F5.large.jpg.

Not sure what to do. I appreciate the help.

deseq2 go terms star rna-seq alignment

As far as I know, DESeq2 will take input in the form of count table like.

TranscriptORgeneId  Sample1ReadCount    Sample2ReadCount SampleNReadCount

If you want to directly give genes in DESeq2 input file instead of transcript ids then you can Vlookup GENCODE Transcript Id to their gene name (GTF file should have that information in 9th column). And after that, you can give such modified input to DESeq2.

The second solution that I would prefer is just run DESeq2 with default input (provided in the question). After getting the differential expression output to add whatever columns you want (as in annotation) from GTF use same Vlookup formula. Where your transcript Ids will be key and rest of the data will be searching array (or you can write a script instead of Vlookup).

Hello Moneeb Bajwa!

We believe that this post does not fit the main topic of this site.

I have created a new post to address the remaining unanswered part to this question: Connect gene names to GO terms/function for loading in DESeq2

For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.

If you disagree please tell us why in a reply below, we'll be happy to talk about it.

Cheers!

1 answer

Get the GTF file and run this code snippet on it:

awk 'FS="\t" {if ($3 == "gene") print $9}' in.gtf | awk -F " " '{print $2,$6}' | awk -F "; " 'OFS="\t" {print $1, $2}' | tr -d '"\;' > ensmus2name.txt

It will give you a 2-column file with both the gene_id and gene_name:

ENSMUSG00000102693.1    RP23-271O17.1
ENSMUSG00000064842.1    Gm26206
ENSMUSG00000051951.5    Xkr4
ENSMUSG00000102851.1    RP23-317L18.1
ENSMUSG00000103377.1    RP23-317L18.4
ENSMUSG00000104017.1    RP23-317L18.3
ENSMUSG00000103025.1    RP23-115I1.6
ENSMUSG00000089699.1    RP23-115I1.1
ENSMUSG00000103201.1    RP23-115I1.5
ENSMUSG00000103147.1    RP23-115I1.2

Once you have this, load it in R together with your counts. Here I assume that the above 2-column table is called ensmus2name and your count matrix is counts:

counts <- read.table("~/counts.txt", sep="\t", header = F)
ensmus2name <- read.table("~/ensmus2name.txt", sep="\t", header = F)

ENSMUS_2_NAME <- function(counts, ensmus2name){

  geneID_IDX <- sapply(as.character(counts[,1]), function(x) match(x, ensmus2name$V1))
  gene_IDs <- as.character(ensmus2name[as.numeric(geneID_IDX),2])

  return(
    cbind(data.frame(gene_IDs),
          counts[,2:ncol(counts)])
  )
}

count_with_geneName <- ENSMUS_2_NAME(counts = counts, ensmus2name = ensmus2name)

It will simply exchange the first column with the gene_ids with the correct gene_name. Should not take more than a minute or so.

EDIT: As OP mentiones in a comment below, using base::merge() probably does the same thing and is easier to use, I was not aware of that command (after two years of working with R :-D )

Thank you greatly. Regarding the GO terms and functional annotation I had mentioned, how would that come into play? The paper I am using mentions that they performed functional annotation with DAVID; my goal is to reproduce a figure similar to Figure D in the link. In the paper, the caption for that Figure D says, "Functional annotations by GO term analysis using 363 AKT common DEGs are shown (log of P values)." Not sure what to do.

EDIT: I have created a new post to address the remaining unanswered part to this question: Connect gene names to GO terms/function for loading in DESeq2

I was not aware of merge() (well, could have guessed that it exists for such a trivial problem). Anyway, merge is probably faster and easier to use, but both ways work and that is the main goal.

Log in to answer this question.