split_cesar_jobs error

[KabitaBaral1]

Hi, I am running into error when I run TOGA at split_cesar_jobs step:

My code is:
./toga.py /work/dk_lab/Cactus-447way/paf2chain/Lemur_catta.chain /work/dk_lab/Cactus-447way/Input_chain/toga.transcripts.bed /work/dk_lab/Cactus-447way/All_fasta/new_Twobit/Homo_sapiens.2bit /work/dk_lab/Cactus-447way/All_fasta/new_Twobit/Lemur_catta.2bit --kt --pn TOGA_Lemur_minimap -i /work/dk_lab/Cactus-447way/Input_chain/toga.isoforms.tsv --nc /work/dk_lab/Cactus-447way/TOGA/TOGA/nextflow_config_files/ --cb 16,32 --cjn 300 --msc 100 --ncf

I tried to edit the number of buckets into --buckets 8,16,32,64,128,256,512,999 based on suggestion on issue #121

My error message stays the same regardless and looks like this:

• selected chain class to annotate transcript ENST00000491939.TXNRD2: ORTH
• selected chain class to annotate transcript ENST00000370655.ANKRD2: ORTH
  split_cesar_jobs: number of transcripts to create CESAR jobs: 37734
  split_cesar_jobs: total number of 104782 transcript/chain pairs
  split_cesar_jobs: skipped total of 3 transcripts
  split_cesar_jobs: out of them, transcripts not intersected by chains: 3
  split_cesar_jobs: assigning MISSING class to 3 transcripts not intersected by any chain
  split_cesar_jobs: creating a list of RAM-limit buckets based on user arguments
  split_cesar_jobs: defined memory limit: 999, RAM-limit buckets: {8: [], 16: [], 32: [], 64: [], 128: [], 256: [], 512: [], 999: []} (to be filled with CESAR jobs)
  split_cesar_jobs: reading bed file /work/dk_lab/Cactus-447way/TOGA/TOGA/TOGA_Lemur_minimap/temp/toga_filt_ref_annot.bed
  split_cesar_jobs: got data for 39664 transcripts
  split_cesar_jobs: reading transcript fragments data from /work/dk_lab/Cactus-447way/TOGA/TOGA/TOGA_Lemur_minimap/temp/gene_fragments.txt
  split_cesar_jobs: got data for 22965 transcripts potentially fragmented in the query genome
  split_cesar_jobs: precomputing query regions for each transcript/chain pair
  split_cesar_jobs: batch size: 37734
  split_cesar_jobs: first, invert gene-to-chains dict to chain-to-genes
  split_cesar_jobs: for each of 42342 involved chains, precompute regions
  Error! Requested range num 0 lies in the chromosome chr1 meanwhile the chain covers chrom ????U in target genome
  Traceback (most recent call last):
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 1054, in
  main()
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 873, in main
  regions, skipped_2, predef_glp = precompute_regions(
  ^^^^^^^^^^^^^^^^^^^
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 505, in precompute_regions
  chain_coords_conv_out.append(raw_ch_conv_out[i].decode("utf-8"))
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  UnicodeDecodeError: 'utf-8' codec can't decode byte 0xa0 in position 6: invalid start byte

Is there anything you'd suggest to fix the issue?
Thank you

[MichaelHiller]

My guess is that the "chrom ????U" is where the error comes from. Likely a non-standard (binary??) character.

Target likely refers to hg38. Could you pls check if
twoBitToFa Homo_sapiens.2bit test.fa
faToTwo test.fa Homo_sapiens.new.2bit
gives the same file size?
And if twoBitInfo on both 2bit files gives identical files (diff them).

And lets look if there is a non-standard name in there. Might be hard to see in a standard editor though.
Google found this

We can use this command to find all non-ASCII characters:

$ grep --color='auto' -P -n "[\x80-\xFF]" sample.txt

[KabitaBaral1]

Hi Michael,

Thank you for the prompt response. Yes, I have confirmed that the file size and contents are identical.

non-standard names: - I ran the script on .fa file, chain and bed files and found no matches.
grep --color='auto' -P -n "[\x80-\xFF]" test.fa

Any other suggestions? I thought the issue might be from the gene_fragements.txt file because its when TOGA gets to that step that I run into the error. I am not sure how to fix it though.

[MichaelHiller]

Hmm, I'll ask Yury but I think we have never seen an error like this, hence I assume that there is a problem in the input somewhere.

Could you pls also check with grep the input transcripts (bed and isoform file) and the chain file?
There shouldn't be any non ASCI as they are all text files.

[KabitaBaral1]

Thank you Michael, I have confirmed that with all the files and there's no non ASCII.

I was able to fix that error by replacing 'utf-8' with 'latin-1' on line 505 of ./split_exon_realign_jobs.py file
But now I get this error:

split_cesar_jobs: got data for 22965 transcripts potentially fragmented in the query genome
split_cesar_jobs: precomputing query regions for each transcript/chain pair
split_cesar_jobs: batch size: 37734
split_cesar_jobs: first, invert gene-to-chains dict to chain-to-genes
split_cesar_jobs: for each of 42342 involved chains, precompute regions
Error! Requested range num 0 lies in the chromosome chr1 meanwhile the chain covers chrom in target genome

• !!transcript ENST00000371370.LRRC42 / chain 146007 pair leads to very short query locus: 48bp whereas limit is 64.35000000000001. CESAR will not be called for this pair. Defining a class automatically.
  Traceback (most recent call last):
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 1054, in
  main()
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 873, in main
  regions, skipped_2, predef_glp = precompute_regions(
  ^^^^^^^^^^^^^^^^^^^
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 559, in precompute_regions
  proj_stat = define_short_q_proj_stat(q_chrom, q_start, q_end, q_2bit)
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 397, in define_short_q_proj_stat
  query_chrom = query_genome_sequence[q_chrom]
  ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^
  KeyError: ''

[MichaelHiller]

Great that the first problem is fixed. I am not familiar with 'utf-8' vs. 'latin-1' which is maybe a system setting thing.

"Requested range num 0 lies in the chromosome chr1 meanwhile the chain covers chrom in target genome"
comes from

modules/chain_coords_converter_slib.c:                    fprintf(stderr, "Error! Requested range num %d lies in the chromosome %s "
modules/chain_coords_converter_slib.c-                                    "meanwhile the chain covers chrom %s in target genome\n",
modules/chain_coords_converter_slib.c-                            i, t_regions[i].chrom, head.tName);

meaning 'chrom' is something that occurs somewhere as a chrom-name in your input but 'chrom' doesn't match a proper chromosome. Could you pls check?

[ReverendCasy]

Hello @KabitaBaral1 ,
From what I see, chain_coords_converter C function returns a byte sequence that Python cannot further decode properly. Could you please specify the following:

1. which OS you are working on;
2. how you obtained your 2bit and chain files and whether you have anyhow modified them;
3. how chromosome names for your reference and query genomes look like?

For now, I can suggest trying specifying other encodings, such as windows-1252, at line 505.
Hope that helps.

Best,
Yury

[KabitaBaral1]

Hi both,

@ReverendCasy ,

1. I am running everything on server - Linux & below are the details:
   cat /etc/os-release
   NAME="Rocky Linux"
   VERSION="8.10 (Green Obsidian)"
   ID="rocky"
   ID_LIKE="rhel centos fedora"
   VERSION_ID="8.10"
   PLATFORM_ID="platform:el8"
   PRETTY_NAME="Rocky Linux 8.10 (Green Obsidian)"
   ANSI_COLOR="0;32"
   LOGO="fedora-logo-icon"
   CPE_NAME="cpe:/o:rocky:rocky:8:GA"
   HOME_URL="https://rockylinux.org/"
   BUG_REPORT_URL="https://bugs.rockylinux.org/"
   SUPPORT_END="2029-05-31"
   ROCKY_SUPPORT_PRODUCT="Rocky-Linux-8"
   ROCKY_SUPPORT_PRODUCT_VERSION="8.10"
   REDHAT_SUPPORT_PRODUCT="Rocky Linux"
   REDHAT_SUPPORT_PRODUCT_VERSION="8.10"

2. 2bit - I converted fa to 2bit using fatoTwoBit, where my input fasta files were obtained from Cactus 447 mammalian genome dataset.
   Chain file - using minimap 2 to convert fa to paf with the following command
   minimap2 -c /work/dk_lab/Cactus-447way/All_fasta/new_fa/Homo_sapiens.fa /work/dk_lab/Cactus-447way/All_fasta/new_fa/Lemur_catta.fa -t 3 > Lemur_catta.paf
   then, paf to chain using paf2chain

3. chromosome names: for query and target in chain and 2 bit files and .bed file.

chromosome_names.txt

Michael,
I don't see "chrom" anywhere in my input files. I have checked all of them.

Thank you.

[MichaelHiller]

Sorry, 'chrom' was wrong.

covers chrom %s in target 

indicates that the string variable %s is empty and that something in your input is incomplete.
Any empty lines?

I mentioned this before, but I am not sure if extracting an entire genome from a Cactus alignment would work. The human genome should be available for download.

[KabitaBaral1]

Hi Michael, I downloaded human 2bit file and .bed file from TOGA's data available. I also compared the files generated and my .2bit and chain files to confirm if any columns are missing. The only thing that is different in my chain file is that scaffold naming for Lemur. My data is from cactus so it uses UCSC nomenclature for scaffold. I am unsure how that might be an issue.
@ReverendCasy is there anything else you'd recommend.

If it helps, this is how the error looks like:
if you check the lines that says Error!, you can see that for chr11 it says range num 0, 1 ... I have noticed that the chromosome name is different everytime I run it but for each chromosome you see numbers ranging as 0,1,2,3,4... upto 6 in some instances.

I can't tell what information is missing though and feel like I have hit the dead ends.
I have attached the input files I used if you could please run it and see if it works on your end.

Much appreciated. Thank you.

• selected chain class to annotate transcript ENST00000621453.FMR1: ORTH
• selected chain class to annotate transcript ENST00000619416.KIAA0586: ORTH
  split_cesar_jobs: number of transcripts to create CESAR jobs: 37734
  split_cesar_jobs: total number of 104782 transcript/chain pairs
  split_cesar_jobs: skipped total of 3 transcripts
  split_cesar_jobs: out of them, transcripts not intersected by chains: 3
  split_cesar_jobs: assigning MISSING class to 3 transcripts not intersected by any chain
  split_cesar_jobs: creating a list of RAM-limit buckets based on user arguments
  split_cesar_jobs: defined memory limit: 32, RAM-limit buckets: {16: [], 32: []} (to be filled with CESAR jobs)
  split_cesar_jobs: reading bed file /work/dk_lab/Cactus-447way/TOGA/TOGA/TOGA_Lemur_minimap/temp/toga_filt_ref_annot.bed
  split_cesar_jobs: got data for 39664 transcripts
  split_cesar_jobs: reading transcript fragments data from /work/dk_lab/Cactus-447way/TOGA/TOGA/TOGA_Lemur_minimap/temp/gene_fragments.txt
  split_cesar_jobs: got data for 22965 transcripts potentially fragmented in the query genome
  split_cesar_jobs: precomputing query regions for each transcript/chain pair
  split_cesar_jobs: batch size: 37734
  split_cesar_jobs: first, invert gene-to-chains dict to chain-to-genes
  split_cesar_jobs: for each of 42342 involved chains, precompute regions
  Error! Requested range num 0 lies in the chromosome chr11 meanwhile the chain covers chrom in target genome
  Error! Requested range num 1 lies in the chromosome chr11 meanwhile the chain covers chrom in target genome
• !!transcript ENST00000614710.MOB2 / chain 104059 pair leads to very short query locus: 20bp whereas limit is 35.25. CESAR will not be called for this pair. Defining a class automatically.
  Traceback (most recent call last):
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 1054, in
  main()
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 873, in main
  regions, skipped_2, predef_glp = precompute_regions(
  ^^^^^^^^^^^^^^^^^^^
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 559, in precompute_regions
  proj_stat = define_short_q_proj_stat(q_chrom, q_start, q_end, q_2bit)
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/work/dk_lab/Cactus-447way/TOGA/TOGA/./split_exon_realign_jobs.py", line 397, in define_short_q_proj_stat
  query_chrom = query_genome_sequence[q_chrom]
  ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^
  KeyError: ''
  Could not split CESAR jobs!

script:
./toga.py /work/dk_lab/Cactus-447way/paf2chain/Lemur_catta.chain /work/dk_lab/Cactus-447way/Input_chain/toga.transcripts.bed /work/dk_lab/Cactus-447way/All_fasta/new_Twobit/Homo_sapiens.2bit /work/dk_lab/Cactus-447way/All_fasta/new_Twobit/Lemur_catta.2bit --kt --pn TOGA_Lemur_minimap -i /work/dk_lab/Cactus-447way/Input_chain/toga.isoforms.tsv --nc /work/dk_lab/Cactus-447way/TOGA/TOGA/nextflow_config_files/ --cb 16,32 --cjn 300 --msc 100 --ncf

The only modifications I made was change "utf-8" to "latin-1" on line 505
and int() to int(float()) on line 518 on split_exon_realign_jobs.py file.

link to input files:
https://drive.google.com/drive/folders/1p6CYw-LfpiczgNGLS13AZeKIPBLzNYXA?usp=sharing

[MichaelHiller]

the chrom names match between the 2bits and the chain.
Lemur has a 'scaffold' while all others are 'scaffold41' or 'scaffold_23'. But 'scaffold' without a number is in the 2bit.

I only noticed that 3 chains are not tab separated (unlikely to a be problem).

chain 255 GL000009 201709 + 153693 158982 scaffold 285822530 + 27375281 27380955 13558
chain 255 GL000205 185591 + 156642 157961 scaffold1 83724376 - 76362388 76363735 163844
chain 255 GL000216 176608 + 109150 109768 scaffold3 68826487 - 4975212 4975840 181474

But here is a different problem.
You have
g chain Lemur_catta.chain -c
286397
chains in your file.

Many of these chains are very short, which is either because your assembly is fragmented (GCF_020740605.2) or the Cactus chains get fragmented during the hal conversion (or however you got them). Also, they lack the proper scores, as all your scores are 255. TOGA filters chains less than 1000 by default.

Here is a browser shot to show the level of fragmentation. I would recommend using make_lastz_chains to compute the chains.

[KabitaBaral1]

I see, thank you Michael. I looked into the scores and apparently the UCSC chain scores have highest possible score value of 255 so I added two flags --msc 100 --ncf.
Is the issue because of how highly fragmented the chain is?

I also downloaded the Lemur.fa file from NCBI (the same one you used for TOGA) and that gave me the same issue as this one. So, I am not sure if the issue is because of how fragmented my genome is.

On a different note,
I tried to use make_lastz_chain but it ran for over 2 days and didn't finish running.

Here's the setting I used:
From my parametrs.json file:
"target_name": "Homo_sapiens",
"query_name": "Acinonyx_jubatus",
"project_dir": "test_out",
"continue_from_step": null,
"cluster_executor": "slurm",
"cluster_queue": "apophis",
"cluster_cpus": 1,
"cluster_time": "48:00:00",
"cluster_memory": "2G",
"params_from_file": "my_awesome_config.json",
"seq_1_dir": "/work/dk_lab/Cactus-447way/LASTZ/make_lastz_chains/Cactus_out/target.2bit",
"seq_1_len": "/work/dk_lab/Cactus-447way/LASTZ/make_lastz_chains/Cactus_out/target.chrom.sizes",
"seq_2_dir": "/work/dk_lab/Cactus-447way/LASTZ/make_lastz_chains/Cactus_out/query.2bit",
"seq_2_len": "/work/dk_lab/Cactus-447way/LASTZ/make_lastz_chains/Cactus_out/query.chrom.sizes",
"lastz_y": 9400,
"lastz_h": 2000,
"lastz_l": 3000,
"lastz_k": 2400,
"seq_1_chunk": 175000000,
"seq_1_lap": 0,
"seq_1_limit": 5000,
"seq_2_chunk": 50000000,
"seq_2_lap": 10000,
"seq_2_limit": 10000,
"chain_min_score": 1000,
"chain_linear_gap": "loose",
"fill_chain": true,
"fill_unmask": true,
"fill_chain_min_score": 25000,
"fill_insert_chain_min_score": 5000,
"fill_gap_max_size_t": 20000,
"fill_gap_max_size_q": 20000,
"fill_gap_min_size_t": 30,
"fill_gap_min_size_q": 30,
"fill_lastz_k": 2000,
"fill_lastz_l": 3000,
"fill_memory": 16,
"fill_prepare_memory": 50,
"num_fill_jobs": 200,
"clean_chain": true,
"chaining_memory": 16,
"chain_clean_memory": 100,
"clean_chain_parameters": "-LRfoldThreshold=2.5 -doPairs -LRfoldThresholdPairs=10 -maxPairDistance=10000 -maxSuspectScore=100000 -minBrokenChainScore=75000",
"job_time_req": "48h",
"keep_temp": false

and the script:
./make_chains.py Homo_sapiens Acinonyx_jubatus /work/dk_lab/Cactus-447way/Input_chain/Human.2bit /work/dk_lab/Cactus-447way/All_fasta/new_Twobit/Acinonyx_jubatus.2bit --pd Cactus_out --force --cfs lastz --chaining_memory 150 --fill_memory 150 --cluster_executor slurm --cluster_queue apophis --params_from_file my_awesome_config.json --job_time_req 48:00:00

For the chain files you generated for TOGA, what were the specifications you used? Do you mind sharing the parameters.json, lastz_config.nf files and the script you used please? Also, how long did each genome take to run for TOGA's chain files?
I ran everything on slurm, but as it is, it is not feasible to use make_lastz_run because of how long it took to run it on one genome.

Thank you.

[MichaelHiller]

Even the very small chains will have scores of 5000 and higher. In principle, considering all chains with a score of 100 should work.

Close species alignments can take some time. Did it crash? Did you repeatModel both the ref and the query and mask with the lib? If the masking comes from NCBI, then I can guarantee that this won't work.

What is the Lemur GCA/GCF ID? Maybe I have the chains already.

The default parameters in make_chains is what we also use.

[KabitaBaral1]

what do you mean by repeatModel ref and query and mask with the lib? I tried looking it up but all I can find is that I might need to use RepeatModeler and then run RepeatMasker for output of RepeatModeler. Is this what you are refering to?
If yes, do I run it on .fa files?

The lemur was just a test run. I need .chain files for all 447 mammals from Cactus 447 mammalian genomes. I compared the list to the ones you have from TOGA and here's the list of the genomes that I need to generate .chain files. It didn't crash but even for cheetah (which was a test run, it ran for over 2 days and took a long time to finish). Do you recommend changing anything from the settings I have or should I expect it to take a few days for primate genomes?

PS: If you have any chain files for the genomes listed in the .txt file attached, could you share them please? Thank you, Michael.
unique_cactus.txt

[MichaelHiller]

"use RepeatModeler and then run RepeatMasker for output of RepeatModeler." --> Yes that is what I meant

2 days is relative, as this depends on how many jobs are running in parallel.

We will start TOGA2 runs for a lot of mammals very soon. If you send me a list of GCA/GCF accession numbers, I can check which ones we don't yet have aligned.

[KabitaBaral1]

I see, Thank you Michael. I have created a spreadsheet with all the genomes and highlighted the ones for which I couldn't find chain files in your data from TOGA paper (https://genome.senckenberg.de/download/TOGA/human_hg38_reference/hg38_chains/)
If you have chain files for them, that'd be great. Thank you again :) I appreciate it.

Genomes_Chain_KB.csv