Address RefSeq transcript misalignments

[holtgrewe]

RefSeq transcripts can align with indels and mismatches to the reference sequence. While mismatches could be argued to be non-critical (assuming the GenBank entries that the RefSeq transcript is based on is from healthy individuals), indels cannot.

For hg19, 884 transcripts in 501 genes are affected.

The following solution will be implemented:

• The default_sources.ini file gets a settings "fixIndels" and "fixIndelsUcsc".
• When parsing the RefSeq transcript database, the Note attribute is analyzed.
  If it contains the substrings "indel" or "substitution" then this is recorded into the built TranscriptModel.
• When fixIndels=true is given then the user also has to provide the path to the reference sequence.
• The file at fixIndelsUcsc is used for providing the UCSC transcript alignments.
  This will be used for the exon and CDS information.
  The sequence will be taken from the reference.

NB: This will create an incompatibility between the databases built before and after Jannovar v0.29.

For each hg*/refseq* entry, a _fixindel variant is added that contains these fix transcripts. This way, the fixed transcripts are strictly opt-in and only supplement those where the indel is not fixed. Variants for both can be reported.

[holtgrewe]

This approach for the resolution will work for now. A proper resolution will be done in #450.

@pnrobinson @julesjacobsen as ENSEMBL is not affected this will probably not be of super big interest to you

@roland-ewald I think this might be important for you to know

A good example is LTBP4:

• UCSC
• ENSEMBL

[roland-ewald]

@holtgrewe Thank you for the heads-up, I'm aware of the issue. BTW besides the Note attribute there is also the Gap attribute in the GFF file which gives the exact position and structure of the mismatch between reference genome and transcript (for example, Gap=M123 I1 M456 means after 123 matching nt the GFF region has one additional nucleotide that is not in the reference genome). It might be safer to use this attribute for 'indel-detection' instead of relying on the Note attribute.

[holtgrewe]

ohhhh! thanks for this info

[roland-ewald]

Sure, you're welcome!

[holtgrewe]

After looking into this some more, I think the correct way is to use the alignments from the RefSeq file that are stored separately from the exons themselves. It is important to both capture indels in the local alignments themselves (the Gap attribute) and the boundaries of the alignments in the Target attribute.

For example, NM_014747.2 does not have any record with a Gap attribute but the first alignment is Target=NM_014747.2 5 263 +, that is 4 bases need to be clipped from the transcript sequence.