Resolve mismatch between ProteinGym and Haplosaurus sequences

[bschilder]

Haplosaurus provides each haplotype as an unaligned sequence, and as an sequence aligned with the reference genome.

I'm using the aligned sequences as a convenient way to grab the right coordinates and insert ClinVar mutations, using the biopython MSA object class.

However, my checks indicate that the ref sequence is often far shorter than the ref sequence provided by ProteinGym in their clinical sub/indel files.

This could be due to:

1. A problem with Haplosaurus (eg truncating ref sequences when the haplotype sequence is truncated).
2. A problem with ProteinGym (e.g. the ref sequence is wrong, or the mutation data is wrong)
3. A mismatch in the reference genome being used to get the sequences between Haplosaurus and ProteinGym.
4. A problem with my gprofiler submodule and how it's handling mapping.

[bschilder]

Potential solutions:

• Create a new pairwise alignment between the ProteinGym ref sequence and each Haplosaurus haplotype.
• Create a single MSA between the ref sequence (from ProteinGym, or ESNP**:REF seq, or Ensembl) and all haplotypes from that protein.

Resources:

• https://biopython.org/docs/1.76/api/Bio.Align.html#Bio.Align.PairwiseAligner
• https://biopython.org/docs/dev/Tutorial/chapter_msa.html
• https://biopython.org/docs/dev/Tutorial/chapter_msa.html#alignment-tools

[bschilder]

• Add unit tests to check whether ESNP**:REF seq is the same as all the aligned seqs

[bschilder]

Added some unit tests to check this
https://github.com/bschilder/VEP_protein/blob/main/tests/get_haplotype_ref.py

Using a sample of 503 transcripts, all of the haplotype ref sequences (from the Haplosaurus alignment) have the expected length (as determined by the ENP**:REF unaligned sequence).

So that eliminates Possibility 1.

[bschilder]

Reached out to ProteinGym team to get NP:ENSP mappings directly from them:
OATML-Markslab/ProteinGym#72

[bschilder]

@mengysun do you have any insight as to how I can better map the ProteinGym NP IDs to Haplosaurus ENP IDs?

[mengysun]

Maybe this can be obtain from sensembl biomart (or their FTP site): https://useast.ensembl.org/Help/Faq?id=294 . If only protein level mapping is needed, one way is through ncbi gene ftp site: https://ftp.ncbi.nih.gov/gene/README , but the mapping probably will be less direct, you will need to leverage multiple tables together, which can have missing data.

Do pay attention to the version number and the difference between transcript id and protein id.

[bschilder]

Thanks @mengysun

Maybe this can be obtain from sensembl biomart (or their FTP site): https://useast.ensembl.org/Help/Faq?id=294 . If only protein level mapping is needed, one way is through ncbi gene ftp site: https://ftp.ncbi.nih.gov/gene/README , but the mapping probably will be less direct, you will need to leverage multiple tables together, which can have missing data.

That's more or less what I'm doing currently, namely using the gprofiler API, which pulls from Ensembl:
https://biit.cs.ut.ee/gprofiler/page/apis

The problem is, for every NP ID there are many corresponding ENSP IDs (1-3000 in some cases).

I'll try using the Ensembl API directly and see if that helps.

Do pay attention to the version number and the difference between transcript id and protein id.

Good point, I think this may help reduce the number of mappings. For some reason gprofiler doesn't provide the version.