This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Why are there multiple RefSeq sequences occurring with the same name?

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?
clinvar refseq htsjdk ucsc
• 5,475 views • link
updated by Ram Ram • written by gresserT • gresserT

3 answers

  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.
• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

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?

• 0 views • link
updated by Ram Ram • written by gresserT • gresserT

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.

• 0 views • link
updated by Ram Ram • written by Devon Ryan 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.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

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.

• 0 views • link
updated by Ram Ram • written by Kim 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)

• 0 views • link
updated by Ram Ram • written by Kim Kim

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

• 0 views • link
written by Biomonika (Noolean) Biomonika (Noolean)

Log in to answer this question.

More posts like this