UTA Query for transcript is non-deterministic

[zealws]

Describe the bug

The query for transcripts against the UTA database is non-deterministic, and returns an arbitrary matching transcript each time it's run.

This is a result of a join between associated_accessions and tx_exon_aln_v which doesn't fully constrain the rows in the tx_exon_aln_v table, resulting in the database choosing one arbitrarily. The fields of this arbitrary row are used to order the results, resulting in a random transcript being selected and returned.

I ran fusor 50 times using a script, and it resulted in a variety of transcripts being selected for the 5' partner (columns are count and transcript ID):

   8 refseq:NM_001320453.1
   6 refseq:NM_001320454.1
   4 refseq:NM_001320454.2
  11 refseq:NM_001320455.1
   3 refseq:NM_001320455.2
  14 refseq:NM_006365.2
   4 refseq:NM_006365.3

Steps to reproduce

1. Download this arriba fusion.
2. Use fusor to annotate this fusion.
3. Check the 5' partner's transcript.
4. Rerun fusor, and check the transcript again.

Expected behavior

The chosen transcript is deterministic and does not change from one run to another.

Current behavior

The chosen transcript is non-deterministic, and returns an arbitrary matching transcript each time it is run.

Possible reason(s)

The query joins the tables associated_accessions with tx_exon_aln_v on associated_accessions.tx_ac = tx_exon_aln_v.tx_ac (this happens transitively through a join on the transcript with the same column).
For a given value of associated_accessions.tx_ac, there may be several matching tx_exon_aln_v rows (example below).

As a result, the join chooses an arbitrary row from tx_exon_aln_v. These rows differ in their tx_start_i and tx_end_i, which is used to order the transcripts, resulting in the transcript order varying each time the query is run.

Example query of matching values of tx_exon_aln_v for a given tx_ac:

uta=> select tx_ac, tx_start_i, tx_end_i from tx_exon_aln_v where tx_ac = 'NM_001320453.1';
     tx_ac      | tx_start_i | tx_end_i
----------------+------------+----------
 NM_001320453.1 |          0 |      173
 NM_001320453.1 |        173 |      248
 NM_001320453.1 |        248 |      344
 NM_001320453.1 |        344 |      444
 NM_001320453.1 |        444 |      794
 NM_001320453.1 |          0 |      173
 NM_001320453.1 |        173 |      248
 NM_001320453.1 |        248 |      344
 NM_001320453.1 |        344 |      444
 NM_001320453.1 |        444 |      794
(10 rows)

Suggested fix

Find a way to make the join from associated_accessions.tx_ac = tx_exon_aln_v.tx_ac deterministic.

Branch, commit, and/or version

main

Screenshots

No response

Environment details

Not relevant

Additional details

No response

Contribution

I can help prepare a PR, but I would like someone with a better understanding of the structure of the UTA database to determine how best to make the join deterministic.

[korikuzma]

@jarbesfeld Are you able to assist with this? I won't have time until next week possibly

[zealws]

After further discussion with @jsstevenson and @jarbesfeld, I think this part of the SQL query is the main problem.

Here's the whole query (logged from a fusor run), for context.

            SELECT AA.pro_ac, AA.tx_ac, ALIGN.alt_ac, T.cds_start_i
            FROM uta_20241220.associated_accessions as AA
            JOIN uta_20241220.transcript as T ON T.ac = AA.tx_ac
            JOIN uta_20241220.tx_exon_aln_v as ALIGN ON T.ac = ALIGN.tx_ac
            WHERE ALIGN.alt_aln_method = 'splign'
            AND T.hgnc = 'MIR9-1HG'
            AND ALIGN.alt_ac = 'NC_000001.11'
            
            
        ORDER BY SUBSTR(ALIGN.alt_ac, 0, position('.' in ALIGN.alt_ac)),
        CAST(SUBSTR(ALIGN.alt_ac, position('.' in ALIGN.alt_ac) + 1,
            LENGTH(ALIGN.alt_ac)) AS INT) DESC,
        ALIGN.tx_end_i - ALIGN.tx_start_i DESC;

