Why are there multiple RefSeq sequences occurring with the same name?

Question

[gresserT]

I want to run some algorithms on splice site mutations. This is what I have done (or at least tried) so far:

1. Download all sequences of RefSeq from category NM_* via UCSC Table Browser in Fasta format.
2. Create index files with Samtools and bwa.
3. Read the files for my Maven Java program with the HTSJDK Library

---

As I run my program I get an Exception from HTSJDK because there are multiple RefSeq entries with the same name:

Exception in thread "main" htsjdk.samtools.SAMException: Contig 'hg19_refGene_NM_001037501' already exists in fasta index.

in the following line of my code:

FastaSequenceIndex faIndex = new FastaSequenceIndex(new File("data/RefSeqSequencesGRCh37_NM.fa.fai"));

These RefSeq entries are on both strands (+ and -) and have different positions. The sequences show some differences in the sequences, too. But usually not far from each other on the same chromosome.

---

The Questions

1. Why are there multiple RefSeq sequences with the same name?
2. Is there a way HTSJDK can handle fasta sequences with the same name?
3. Am I doing something completely wrong or inconvenient?

Answer 1

[Devon Ryan]

1. If a gene apparently has mulitple copies throughout the genome then they can get the same ID. This is a reason to use Ensembl's database, since they'll have unique IDs there.
2. I don't use the JDK (I use HTSlib), but I would guess not. If you never need to actually get the sequence of any of these then it might be worthwhile coding around this. However, I suspect that you do need the sequence, in which case it's often ambiguous what the correct sequence is.
3. No, just use Ensembl's database instead and save yourself the headaches.

[gresserT]

Since I am using ClinVar, I do only have RefSeq accession numbers.

I can't use HTSlib because it is for C and I want to extend a bigger Java project.

Shouldn't a Variant have a reference to a unique transcript?

[Devon Ryan]

BTW, since you presumably do need the sequence, the FAI files are actually pretty simple to parse and retrieve sequence from. I suspect that you'll have to write your own parser/sequence extractor that will allow iterating over all instances of a gene in the file and either report all unique sequences or all compatible sequences.

[Devon Ryan]

Ah, I had missed the Clinvar tag in your post. I have to admit that I'm not familiar enough with Clinvar to offer much guidance. Realistically, those variants were most likely originally mapped against the genome and then simply annotated with gene information. Whether those original genome mappings are on Clinvar or not I don't know, however.

Answer 2

[Kim]

This is because the UCSC RefGene track is merely an alignment track of NCBI's known RefSeq dataset, not the genome placement of those records as is provided by NCBI. It is not at all surprising that close paralogs will have more than one alignment placement.

Answer 3

[Kim]

NCBI's placement of RefSeq transcripts is available as a GFF file: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_mammalian/Homo_sapiens/reference/GCF_000001405.28_GRCh38.p2/GCF_000001405.28_GRCh38.p2_genomic.gff.gz

Genomic and transcript FASTA files are also available at this location. (Note an update is coming soon)

[Biomonika (Noolean)]

This looks very nice, but do you by any chance know about similar file with unspliced genes, not transcripts (I want introns, too)?

[gresserT]

I haven't found anything so far. Looks like getting the whole hg19 (GRCh37 or GRCh38) (from UCSC) is the way to go.

I took the whole chromosome-wise GRCh37 as Reference chromFa.tar.gz from UCSC bigZips and RefGene-File to get the Transcripts positions and the Exons positions.

Edit: I found a similar question: Convert Nm_ Mrna Position Into Corresponding Grch37 Genomic Dna Position?