Applying singleM to single cell sequencing data

[lingrongjin]

I'm interested in applying singleM to single cell amplified genomes (SAGs) to get their taxonomic classifications. I have a few questions that I would like to hear your suggestions:

1. SAGs have very uneven read coverage distribution (some region of the genome can be unsequenced) due to the MDA process; I'm wondering if this will affect singleM's marker gene based approach for taxonomic classification?
2. Would it be possible to add custom MAGs to the existing reference database and run singleM profiling on the combined new database?

[wwood]

Hi @lingrongjin

1. The standard -p taxonomic profile output from singlem pipe uses an algorithm ('condense') which assumes equal coverage across the genome, so you might get some strange results there.

The --otu-table (possibly with --output-extras) will give you an taxonomic assignment for each gene. So that might give you an answer which is largely unaffected by the ups and downs of SAG coverage, at the expense of being a bit harder to interpret (you get one or more taxons for each marker genes, so you have to bring them together yourself by e.g. interpreting a spreadsheet by eye, or by writing some small processing script).

2. Yes, that is what singlem supplement does. Please see the help. There is also a half-supported snakemake script which helps prepare inputs for supplement when you are adding very many genomes at a time at https://github.com/wwood/singlem/tree/main/extras/prepare_for_supplement

HTH, ben

[lingrongjin]

Hi Ben,

Thanks for your reply and suggestions! The --otu-table command looks helpful. Looks like I can use it to identify single species SAGs with consistent taxonomy of the marker genes.

I tried singlem supplement but ran into some issues during the gtdbtk taxonomic classification process: extern.ExternCalledProcessError: Command gtdbtk classify_wf --cpus 32 --batchfile /tmp/tmp_29c89cr/tmp/tmpdxttuqww --out_dir /tmp/tmp_29c89cr/gtdbtk_output --mash_db /tmp/tmp_29c89cr/gtdbtk_mash.msh returned non-zero exit status 127.
STDERR was: b'bash: gtdbtk: command not found\n'STDOUT was: b''

I'm wondering if you have any suggestions on how to solve this issue (below is the code I used):
singlem supplement --new-genome-fasta-files-list n620_genome_list.txt
--input-metapackage S4.3.0.GTDB_r220.metapackage_20240523.smpkg.zb
--output-metapackage gtdb220_n620.smpkg.zb
--no-quality-filter
--no-dereplication
--no-taxon-genome-lengths
--threads 32

[wwood]

Ah. That occurs because the python version needed by gtdbtk would constrain the python version of all singlem's code, if it were included as a conda dependency. There's a few ways to solve this, but might be easiest to run gtdbtk yourself on the genomes, and then provide --gtdbtk-output-directory to supplement.

Or you could use the Snakemake, or change your PATH so that gtdbtk is on it as well.

[lingrongjin]

I provided --gtdbtk-output-directory, but it's producing this error now:
Exception: The taxonomy d__Bacteria;p__Pseudomonadota;c__Gammaproteobacteria;o__Burkholderiales;f__Burkholderiaceae_B;g__Aquabacterium for genome 24AcMG_bin.1 (originally d__Bacteria;p__Pseudomonadota;c__Gammaproteobacteria;o__Burkholderiales;f__Burkholderiaceae_B;g__Aquabacterium;s__) is not a known taxonomy in the current metapackage

Is it because my MAG is not classified at species level and hence the taxonomy is unknown? Would singleM be able to handle novel genomes like this?

[wwood]

I don't think that's why. More likely a mismatch in gtdbtk and singlem database version. Are they both using gtdb R220? -\------------- Ben Woodcroft Group leader, Centre for Microbiome Research, QUT

…

