Entering edit mode
8.5 years ago
Parham
★
1.6k
Hi,
I need to prepare an object that holds transcripts lengths and geneIDs. In order to do that, I create a TxDb object from a gff3 file and then from TxDb I group transcripts as I show in below. The gff3 file is here.
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("../schizosaccharomyces_pombe.chr.gff3", circ_seqs = "MT", organism = "Schizosaccharomyces pombe", format = "gff3")
txsByGene <- transcriptsBy(txdb, "gene")
tail(txsByGene)
GRangesList object of length 6:
$zrt2
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] I [5317857, 5319528] + | 1487 zrt2
$zta1
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
[1] III [1794719, 1796107] - | 6779 zta1
$zuo1
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
[1] II [3094849, 3096383] - | 5286 zuo1
The problem I have is it fetches geneNames (zrt2, zta1, zuo1) while what I need is geneIDs, e.g, SPAP8A3.03 instead of zrt2.
I appreciate any advice on how to overcome this problem. I am slightly familiar with awk as well if someone can show how to do it directly on gff3 file.
Thanks in advance.
First of all thanks for a very thorough answer. I appreciate that. I ran all the options you recommended on my gff3 file. It seems that there is an issue with the gff3 file that txdb object is not created properly out of it. I have compared a txdb object created from gff3 file to a txdb created from biomaRt. Below are the outputs.
As you can see
TXNAMEandGENEIDsare what I expect and they are not gene names. Also, the length oftranscriptsis only 1 more thantranscriptsBy.Continue below:
Continued from above:
Now if I do the same with gff3 file the output is very different. Compare the
TXNAME,GENEIDand objectlengths:I have the link to gff3 file in the original question. I appreciate if you or anyone else can have a look into the file. Probably it's not only my problem.
Where does your gff file comes from? It seems that the file itself contains gene symbols in the GENEID column, perhaps that's the source of the problem.
I would do a couple of things:
1) mergeByOverlaps of the two dataframes (the one from ensembl, and your file), and check whether the TXIDs are the same for each range. Or simply subset the two dataframes on the same region, to be sure that the TXID correspond. This cannot be seen from the examples you have posted, because the second dataset is not sorted.
2) You can easily get a dataframe with the correct gene names from the ensembl dataframe, and use it to replace the GENEID column in the other (e.g. mcols(genes(txdb.ensembl, columns=c("TXID", "TXNAME", "GENEID"))) etc...
The reason why you get different lengths in the second example is because the GENEID column is not unique, which is not surprising considering that it includes gene names.
I get it from pombase . The reason I insist on using gff3 is that because I used that file for feature counting, therefore I wanted to use the same annotation file for downstream analysis. Thanks for the help.