Exon-level orthology

[zkliesmete]

Thank you for the great tool and annotations!

I am interested in the exon-level orthology between species, however I could not find how the exons correspond. I am using the ferret files from
https://genome.senckenberg.de//download/TOGA2/TOGA2/reference_human_hg38/Mustela_putorius_furo__domestic_ferret__HLmusPutFur2__GCF_011764305.1/

I read in the wiki documents that the reference (human) exon IDs are in the exon_aln.fa.gz file, however if I extract the coordinates from the sequence names and intersect those with the coordinates in query_annotation.gtf.gz , I only get around 2.5k intersecting regions. Why is that and how do I gather these correspondences?

Also a side question - which hg38 genome annotation version (gtf) was used for creating these files exactly?

Thanks in advance!

[ReverendCasy]

Hello @zkliesmete ,
Sorry for the delayed reply. Could you please:

1. Specify how you extract exon coordinates from exon_aln.fa.gz and which exon IDs you use;
2. Specify how you intersect the resulting coordinates to those in the GTF file.

So far I’m a bit lost whether you’re looking for reference or query exon coordinates. I also tried assessing the total exon correspondence myself and found the number to be much higher than what you report:

zcat hg38/TOGA2/vs_HLmusPutFur2/exon_aln.fa.gz | \
    sed 's/\s//g' | \
    awk -F'|' -vOFS='\t' '($10 == "query_exon"){split($4,a,":"); split(a[2],b,"-"); if(b[1] == "" || b[2] == ""){next}; print a[1],b[1],b[2]}' | \
    sort -k1,1 -k2,2n > from_fa.bed3
zcat  hg38/TOGA2/vs_HLmusPutFur2/query_annotation.gtf.gz | \
    awk -F'\t' '($3 == "exon" || $3 == "CDS"){print $1"\t"$4"\t"$5}' | \
    sort -k1,1 -k2,2n > from_gtf.bed3
overlapSelect -inFmt=bed -selectFmt=bed  -overlapThreshold=1.0 from_fa.bed3 from_gtf.bed3 stdout | wc -l
# 1260977

Also, it might be much better to extract the exon coordinates from meta/exon_meta.tsv.gz . The only caveat is that the file lists exon data with regard to reference exon, so precise intron deletions either in the query or in the reference will not be reflected in the file.

Also a side question - which hg38 genome annotation version (gtf) was used for creating these files exactly?

We used GCF_000001405.38 + Ensembl + MANE .

Best,
Yury

[zkliesmete]

Hi Yury,

Sorry that my message was confusing. So essentially I am looking for the orthologous exon annotation across species to do alternative splicing analyses later on. I tried your way but it does not seem to give me something like:

transcript_id exon_id_ferret coord_ferret exon_id_human coord_human

Could you be so kind and indicate how I get that?
Also, where is this custom human reference gtf stored (GCF_000001405.38 + Ensembl + MANE)? I would like to map the human data to it.

Many many thanks!
Best,
Zane

[ReverendCasy]

I see. TOGA2 does not produce any exon level orthology files by default. Assuming you’re interested in 1:1 orthologs, consider something like the following:

awk -F’\t’ ‘($5 == “one2one”){print $4}’ orthology_classification.tsv | \
    zgrep -f- -Fw meta/exon_meta.tsv.gz | awk -F’\t’ -vOFS=“\t” ‘{if ($8 == “I"){print $4,$5,$6,$1,$2,$7}}'

This will give you a BED6 file of all non-missing, non-deleted exons from the 1:1 orthologs annotated by TOGA2 in the selected query, exon numbers stored in the score column. If you want to get orthologous exons from a larger array of queries, consider using toga2.py orthogroups beforehand.

Also, where is this custom human reference gtf stored (GCF_000001405.38 + Ensembl + MANE)? I would like to map the human data to it.

See reference_annotations/hg38 in this repo.

[zkliesmete]

Thank you! I understand that this gives me the exon numbers from human in the ferret genome space.

However, if I intersect the exact coordinates in exon_meta_tsv.gz with those in query_annotation.gtf, I only get a few thousand non-unique correspondences, so I still do not see a simple solution to this exon_x_human == exon_y_ferret orthologue search that I am trying to do..

Should I then just match the exons by the largest mutual overlap in the ferret genome space between exon_meta_tsv.gz and query_annotation.gtf?

Thank you again
Zane

[ReverendCasy]

However, if I intersect the exact coordinates in exon_meta_tsv.gz with those in query_annotation.gtf, I only get a few thousand non-unique correspondences

Oh. Could you please share how you intersect these two files?

[zkliesmete]

sure, I used R here - I hope it's fine

exon_meta<-readr::read_tsv("exon_meta.tsv.gz")
gtf_ferret<-plyranges::read_gff("query_annotation.gtf") %>% as_tibble()
exact_matches<-gtf_ferret %>% dplyr::inner_join(exon_meta, by= c("seqnames" = "chrom", "start", "end", "strand")) #2593 many-to-many matches

Thanks for your help!

[zkliesmete]

Hi Yuri,

Did you get to check out this coordinate mismatch?
Also, if you can point me to a script where you combine and filter the reference gtf (GCF_000001405.38 + Ensembl + MANE), that would be great.

Thanks
Zane

[ReverendCasy]

Hi @zkliesmete ,
Sorry for the delay. I will try to get my hands on the issue in the following days.

Also, if you can point me to a script where you combine and filter the reference gtf (GCF_000001405.38 + Ensembl + MANE), that would be great.

The scripts used to create the input annotation are not yet on GitHub; I will ping the colleagues to share them. The procedure is described in the preprint in "Materials and Methods" -> "Evolutionarily informed transcript selection for input annotations” subsection.