________________________________ From: Lingrong Karen Jin ***@***.***> Sent: Friday, January 17, 2025 6:57:33 PM To: wwood/singlem ***@***.***> Cc: Ben J Woodcroft ***@***.***>; Comment ***@***.***> Subject: Re: [wwood/singlem] Applying singleM to single cell sequencing data (Issue #206) I provided --gtdbtk-output-directory, but it's producing this error now: Exception: The taxonomy d__Bacteria;p__Pseudomonadota;c__Gammaproteobacteria;o__Burkholderiales;f__Burkholderiaceae_B;g__Aquabacterium for genome 24AcMG_bin.1 (originally d__Bacteria;p__Pseudomonadota;c__Gammaproteobacteria;o__Burkholderiales;f__Burkholderiaceae_B;g__Aquabacterium;s__) is not a known taxonomy in the current metapackage Is it because my MAG is not classified at species level and hence the taxonomy is unknown? Would singleM be able to handle novel genomes like this? ― Reply to this email directly, view it on GitHub<#206 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAADX5EPTLBFWUK5EQWRH732LDAX3AVCNFSM6AAAAABVIS37MKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKOJXG4ZTQMRZGM>. You are receiving this because you commented.Message ID: ***@***.***>

[lingrongjin]

Hi Ben, you are right. I changed the database version for gtdbtk and it ran without errors. I have one additional question (thanks again for your explanations and suggestions so far!) about the use of the otu-table output to infer SAG taxonomy. I assume most of my SAGs represent single species bacterial/archeal genomes, but I found that the otu tables for many SAGs contain marker genes from different taxa. I'm wondering if you think it's reasonable to classify the SAGs using some rule like if >50% of the marker genes match a specific taxonomic level (e.g. a genus), we could say the SAG belongs to the specific taxa?

[wwood]

Glad it worked out - I've modified the error message to suggest this solution automatically.

Are you sure that the SAG is uncontaminated? The procedure is imperfect in my experience, which can result in multiple cells in the same well. OTOH, sometimes there is just noise. Can you provide a concrete example either here, or if you are worried about data privacy, over email please?
Thanks.

[lingrongjin]

Hi Ben, you are right, contamination is an issue of concern. Ideally, I would like to discard SAGs with multiple cells from different species and keep those with single species. However, I ran singlem over 10000+ SAGs and about 70% of them contain markers from different taxa at genus level or above. I attached two examples of SAGs with markers from multiple taxa. It appears to me that in the first case (Ac00042), the majority of the marker genes match to a single species, whereas in the second case (Ac00260), the taxonomy of the marker genes is everywhere, which may indicate contamination. I'm wondering if you have any suggestions on distinguishing contamination (multiple cells) vs noise in the marker gene detection process? Thanks!

sag_w_multiple_markers.csv

[wwood]

Hi again,

Thanks for the concrete data.

I agree the first is likely OK, but the second it is hard to tell I think, because (1) the isn't many markers, but also (2) none of the markers are assigned to the species level. There's quite a difference between species level and genus/above - species level assignments mean the sequence of the OTU is 96.7%+ nucleotide ID to the reference, whereas for genus level and above it is based on best DIAMOND blastx hit(s). So for individual genes from novel lineages, the individual OTU taxonomy strings can be a bit all over the place - the condense procedure which converts these tables into taxonomy profiles applies some methods to try to account for this.

One of the things is that it requires 3 or more genes from a taxonomy in order to be included in the final output. So for that second SAG, only c__Bacteroidia and d__Bacteria will have non-zero coverage. I'm not convinced there is enough evidence to conclude that the second is contaminated.

GIven unlimited time, perhaps the best solution would be to develop a singlem metapackage with an expanded set of markers (in the order of hundreds/thousands rather than the 35 Bacterial and 37 Archaeal ones in the current metapackage). We don't have plans to make this kind of "sensitive" metapackage, but would be supportive if you wanted to give it a go.

Another option might just be to assemble/bin them, and then only use SAGs that make it through that process? Maybe I'm getting lost trying to understand what you've done so far.

Thanks.
ben

[lingrongjin]

Hi Ben,

Thanks for the detailed explanations! It makes much more sense to me now and based on my understanding (please correct me if I'm wrong):

1. if a marker gene (or OTU?) is classified at the species level, the taxonomy of the OTU is more reliable compared to when an OTU is only assigned at genus level or above based on DIAMOND blastx hits.
2. The fact that OTU taxonomies are all over the place may not necessarily be an indication of contamination, but can also be due to the fact that the SAG comes from a novel lineage not covered by the marker gene database.

The idea of an expanded set of markers sounds nice, but I'm a bit rushed on time so I'm wondering if there is a way to get reliable taxonomic classification for SAGs based on presence of marker genes. For example, if >50% of the marker genes match to a specific species/genus, we could say the SAG belongs to the species/genus (or following the logic of 'condense' algorithm, maybe I should also consider 'number of hits' and/or minimum number of marker genes to be present?)

For your suggestions, when you say "use SAGs that make it through the process", do you mean to assemble the SAGs first then classify them using some tool like gtdb-tk, or do you mean instead of using raw reads as input, we could use assembled contigs from each SAG as input to singlem pipe? For the former, I believe many of the SAGs have very low completeness (e.g., < 10%) so it may be difficult to get taxonomic classification using gtdb-tk, which is why I'm looking for other methods to get their taxonomies.