5
$\begingroup$

I'm not sure what I'm missing, but I'm struggling to find an official hg38 GTF file with RefSeq annotations. I'd like to provide the GTF to Salmon to get gene-level annotations.

Here's Salmon's help info for --geneMap:

File containing a mapping of transcripts to genes. If this file is provided Salmon will output both quant.sf and quant.genes.sf files, where the latter contains aggregated gene-level abundance estimates. The transcript to gene mapping should be provided as either a GTF file, or a in a simple tab-delimited format where each line contains the name of a transcript and the gene to which it belongs separated by a tab. The extension of the file is used to determine how the file should be parsed. Files ending in '.gtf', '.gff' or '.gff3' are assumed to be in GTF format; files with any other extension are assumed to be in the simple format. In GTF / GFF format, the "transcript_id" is assumed to contain the transcript identifier and the "gene_id" is assumed to contain the corresponding gene identifier.

Here's what I've tried:

  1. NCBI seems to only provide GFF files. I know you can convert between the two, but that hasn't worked well for me in the past. Both file formats allow a lot of freedom, which makes conversions sloppy. Salmon, specifically, is looking for the gene_id and transcript_id annotations.
  2. UCSC provides GTF files from RefSeq, but the gene_id annotation is identical to the transcript_id annotation (i.e., it's the NM number). Or maybe there's an option I'm missing.
  3. Ensembl provides exactly the format and all of the information I could hope for, except the transcript_id is the Ensembl ID (ENST#), not RefSeq (NM#).

Can anyone tell me how to get a GTF file with hg38 RefSeq annotations?

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ refseq IDs are per-transcript, so it's expected that the gene and transcript IDs will be the same. Given that, what would you like as the gene_id? $\endgroup$
    – Devon Ryan
    Commented Sep 21, 2017 at 7:17
  • $\begingroup$ @Devon Ryan I'm guessing the EntrezGene ID or HUGO symbol? $\endgroup$
    – Ian Sudbery
    Commented Sep 21, 2017 at 14:25
  • $\begingroup$ @Mark Ebbert: What are your down stream applications? The as far as I'm aware, the salmon output doesn't include the gene_id, so it shouldn't matter that it is the same as the transcript_id. $\endgroup$
    – Ian Sudbery
    Commented Sep 21, 2017 at 14:29
  • $\begingroup$ I appreciate your help. I added Salmon's help on --geneMap. They expect gene_id to contain a "gene identifier". I understand that could be several things, but I would prefer the HUGO symbol, as @IanSudbery suggested. In my perusing, I thing Salmon should probably target the gene_name, at least how Ensembl formats the GTF. In either case, I'm surprised by @DevonRyan's comment that the gene_id would equal the transcript_id. Not arguing that's the case, but what good does that do anyone?? :-) $\endgroup$
    – Mark Ebbert
    Commented Sep 21, 2017 at 15:19
  • $\begingroup$ @IanSudbery, I ultimately want to aggregate transcript expression into the parent genes. I found last night that one of the Salmon authors now recommend tximport (github.com/COMBINE-lab/salmon/issues/98), which I will use. But this isn't the first time I've "needed" a RefSeq GTF and couldn't find one. Just surprises me, but maybe that indicates I need to assess my approach? $\endgroup$
    – Mark Ebbert
    Commented Sep 21, 2017 at 15:26

1 Answer 1

Reset to default
5
$\begingroup$

You will probably be interested in the following UCSC wiki page, which explains how to go from most of the UCSC tables to GTF/GFF:
http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format

The basic gist is that UCSC doesn't store any data internally as GTF or GFF, and so you will need to use our genePredToGtf utility in order to convert from our internally stored genePred format into GTF/GFF3.

This is the command I usually give when we get this question on our mailing list:

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select * from refGene" hg19 | cut -f2- | genePredToGtf -source=hg19.refGene.ucsc file stdin stdout

Change stdout to the output filename you want in the last command to get an hg19 refGene GTF file:

chr1    hg19.refGene.ucsc   transcript  11869   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357";  gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    11869   12227   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "1"; exon_id "NR_148357.1"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    12613   12721   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "2"; exon_id "NR_148357.2"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    13221   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "3"; exon_id "NR_148357.3"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   transcript  11874   14409   .   +   .   gene_id "DDX11L1"; transcript_id "NR_046018";  gene_name "DDX11L1";
...

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ It would be rather nice if the "GTF" format output from the UCSC table browser simply returned exactly this (including the gene entries, which are always missing). This is the sort of thing people are actually expecting from it. $\endgroup$
    – Devon Ryan
    Commented Sep 21, 2017 at 18:35
  • $\begingroup$ I agree that this would be ideal. The problem is that the gene_id's returned may not be what people want. For example, if you were trying to get a GTF output for the wgEncodeGencodeBasicV27 table, would you want a gene_id like DDX11L1 or a gene_id like ENSG*? The other issue is the plain genePred table format doesn't even have a name2 field, so what do we assign as gene_id in this case? What if the table is some other variation of genePred that also doesn't have a name2? We use many different table variations, and it is hard to support a GTF output for all of them. $\endgroup$
    – Christopher Lee
    Commented Sep 21, 2017 at 19:14
  • 1
    $\begingroup$ @ChristopherLee, that's perfect, thank you. Regarding your response to @DevonRyan's comment, I can appreciate that it's difficult (and maybe unrealistic) to support all possible variations, but I'd also argue that HUGO is probably what most people want. In either case, I'm certain that reporting the transcript ID for both gene_id and transcript_id is not ideal! :-) $\endgroup$
    – Mark Ebbert
    Commented Sep 22, 2017 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.