Unable to index fasta file

[guillaumecharbonnier]

I am unable to get an index for any fasta I try. Here is below what I get if I try IKZF1 as a test:

out/wget/https/github.com/Oshlack/Toblerone/releases/download/v0.0.8/tinyt_amd64 index -i out/toblerone/build_index_fa-genome-hg38_IKZF1/toblerone_transcriptome.tidx out/toblerone/build_index_fa-genome-hg38_IKZF1/toblerone_transcriptome.fa
 2024-07-01T17:59:29.505 INFO  tinyt > Building index from fasta: out/toblerone/build_index_fa-genome-hg38_IKZF1/toblerone_transcriptome.fa
 2024-07-01T17:59:29.511 INFO  tinyt::utils > Reading transcripts from Fasta file
thread 'main' panicked at 'called `Option::unwrap()` on a `None` value', src/utils.rs:115:34
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

I tried to use both v0.0.8 and 0.0.9 with the same issue.
Here is attached the fasta I try to index.

toblerone_transcriptome_IKZF1.zip

Can you share a working fasta so I can understand what is wrong with mine?

[guillaumecharbonnier]

I have inspected https://github.com/Oshlack/Toblerone/blob/master/src/utils.rs and found the fasta header was not in the expected format. I have suffixed each header with a " gene=IKZF1" and the error disappeared.

However, I get no assigned counts with my tidx whereas I get some counts if I use the IKZF1.tidx file in the github repo.

==> test_sample_VS_my_IKZF1_tidx.csv <==
Gene	Deletion	Count	Total	GeneLength	ReadLength	ScaleFactor	Proportion	ScaledProportion
IKZF1	ENST00000331340.8_del2	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del2_3	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del2_3_4	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del2_3_4_5	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del2_3_4_5_6	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del2_3_4_5_6_7	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del3	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del3_4	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del3_4_5	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del3_4_5_6	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del3_4_5_6_7	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del4	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del4_5	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del4_5_6	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del4_5_6_7	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del5	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del5_6	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del5_6_7	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del6	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del6_7	0	1483	6255	0	0	0	NaN
IKZF1	ENST00000331340.8_del7	0	1483	6255	0	0	0	NaN

==> test_sample_VS_Toblerone_IKZF1_tidx.csv <==
Gene,Deletion,Count,Total,GeneLength,ReadLength,ScaleFactor,Proportion,ScaledProportion
Gene	Deletion	Count	Total	GeneLength	ReadLength	ScaleFactor	Proportion	ScaledProportion
2	ENST00000331340.8_del2	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del2_3	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del2_3_4	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del2_3_4_5	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del2_3_4_5_6	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del2_3_4_5_6_7	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del3	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del3_4	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del3_4_5	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del3_4_5_6	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del3_4_5_6_7	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del4	8	1483	6255	0	0	0.005394471	inf
2	ENST00000331340.8_del4_5	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del4_5_6	1	1483	6255	0	0	0.000674309	inf
2	ENST00000331340.8_del4_5_6_7	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del5	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del5_6	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del5_6_7	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del6	7	1483	6255	0	0	0.004720162	inf
2	ENST00000331340.8_del6_7	0	1483	6255	0	0	0	NaN
2	ENST00000331340.8_del7	0	1483	6255	0	0	0	NaN

[lonsbio]

The FASTA library we use is quite bit specific about the header format. I am looking to relax it a bit. The helper script at https://github.com/Oshlack/Toblerone/blob/master/scripts/create_index.bash is designed to help massage the input headers into the correct format for now.

I realise it may not be documented well. Can you try with that script? Or at least look at the lines

sed -e '/^>.*/s/|/~/g' -e "/^>.*/s/$/|${gene_name}|3|4|5|6|7|8|9/g" ${gene_name}_toblerone_transcriptome.fasta > ${gene_name}_toblerone_transcriptome_mod.fasta
sed -e 's/::chr.*:[0-9]*-[0-9]*|/|/g' ${gene_name}_toblerone_transcriptome_mod.fasta > ${gene_name}_toblerone_transcriptome_input.fasta  

[guillaumecharbonnier]

Thanks for your answer, I came up with something similar:
https://github.com/guillaumecharbonnier/mw-lib/blob/0d7bfa84de0f0caa2dd86324f20b0b21a6377750/src/snakemake/rules/toblerone.smk#L89

There is a typo in your shebang:

Toblerone/scripts/create_index.bash

Line 1 in 971ef2a

#!/bin/env hash

I also wrote a script to retrieve bed12 from UCSC from gene symbol that may be of interest for you:
https://github.com/guillaumecharbonnier/mw-lib/blob/0d7bfa84de0f0caa2dd86324f20b0b21a6377750/src/python/retrieve_bed12.py