Hello,
Background:
I am analyzing data for an RNA-seq experiment made on a bacterial strain. I have a previously assembled genome that was annotated with PGAP. I processed my raw RNA-seq reads, aligned them to the genome with STAR, and quantified the alignments with HTseq. I fed the HTseq read counts to DESeq2 and now I have a table with the results, though I end up with a gene_id in the first column of the table and DESeq2 results.
Problem:
The gene_id is not very informative. I want the gene_id, the CDS start and end position (to control for gene length bias), and also the transcript product (protein) or function. I have this information in a GTF annotation file of course file but I can't find any sort of tools or code to add this information to my DESeq2 results.
I thought it would be a simple task, but either it is not, or I am not understanding something properly. I have seen some tutorials using Biomart or other tools like that, but it seems like it's just for eukaryotes and it's using online databases. I guess I don't understand why one would not simply use the GTF annotation file?
My plan going forward is to annotate with GO terms or KEGG to do an enrichment analysis. The step I am stuck on seems easy but for some reason, I can't find a solution.
Any help is welcome, Thank you.
•
written
by
Miguel •
2 answers
I'd say this is more a question for Bioconductor Support https://support.bioconductor.org/ but I am going to answer this on a general basis. Yes, it is a simple task if simplicity is defined using existing packages in R. The only difficulty is to know the right package to use. A search for "R Import GTF file" yields https://www.biostars.org/p/272889/ as top hit. That gives you your answer for how to import a GTF file.
Once you have rtracklayer installed, check the help on width, start, and end. Using R function merge you can join your GTF annotation data or any other annotation with the DEseq results as long as the annotation also contains the same gene ids. Described here: https://stackoverflow.com/questions/1299871/how-to-join-merge-data-frames-inner-outer-left-right
In addition, you might need to use as.data.frame to convert the classes to generic data.frames. This function works on most things in R (with some mixed results), and works very well in this case.
Log in to answer this question.
-
How to load a galaxy DESeq results table into R so I can continue my workflow there
written by Nicolas •
-
featureCounts: ERROR: failed to find the gene identifier attribute
written by v.berriosfarias •
-
Weird looking dispersion & MA plot?
written by Zygpy •
-
Nanopore sequencing bacteria and problems with NCBI PGAP annotation
written by GGG_Alex •
-
BadProteinStart errors when uploading PGAP .sqn file to NCBI
written by yarmda •
-
Setting Up Model in DeSeq2
written by lkw222 •
-
HTSeq errors in the command line
written by Morris_Chair
-
Building RNAseq pipeline STAR-RSEM-tximport-DESeq2
written by jkim •
-
Obtaining a GTF with a required set of gene IDs
written by spacer •
-
BEDOPS gtf2bed conversion error with Ensembl GTF
written by bioinformatics.cancer
This information should be in the GTF/GFF file that you used. Can you show us a few lines of that file?
Yes, it is in that file. I just did not know how to join it with DESeq2 results. I will try the solution from post below and see if it works.
Here's a line:
I don't know how you did the read counting but
featureCountscan assign the gene name using the ID you specify. You should have usedNameas the id.I used HTseq which is integrated into STAR. Perhaps I could try featureCounts yes.
With
htseq-countyou could have used-i Nameor--additional-attrto add a column with names.