Baiting for ITS1-5.8S-ITS2 ribosomal regions

[max-baer]

Hi Kinggerm,
I am using your software to try and bait for ITS regions in WGS data of a trichuris species. We conducted WGS with approx. 6x coverage, paired end (mean read length 185 bp). I used the following command line: (I adjusted the word size to 110 given the long read lengths)

get_organelle_from_reads.py -s Trichuris_ITS_reference_all.fasta --genes trichuris_ITS.label.fasta -1 Sample1_forward_paired.fq.gz -2 Sample1_reverse_paired.fq.gz -o Sample1_ITS_trimmed -F anonym --reduce-reads-for-coverage inf --max-reads inf -R 10 -w 110 -t 4 -P 0

Some elaboration:
-s <Reference Trichuris ITS sequences as seed (See attachment)>
--genes <Label file: 2 of the previously assumed closest related species>
-F anonym (as advised for custom assamblies)

I have a bunch of questions regarding the software and I am a bit new to the field (and super interested). However, I understand that your time is probably super limited already, and am happy for any answers I can get. Thanks so much for the open access of your software and any responses you might be able to provide :):)

1. The ITS regions are not circular is that a problem?

2. Is it a problem that my target region is shorter than 1000bp?

3. I looked at the K115 output in bandage (see pictures attached) and in one image there is a super chaotic picture but you see the labels. I then manually cleaned this up. Am I correct in the assumption that this is exactly step 4 in the getorganelle software? (In the final anonym.K115.contigs.graph1.selected_graph.gfa I get one single 450 bp long sequence without any nodes, see final)

4. I had several warning messages (log file attached, most of them failing Disentangling). Should I be concerned about this or can I use the output file as is (anonym.K115.scaffolds.graph1.1.path_sequence.fasta)?

5. What metric can I use to assess the “quality”?

6. Am I correct in assuming that this software is not compatible with 3rd generation (nanopore) sequencing, because it requires forward and reversed inputs? (SPAdes dependency)

7. I constructed I phylogenetic tree (see attached) which is how I expect it to be, however I used only the closest related ITS regions in the label file. Does this bias my output towards the sequences contained in the label file?

Greetings from Basel, Switzerland,
Max

get_org.log.txt
Trichuris_ITS_reference_all.txt
trichuris_ITS_label.txt

[Kinggerm]

Hi Max,

Thanks for reaching out with detailed information.

• (1) It's expected to be either circular or linear. For more information, please see here https://github.com/Kinggerm/GetOrganelle/wiki/FAQ#why-does-getorganelle-generate-a-circular-genome-or-not-for-embplant_nrfungus_nr.
• (2&3&4). A shorter length is problematic. It results from the conservative shortcut that the current GetOrganelle version took, which only assigns the unique label to only one contig, which should work with a complete circular result but not the other way around. I have improved this in the next version (available at the intermediate_graph branch); however, it will take a while for me to test the new version comprehensively before releasing it.
  I believe you are doing step 4 manually and correctly. The filtered*fastg file can be helpful too. I usually check it rather than the original fastg file. You can also do step 4 using a semi-manual (combining "slim_graph.py") approach.
• (5) Probably using evaluate_assembly_using_mapping.py and/or building a phylogeny as you did.
• (6) Yes and no. You can input single-end short read using -u. It's not compatible mainly because using long error-prone reads alone can be problematic. I'm developing another tool (Traversome) for combining short- and long- reads to solve a more generic problem.
• (7) It depends on how diverged the sequences are. If it were plastomes within angiosperm scope, which evolved slowly, I would say no influence at all. I don't know much about animal ITS. But you can easily confirm this by using a more related label.

Extra tip, given the depths & the graph, you can speed up each run by using even a larger -w and fewer rounds.

Best,
Jianjun

[max-baer]

Hi Jianjun,

This is an awesome feedback thanks so much, really appreciate it!

Regarding (2&3&4): I see, so the software in the end "picks" one contig out of the many in the above picture which were labeled. Is this done based on an optimality metric, or out of a subset of optimally fitting contigs one is chosen at random, or the first best is chosen?

I will play around with the software a bit more and also with the current intermediate branch to see a bit what comes out.

Again thanks so much for the quick response!

All my best,
Max

[Kinggerm]

It's calculating a customized hit weight value (see more in the manuscript) as the metric, which is mainly for identifying the correct graph component in the complete circular case and works poorly for incomplete results.

The latest branch should be an improved solution. However, blast hits have many flaws too. The best solution I can envision should be conducting annotation along the assembly graph, which I haven't seen yet.