Questions about using GetOrganelle to assemble/bait Angiosperms353 loci sequences

[CharlesLanceolata]

First of all, thank you very much for this really convenient tool! I have successfully assembled plastome sequences and nrDNA cistron sequences for my research. :D

I found out in Chen et al. (2023) that GetOrganelle could be used to assemble Angiosperms353 loci sequences, and the assembled sequences yield better downstream analysis results than Easy353. However, since assembling one loci at once is too time & resource consuming (even if GNU parallel is used), I tried to generate my own label database and am able to assemble all loci sequences in the same round.

But now I face another issue: the sequences in *.scaffolds.graph1.1.path_sequences.fasta do not have headers containing loci names, and manually BLASTing these sequences or identifying them via Bandage is laborious as well. Therefore, can you kindly suggest a convenient/lazy way to do this? I'm really sorry if I overlooked something in the Wiki and Discussion, leading to me asking these ignorant questions.

Anyways, thank you for your time and patience regarding to this matter. Looking forward to hearing from you! ^^

[JianjunJin]

I am pleased to learn that GetOrganelle has benefited your research, and I appreciate you bringing this issue to our attention.

While there have been inquiries about assembling custom loci using GetOrganelle, and I am aware of the Corydalis study you mentioned, I must admit that I had not previously delved into their methodology. It is indeed intriguing to see that their approach results in improved downstream analysis compared to Easy353, and their rationale appears sound. However, I have not yet incorporated 353 data into my own research.

Regarding the current version of GetOrganelle, it is not specifically optimized for this particular purpose, which means there is no straightforward one-command solution available. Nonetheless, since you have successfully created a customized label database and utilized it for assembly, the subsequent steps should be more manageable. Here is what you should do next: Examine the accompanying *.csv file (which is actually in TAB format) produced by GetOrganelle. This file contains all the gene marker labels for the target contigs. In your downstream process, you will need to create and execute a custom script that 1) parses the *.csv TAB file and creates a map of contig-to-gene associations, and 2) uses this map to extract the relevant contigs from the current FASTA output, storing them as individual gene.fasta files.

Should you encounter any difficulties, please do not hesitate to contact me. If you need assistance with scripting, feel free to send me one of your current outputs (*.csv and *.fasta), and I may be able to help in crafting a Python script.

[CharlesLanceolata]

Wow, why didn't I thought of that, I definitely will first try writing a script that does what you've suggested, if anything goes wrong then I'll come back for help.
Marking your comment as answer for now, thank you very much Dr. Jin! 😀

[CharlesLanceolata]

Hi Dr. @JianjunJin,

I'm not sure if it is recommended to open another discussion, but since this following question is still related to "using GetOrganelle to assemble Angiosperms353 loci sequences", and in order to make others easy to locate my questions which they may also encounter, I've decided to edit the title of this discussion. Please let me know if I should split this comment apart.

After playing around with my data, I found out that there were two "edges" that has more than one BLAST hits (as seen in the attached *.csv file), although only one of the edges were selected into the *.scaffolds.graph1.1.path_sequence.fasta file, I'm not sure which treatment would be better:

1. Discard the "chimeric" scaffolds, as they may be a result of incorrect assembly, and proceed with downstream analysis.
2. Try to select another edge with only one loci hit, and has the second highest coverage.
3. Try other seed and label dataset combinations (with the remaining options unchanged, as in the log file attached), see if they do any differences.

Eventually, I've chosen to test the third option, and discovered that the number of loci sets included in the seed file (e.g., if I take Angiosperms353 loci from 2 samples as seed, then there will be 2 sets of loci) affected the resulting number of scaffolds, and also the presence of "chimeric" scaffolds in the *.scaffolds.graph1.1.path_sequence.fasta file, as listed in the following table:

Attempts	Sets of loci in the seed file	Sets of loci in the label dataset	Final no. of scaffolds	Edges with more than 1 Angiosperms353 loci blast hits
1	2	many	321	1859 (loci hits: 7141,7279) 2053 (loci hits: 6448, 6544)
2	many	many	284	n/a
3	1	many	317	n/a
4	2	2	324	1869 (loci hits: 7141,7279) 2071 (loci hits: 6448, 6544)
5	1	2	319	n/a
6	1	1	319	n/a

For the seed & label files that included "many" sets of loci, I've inputted the mega353.fasta as suggested in HybPiper (the file can be downloaded here), they are more distantly related to the taxa I'm interested in. As for those that included 1 or 2 sets of loci, I've selected closely related taxa, maybe this is the reason that made the differences? Another thing is that, I've tested the label dataset with only 1 set of loci included due to your statement in the Wiki page:

But in the FAQ, you've mentioned that the seed file and label dataset can be the same.

I'm not quite sure how to interpret this result and whether I can use these assembled sequences in further analyses, since I'm afraid that there might be some aspects that I've been missing (e.g., assembly errors, wrong associations between scaffolds and loci). Therefore, any suggestions are appreciated!

Thank you Dr. Jin, if any other files or record is needed please let me know.

K01-353extended_K77.assembly_graph.fastg.extend-label.csv
K01-353get_org.log

[JianjunJin]

Thank you for providing the details.

Firstly, I must clarify that what appear to be two hits do not represent two distinct genes. They are more likely two conserved regions within a single gene, separated by a variable region that fails to match any entries in the database through BLAST. So, I wouldn't claim these to be chimeric, given the current information. To do more, you'd better follow the next point.

Secondly, based on the .csv results, it appears that the label file you created does not adhere to the recommended format. Specifically, the sequence header format does not align with the guidelines provided here or follow the example of the default mitochondrial database. Consequently, your generated .csv file lacks gene information. By preparing the file correctly, you will be able to directly identify the contigs and their associated genes from the .csv file.

Only after understanding these points should you consider the possibility of a contig genuinely hitting multiple genes, which might occasionally occur due to the similarity between two loci. This issue can often be resolved by adjusting the BLAST e-value option. For example, you can use --slim-options '-e 1e-50' with get_organelle_from_assembly.py to refine your search (note that this option is not currently available in get_organelle_from_reads.py, but you can process the assembly based on your results, which should be relatively quick).

Let me know.