The ORDER BY ... ALIGN.tx_end_i - ALIGN.tx_start_i DESC is the main issue.
The ALIGN table is joined as: JOIN uta_20241220.tx_exon_aln_v as ALIGN ON T.ac = ALIGN.tx_ac

I believe the ALIGN.tx_end_i - ALIGN.tx_start_i is intended to sort the returned transcripts by transcript length, but the rows in the tx_exon_aln_v view are exons, not transcripts:

uta=> select tx_ac, alt_ac, ord, tx_start_i, tx_end_i, (tx_end_i - tx_start_i) as tx_length, alt_start_i, alt_end_i, (alt_end_i - alt_start_i) as alt_length from tx_exon_aln_v where tx_ac = 'NM_001320453.1' limit 3;
     tx_ac      |    alt_ac    | ord | tx_start_i | tx_end_i | tx_length | alt_start_i | alt_end_i | alt_length
----------------+--------------+-----+------------+----------+-----------+-------------+-----------+------------
 NM_001320453.1 | NC_000001.10 |   0 |          0 |      173 |       173 |   156399167 | 156399340 |        173
 NM_001320453.1 | NC_000001.10 |   1 |        173 |      248 |        75 |   156396589 | 156396664 |         75
 NM_001320453.1 | NC_000001.10 |   2 |        248 |      344 |        96 |   156386560 | 156386656 |         96

Furthermore, there are more than 1 exon for each transcript, and so joining to the exon view by the transcript ID alone results in the database choosing a matching exon arbitrarily, and the length of that exon is taken as the transcript length to sort the transcripts.

If the purpose of the join to tx_exon_aln_v was only to get the transcript length, then I think it needs to be modified.

One approach is to use a subquery like select max(tx_end_i) where tx_ac = '%s'; which would return the end coordinate of the final exon, which should also be the sum of the lengths of all exons. I don't know if there's a different location in UTA where transcript length is recorded. But this provides a proxy for the transcript length that could be substituted into the existing query.

A different approach would be to query first for the transcripts for a gene, and then perform a second query for the exons for each transcript, and then sort the transcripts in python based on the exon information.

[zealws]

Since #439 merged, this is not user-visible in fusor, and so this is not as high-priority anymore. Downgrading the priority to medium.

[korikuzma]

@zealws I think the docstring is outdated. This is the helper method we use to get the prioritized transcripts (which gets the transcript by length)

cool-seq-tool/src/cool_seq_tool/mappers/mane_transcript.py

Line 650 in 6c8cbad

def _get_prioritized_transcripts_from_gene(self, df: pl.DataFrame) -> list:

cool-seq-tool/src/cool_seq_tool/mappers/mane_transcript.py

Lines 683 to 689 in 6c8cbad

	copy_df = copy_df.with_columns(
	copy_df.map_rows(
	lambda x: len(self.seqrepo_access.get_reference_sequence(x[1])[0])
	)
	.to_series()
	.alias("len_of_tx")
	)

[zealws]

In that case, I'm confused about the desired behavior for UTADatabase.get_transcripts?

[korikuzma]

This method is intended to get transcript info for a given gene or genomic accession on a given position, if provided.

We might want to change the ordering behavior for the transcript part to follow what's done for the genomic accession (and update the docs accordingly).

ORDER BY 
  SUBSTR(ALIGN.alt_ac, 0, position('.' in ALIGN.alt_ac)),
  CAST(SUBSTR(ALIGN.alt_ac, position('.' in ALIGN.alt_ac) + 1, LENGTH(ALIGN.alt_ac)) AS INT) DESC,
  SUBSTR(ALIGN.tx_ac, 0, position('.' in ALIGN.tx_ac)),
  CAST(SUBSTR(ALIGN.tx_ac, position('.' in ALIGN.tx_ac) + 1, LENGTH(ALIGN.tx_ac)) AS INT) DESC;