mKmer: An unbiased K-mer embedding of microbiomic single-microbe RNA sequencing data

mKmer: An unbiased K-mer embedding of microbiomic single-microbe RNA sequencing data | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Method Article mKmer: An unbiased K -mer embedding of microbiomic single-microbe RNA sequencing data Fangyu Mo, Qinghong Qian, Xiaolin Lu, Dihuai Zheng, Wenjie Cai, and 10 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-5748035/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract The advanced single-microbe RNA sequencing (smRNA-seq) technique addresses the pressing need to understand the complexity and diversity of microbial communities, as well as the distinct microbial states defined by different gene expression profiles. Current analyses of smRNA-seq data heavily rely on the integrity of reference genomes within the queried microbiota. However, establishing a comprehensive collection of microbial reference genomes or gene sets remains a significant challenge for most real-world microbial ecosystems. Here, we developed an unbiased embedding algorithm utilizing K -mer signatures, named mKmer, which bypasses gene or genome alignment to enable species identification for individual microbes and downstream functional enrichment analysis. By substituting gene features in the canonical cell-by-gene matrix with highly conserved K -mers, we demonstrate that mKmer outperforms gene-based methods in clustering and motif inference tasks using benchmark datasets from crop soil and human gut microbiomes. Our method provides a reference genome-free analytical framework for advancing smRNA-seq studies. Bioinformatics K-mer msm RNA-seq HCK reference genome-free K-motif Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Background Recently, we developed a high-throughput single-microbe RNA sequencing (smRNA-seq) technique for microbiome samples [ 1 ], which can generate RNAs from over 5,000 single microbes of a microbial community. This technique can effectively solve the problem of significant cell heterogeneity among bacterial populations, and thus achieve complete functional characterization of host-related microorganisms. However, there is a key challenge in the microbiomics smRNA-seq (msmRNA-seq) data analysis: constructing a high-quality gene expression matrix for downstream analysis [ 2 ]. Typically, the construction of a gene expression matrix requires a reference genome [ 3 ]. While the human reference genome is largely complete and accurate, this is not the case for other organisms. In particular, many species in a microbiome sample lack reference genomes or do not have a high-quality annotated gene set [ 4 ]. These issues are especially pronounced when sequencing data encompasses multiple unclassified species (e.g., microbiome data). Any attempt to predefine which species should be included in the reference genome(s) inevitably introduces bias. K -mer refers to short sequences of a specific length ( K ) and any genomic or RNA sequences are composed of different K -mers. K -mers have been widely used in bioinformatics analysis, including genome survey and assembly, and also single-cell omics data [e.g. 5, 6, 7]. Theoretically, RNAs from a microbiome sample can be characterized by a specific length of short sequence (i.e., K -mers). In this study, we developed a new frame of msmRNA-seq analysis (named mKmer) based on high-frequency conserved K -mer (HCK) rather than a gene expression matrix. Benchmark tests on seven datasets from soil and gut demonstrate that mKmer significantly improved species identification compared to the cell-by-gene matrix. To demonstrate applications of our method, we used a clinical msmRNA-seq data from the gut microbiome of colorectal cancer patients before and after treatment. Results Overview of mKmer method mKmer is a tool for msmRNA data analysis by constructing a cell-by-HCK matrix to achieve efficient cell classification, species annotation, and functional analysis. The tool is currently divided into seven analysis modules, including KmerRank and KmerCell , aiming to provide personalized services for analysts (Fig. 1 ). mKmer extracts biological information from raw sequencing data by testing different K -mer lengths to determine the optimal K (Fig. 2 ). This optimal K can identify key conserved sequences while effectively distinguishing noise from non-critical gene sequences. In our analysis framework, we successfully identified species and performed functional analysis for soybean soil samples and human fecal samples by selecting these HCKs before the inflection point. This demonstrates that the HCKs already densely contain a large amount of hierarchical information in microbial taxonomy, enabling precise and efficient classification at various levels, including domain, phylum, class, order, family, genus, and species. On the other hand, the low-frequency K -mers contain sparse taxonomic information, and discarding them does not affect the classification results. Therefore, the HCKs means that may be translated into amino acid sequences (i.e. protein motifs) accordingly. Additionally, HCKs address issues of computational time and high matrix dimensionality caused by the excessive variety of K -mers, which significantly enhanced mKmer's performance and practicality. Because it does not require any reference genome, the reproducibility of sequencing reads reached 100%. K-mer scanning and rank plot For a raw msmRNA sequencing dataset, cells were selected as usual (such as with UMI-tools), and reads from the selected cells were filtered to remove duplicates before the downstream K -mers scanning (Fig. 1 ). We tested frequency distributions of different K -mers sizes for the seven msmRNA datasets and found a change of distribution curves happening between 12-mer and 13-mer in all samples (Fig. 2 , other five samples see Fig. S1). Due to the high-throughput smRNA-seq technique using random primers to capture RNA of individual cells, the combined amplification and release yield an average RNA sequencing depth of 1× coverage. Therefore, low-frequency K -mers with a frequency of 1, where the K value is at its peak, are considered optimal. The 12-mer was therefore used as the default size for K -mers scanning. We further ranked all 12-mers or 13-mers scanned from the msmRNA data by count per cell (Fig. 3 , other five samples see Fig. S2). The K -mer rank plot (from highest to lowest K -mer depth) is an interactive plot that shows all K -mers detected in a microbiome sample or a msmRNA-seq dataset. Identification of HCKs The overall shape of the K -mer rank plot (Fig. 3 , left) is similar to the barcode rank plot (Fig. 3 , right). Typically, a “cliff-and-knee” shape can be observed in the K -mer rank plot of a microbiome sample. In this case, the steep cliff, followed by the plateaued knee, demonstrates that the K -mer calling algorithm was able to distinguish feature K -mers from others. HCKs mainly come from the evolutionary conserved regions (e.g., motifs in protein domains and DNA-binding sites) of bacterial species in a microbiome sample. As an example, the region of the genus Staphylococcus HSP60 gene encodes the conserved NdhRMIQE motif [ 8 ], and a high number of K -mers could be counted within this region (Fig. 4 ). The conserved K -mers exist in a wide variety of microbe species in a microbiome. From a microbial taxonomy perspective, HCK corresponds to the lowest common ancestor (LCA) sequence at the taxonomic level. The use of K -mers for microbial classification has been demonstrated by many classical K -mer-based metagenomic taxonomy annotation software, such as Kraken [ 9 ], Centrifuge [ 10 ], and Kaiju [ 11 ]. HCK was first discovered and successfully applied in the identification of taxa in single-cell data. Identificaiton of maker K-mers Based on the cell-by-HCK matrix and routine cell clustering and dimension reduction (as shown in Fig. 1 ), a visualization result by uniform manifold approximation and projection (UMAP) of a msmRNA sample and species annotation can be obtained (examples shown in Fig. 5 , right, and the cell-by-gene matrix clustering results of this sample are shown in Fig. 5 , left). A good clustering of the same species/cells was observed in the UMAP plot. Further, marker K -mers can be identified among the different clusters (species or subspecies) using routine approaches (same as those for marker genes) such as the Seurat function ( FindAllMarkers ). Function annotation by K-motifs The marker K -mer can be used for functional annotation based on their encoding motifs as mentioned above. Based on motif and domain databases (MEME and Pfam), the DNA motifs and protein motifs in domains (termed the K -mer-contained motifs as K -motifs) can be identified for gene ontology (GO) annotation by mKmer functions ( KmerGOn and KmerGOp ), respectively. At the protein level, the longest translated amino acids (AA) (4-mer AA) for marker K -mers (12-mer nt) in a microbiome sample should have significant sequence similarity to the protein domain’s motifs in Pfam. The motif-contained K -mers are those highly conserved K -mers which transcript from the motif-contained regions of orthologous genes of different microbe species in a microbiome sample. Using the K -motifs and their GO IDs identified, function analysis such as GO and pathway enrichment can be done as usual (Fig. 6 D). Benchmark test with cell-by-gene matrix To compare the performance of mKmer with the traditional gene matrix-based method, we generated a msmRNA-seq dataset from soybean ( Glycine max ) soil and collected four publicly available msmRNA-seq datasets from human guts. Firstly, when comparing the dimensionality reduction and clustering results using the traditional cell-by-gene matrix (two examples are shown in Fig. 5 , left), regardless of whether the msmRNA-seq data came from soil (Fig. 5 A) or human gut (Fig. 5 B), the cell-by-HCK matrix (Fig. 5 , right) was more distinct and accurate (the other three examples of human gut shown in Fig. S4.). Secondly, for the same msmRNA-seq dataset, the number of species identified using the cell-by-HCK matrix was significantly higher than that identified using the cell-by-gene matrix, increasing by more than 4.5 times. At the same time, the number of each species in both the soil and the gut also increased significantly. This is particularly evident for five species with low or moderate abundance (number < 300), including Bordetella pertussis , Flavobacterium sp. CJ75 , and Labrys sp. KNU-23 in the soybean soil sample. The newly identified strains accounted for 5/21 of the original strains. In the human gut samples, nine types of microorganisms were newly identified, including Agathobacter rectalis , Faecalibacterium prausnitzii , Parabacteroides distasonis , Phascolarctobacterium faecium, and Sutterella wadsworthensis . The number of newly identified bacterial species in the human gut (9 species) even exceeds the original number of species (7 species). Therefore, mKmer shows significant advantages in species identification in both complex soil environments and gut environments. This advantage is particularly pronounced in datasets where the original results were not very good and the species diversity was relatively low. Numerous studies have shown that the five bacterial species only identified by mKmer are commonly found in soil. Sphingopyxis terrae and F. sp. CJ75 are capable of degrading complex organic compounds [ 12 , 13 ]. Mucilaginibacter mallensis can produce mucilaginous polysaccharides, thereby improving soil structure and fertility [ 14 ]. L. sp. KNU-23 is widely present in organic matter-rich soils [ 15 ]. Surprisingly, B. pertussis , primarily known as a human pathogen transmitted through the air, is not commonly found in soil environments, and its survival in soil has been seldom studied ( 10.1038/nrmicro886 ). The microbiome of the human gut is being studied more thoroughly. The species identified by mKmer, such as A. rectalis , F. prausnitzii , Mediterraneibacter gnavus , Odoribacter splanchnicus , P. distasonis , P. faecium , Phocaeicola coprophilus , Roseburia intestinalis , and S. wadsworthensis , are common human gut bacteria based on the literature [ 16 – 24 ]. To further validate the reliability of our results, we applied the mKmer functions KmerGOn and KmerGOp to in-depth explore the functions of the mKmer identified species. The B. pertussis in soybean soil exhibits a unique ability to bind iron ions in soybean soil (Fig. S3A, left). Iron is an essential micronutrient for plant growth, affecting soybean health and yield. Microorganisms can inhibit the growth of pathogens by competitively adsorbing iron in the soil, thereby reducing the occurrence of diseases. Compared to other microorganisms, acid phosphatase activity is significantly enriched in Bordetella pertussis (Fig. S3A, right), which can help decompose organic phosphorus compounds in the soil, releasing inorganic phosphorus that plants can absorb, thereby promoting phosphorus uptake by soybeans and enhancing crop growth. For human gut msmRNA-seq data, the A. rectalis identified by mKmer revealed processes related to the metabolism of acids, including aconitate hydratase activity and the dicarboxylic acid metabolic process using KmerGOn (Fig. S3B, left), consistent with findings by Abdugheni et al. [ 16 ]. In addition, several processes related to the biosynthesis of amines have been found by KmerGOp , such as 6-pyruvoyltetrahydropterin and tetrahydrobiopterin biosynthesis (Fig. S3B, right). Taken together, mKmer doesn’t depend on reference genomes, and can unbiasedly and effectively parse the complex biological information in msmRNA-seq data. A case study using mKmer To demonstrate the practical applications of mKmer, we collected fecal samples from a colorectal cancer patient before and after immunotherapy for single-microbe sequencing. Using mKmer to analyze the two msmRNA-seq datasets, we identified 19 microbial species in the pre-treatment fecal sample (Fig. 6 A) and 27 species in the post-treatment sample (Fig. 6 B, and the cell-by-gene matrix clustering results of these sample data are shown in Fig. S5). The microbial richness increased by over one-third, which is conducive to the restoration of a healthy gut microbiota [ 25 ]. Through the analysis of these upregulated marker K -mers with K -motifs, we found that hormonal regulatory activity and carbohydrate metabolic processes were enriched in the post-treatment sample (Fig. S6). Further investigation into the shared species between the two samples, such as Phocaeicola dorei (Fig. 6 C), revealed that its populations in the pre- and post-treatment samples did not cluster together completely, indicating significant differences in gene expression. Consequently, we performed an in-depth analysis on functional changes in P. dorei between these two samples. GO enrichment results for marker K -mers in P. dorei (Fig. 6 D) showed that functions related to polysaccharide metabolism and outer membrane-associated defense responses were significantly upregulated in the post-treatment sample. Studies have shown that short-chain fatty acids produced by polysaccharide metabolism suggest efficient therapy and good prognosis for colorectal cancer [ 26 ]. In addition, outer membrane binding and periplasmic space can promote biofilm formation, enhancing the host’s immune defense [ 27 ]. Therefore, the innovative mKmer method can help researchers more comprehensively analyze the dynamic changes of intestinal microecology during immunotherapy, and identify potentially beneficial microorganisms or their metabolites, which is expected to improve the immunotherapy efficacy of solid tumors. Discussion Our study presents mKmer, a novel reference genome-free approach for analyzing msmRNA-seq data. The use of K -mer for species taxonomic identification has been demonstrated by many classical K -mer-based metagenomic taxonomy annotation software [ 9 , 10 , 11 ]. Although these tools classify species based on DNA, many conserved gene sequences remain highly consistent within species when DNA is transcripted into RNA. This implies that RNA sequences also contain species-specific conserved regions, and are useful in K -mer analysis for species identification. For instance, rRNA and tRNA are widely used in taxonomic studies due to their significant conservation and variation among species [ 28 – 29 ]. mRNA, on the other hand, reflects gene expression, which varies significantly between species. By analyzing high-frequency K -mers in mRNA, species-specific expression characteristics can be captured. In this study, we discovered the presence of HCKs in every single-cell sequencing dataset, and further used them as genic sequences for downstream msmRNA-seq analysis. By leveraging the strong correlation between marker K -mers and K -motifs, we further explored and obtained reliable results on the functions of the microorganisms. Compared to well-known tools such as Cell Ranger and STAR [ 3 ], mKmer overcomes the limitations of incomplete or poorly annotated reference genomes by utilizing a cell-by-HCK matrix instead of the traditional cell-by-gene matrix. Benchmark tests with soybean soil and human gut msmRNA-seq data demonstrated that mKmer captures more data than those available tools and achieves clearer species clustering. This is understandable, as both Cell Ranger and STAR align the obtained msmRNA-seq data to the available reference genomes which have been sequenced. This process inevitably introduces biases. Considering the rapid evolution and variation of microorganisms, using single reference genomes per species does not align with established facts; alignment failures due to genetic variations would result in discarding a significant amount of valuable biological information obtained from the msmRNA-seq. In contrast, mKmer, which operates without the need for reference genomes, provides a more comprehensive and unbiased analysis of microbiome samples. It is well known that scRNA-seq data contain numerous empty droplets and some doublets, as well as other impurity-containing droplets, which can severely impact the quality of sequencing results. Therefore, removing impurity information is crucial for single-cell analysis techniques. By examining the distribution of UMIs and barcodes, high-quality cells can be effectively filtered. Additionally, aligning to a reference genome to create a cell-by-gene matrix is an effective approach to exclude impurities from sequencing results. We know that the longer the K -mer, the higher its specificity; thus, longer K -mers are more efficient in detecting impurities. Conversely, shorter K -mers have higher conservation, which increases information utilization when identifying the same species. mKmer extracts biological information from raw sequencing data by testing different K -mer lengths to determine the optimal K. This optimal K can identify key conserved sequences while effectively distinguishing noise from non-critical gene sequences. As K -mer gets longer,, the sequencing results from each sample showed a consistent trend (e.g. Figure 2 ). Specifically, at a certain K , the number of K -mers with a frequency of 1 was the highest among all K -mer frequencies. We interpret K -mers with a frequency of 1 as sequences that are useless for species clustering and may even be impurities. To effectively filter out useless sequence, we ranked each K -mer in descending order of frequency and observed a distinct inflection point (Fig. 3 A and Fig. 3 C). Based on this, we consider the high-frequency K -mers before the inflection point as conservative K -mers (i.e., HCKs), and cells containing these K -mers are likely derived from a common ancestor (i.e., LCA) [ 30 ]. Therefore, HCKs exhibit high species recognition. The low-frequency K -mers after the inflection point may indicate that the taxonomic information contained in these K -mers is sparse, and they may even interfere with species identification and differentiation. It is not recommended to use low-frequency K -mers in the downstream dimensionality reduction and clustering process. It is undeniable that there may be errors near the inflection point, where some highly specific conserved K -mers may be misclassified as impurities. However, this has minimal impact on species identification across the entire cell set. In our analytical framework, we select these high-frequency conserved K -mers at inflection points for species identification and functional analysis. At the same time, selection of HCKs address issues of computational time and high matrix dimensionality caused by the excessive variety of K -mers, significantly enhancing mKmer's performance and practicality. Additionally, the high frequency of conserved K -mers means that they are more likely to come from regions that can translate conserved protein motifs. These protein motifs have been used for cross-species functional annotation of single-cell RNA sequencing (scRNA-seq) [ 31 ]. We further hypothesize that these HCKs may correspond to motif fragments of certain gene families within microorganisms. Therefore, by functionally annotating these gene motifs, mapping them onto the microbial communities, and performing enrichment analysis, we can infer the specific functions of the species in the sample.The key innovations of mKmer, including HCKs, marker K -mers, and K -motifs, enhance species identification and distinction. This method allows for a holistic view of microbial communities, advancing our understanding of microbial ecology and functional roles. Despite these strengths, mKmer still has room for improvement. For example, species annotation of Kraken 2 [ 32 ] can be corrected based on the results of clustering. Future research and development will further refine mKmer, for example, by expanding its functional modules, and optimizing its performance to provide more powerful support for microbiology research. Conclusions mKmer is a reference genome-free approach for msmRNA-seq analysis and allows studies on cellular heterogeneity, marker motif discovery, and efficiency functional annotation. In this study, we discovered and defined HCK. More accurate clustering results confirmed that HCKs densely encapsulate a large amount of hierarchical information from microbial taxonomy. In benchmark tests on soybean soil and human gut msmRNA-seq datasets, mKmer can use more msmRNA-seq data than the traditional annotated gene-based methods, achieving more and clearer species clustering for subsequent comprehensive functional analysis. Our method therefore provides an unbiased way to analyze all species in msmRNA-seq samples and allows diverse microbiomic single-cell problems to be formulated in a unified way. Methods msmRNA-seq data collection and generation A total of seven msmRNA-seq datasets, four from healthy donors by our previous study [ 1 ] (PRJCA017256), two from a patient and one from soybean soil generated by this study, were used for performance and benchmarking of mKmer. These two fecal samples were collected from the same patient with colorectal cancer both before and after immunotherapy. The study protocol was approved by the Ethics Committee of the First Affiliated Hospital, Zhejiang University School of Medicine, China (2021IIT A0239). The protocols for sample treatment for high-throughput msmRNA-seq followed our previous study [ 1 ] and the msmRNA-seq data generated by M20 Genomics were deposited at the NGDC database ( https://ngdc.cncb.ac.cn/ ) under accession code PRJCAXXX (publicly available as of the date of publication). Quality control of raw data UMI-tools (v1.1.4) [ 33 ] were used to process the unique molecular identifiers (UMIs) of our msmRNA-seq data. We utilized the umi-tools whitelist for quality control of the raw data to estimate the number of cells accurately. Given that the raw data file R1 is approximately 1GB, we specified an expected cell number before determining the actual count. Therefore, the --expect-cells parameter was set to 10,000. The raw data had a barcode length of 20bp and a UMI length of 8bp. To obtain the whitelist, the --bc-pattern was set to CCCCCCCCCCCCCCCCCCCC NNNNNNNN , and the --set-cell-number parameter was set to 7,000, which corresponded to the cell number at the inflection point in the barcode rank plot. We then used the umi-tools extract to filter the raw data files R1 and R2 based on the obtained whitelist, resulting in cleaned raw data. During the PCR process of smRNA-seq, some molecules may have been disproportionately amplified due to sequence characteristics (e.g., GC content) or random factors, resulting in multiple reads with the same barcode and UMI. To address this, we employed the RemoveDuplicates within UMI-tools to retain only the highest quality read, as determined by the Phred quality scoring system, among those with the same barcode and UMI. The sequencing data distinguishes each read in the form of 20bp_8bp, so the RemoveDuplicates defaults to the last 29 characters of the sequence information line in the FASTQ file as a unique identifier. K -mers scanning and counting Jellyfish (v2.2.10) [ 34 ] was used for fast, memory-efficient counting of K -mers in DNA sequences. In our experiment, the --m parameter of the jellyfish count was set to the default value of 10, to determine the K -mer length ( K ). This K corresponds to the smallest K where the peak in K -mer frequency distribution occurs at x = 1. The --s parameter was set to the default value of 10M. To ensure detection of all high-frequency K -mers, the --h parameter was set to 100,000,000 based on experimental testing. To confirm the selected K value was reasonable, we observed the distribution of peaks with different K values by drawing KmerFrequency plots. The --put of the KmerFrequency was three histo files with different K values specified, and the --out argument specified the path to the output KmerFrequency plot file. The histo file generated with the selected K was used as the input for the KmerRank to create a K -mer rank plot, where the x value at the inflection point indicates the number of HCKs. The jellyfish dump was then used to convert the jf file into a readable format for extracting the top-counts K -mers. HCK = select_top {sort [count (K-mer, C)]} M ij = count (HCK i , C j ) In the frequency matrix M , a specific element M ij represents the occurrence count of the i -th HCK in the j -th cell among the selected HCKs. Generation of cell-by-HCKs matrix After detecting each K -mer, they were sorted by detection depth in descending order. HCKs were then selected based on this sorted list. Using the generated HCK list, K -mers were read sequentially from the cleaned R2 reads. Each K -mer in every cell was counted to generate the cell-by-HCKs matrix. To minimize memory usage during execution, the program generated cell-by-HCKs matrix for every 1,000 cells read and then merged these matrices. For ease of use, we integrated the entire process of cell-by-HCKs matrix into a single command named KmerCell . The --kmercount argument required the kmer_counts_dumps.fa file output from jellyfish dump (the file suffix must be _counts_dumps.fa ); --fastq required the clean R2 FASTQ file; --topkmer specified the number of HCKs; and --k specified the selected K . Identification of microbial species We employed a K -mer-based root-to-leaf classification strategy for microbial species annotation, which is integrated in our software under the name smAnnotation . We first applied Kraken 2 (v 2.0.7-beta) [ 32 ], a K -mer-based read classification method, on every read in each barcode based on standard refseq of Kraken 2(Refseq archaea, bacteria, viral, plasmid, human1, https://benlangmead.github.io/aws-indexes/k2 ). After all the reads were assigned into each node of different taxonomic levels (e.g., order, family, genus, species), we calculate the sum of reads in each node from the leaf to the root. Then we performed the taxonomic classification from the root to leaf taxonomic levels. In the root taxonomic level, we ranked all the nodes from the highest to lowest based on the number of reads of the nodes, and selected the node with the most reads as a potential annotation candidate. Based on the annotation results, then we performed the same annotation process in the next lower taxonomic level, until the leaf nodes (species level). Then Bracken (v 2.5) [ 35 ] was used to count the fraction_total_reads of the species classified into each cell, and the species with the largest value was selected as the final annotation result. For Kraken 2, the comparison database was the NCBI standard database by default (Archive size: 60GB) and resolution parameter --r was set to 100. For smAnnotation , clean R2 as the specified file for --input , and the output file named smAnnotation.report was placed in the current working directory by default. Visualization and clustering To visualize the data, we further reduced the dimensionality of all filtered cells using Seurat (v4) [ 36 ] and used UMAP to project the cells into 2D space. The annotation results of Kraken 2 and Bracken were mapped to the Seurat object; only the annotation results with fraction_total_reads values greater than 0.5 were retained, and strains with abundance less than 0.1% were filtered out. The steps include: (i) Using the LogNormalize method of the NormalizeData of Seurat to calculate the expression values of K -mers. The scale.factor argument is set to the default 10000, nfeatures to 6000, and the ScaleData object to all genes; (ii) PCA was performed using the normalized expression value; among all the principal components, the top 30 principal components were used to do clustering and UMAP analysis; (iii) To find clusters, a weighted graph-based clustering method, Shared Nearest Neighbour (SNN), was selected, and the resolution is set to 0.5. Marker genes for each cluster were identified with the MAESTRO test with default parameters via the FindAllMarkers in Seurat and the min.pct parameter was set to 0.25. GO annotation by marker K -mers and K -motifs Before performing functional enrichment analysis on sequencing data from different samples, we first need to filter out rRNA from the raw transcriptome data. In this study, we used SortMeRNA (v4.2.0) [ 37 ] for filtering, with the reference dataset including multiple species, such as bacteria, eukaryotes, archaea, and various sequencing databases, for 5S, 16S, 23S, and other rRNA types. After using the FindAllMarkers , we obtained a list of marker K -mers, which were used for functional enrichment analysis of clusters or species of interest. Marker K -mers are considered to be identified from highly conserved sequences, which are likely to represent individual motifs. The functional analysis of bacterial species based on their specific motifs is reliable. The MEME (Multiple Em for Motif Elicitation) suite (v5.0.5) [ 38 ] is a comprehensive resource for discovering and analyzing sequence motifs in DNA, RNA, and protein sequences. Memes [ 39 ] is an R package that provides a seamless R interface to a selection of popular MEME Suite tools. By analyzing the conserved sequences of each strain, we aimed to elucidate the specific functions of the strains. To obtain the specific functions of bacterial species, we designed two partitioning workflows to conduct Gene Ontology (GO) enrichment analysis on K -motifs. Nucleotide motif analysis (KmerGOn) The first analysis workflow involved converting each marker K -mer into a motif file in MEME format. These motifs were then compared against motifs in the microbial nucleotide motif database using the tomtom [ 40 ] integrated into the MEME suite. This step identified the best-matching known motifs. Subsequently, ama and gomo [ 41 ] in the MEME suite were used to compare the identified known motifs against the Escherichia coli database, obtaining GO functions for each motif. Finally, GO functional enrichment and visualization were performed on the clusters or species of interest. Protein motif analysis (KmerGOp) The second analysis workflow utilized SeqKit (v2.8.2) [ 42 , 43 ] to translate each marker K -mer into amino acids using six reading frames (since K = 12, only sequences with 4 AA were retained). These motifs were then compared against motifs in the all-species motif database using the tomtom integrated into the MEME suite. Each protein motif was further analyzed using InterProScan (v5.47-82.0) [ 44 ] to search domain databases (e.g., Pfam, PROSITE, PRINTS, etc.) and obtain GO IDs. Finally, used select of the AnnotationDbi (v1.64.1) [ 45 ] package to match the corresponding term and GO ID from the GO.db (v3.18.0) [ 46 ] database. For both workflows, the output list of marker K -mers from the FindAllMarkers served as the input file. The -- cluster was set to the target cluster, and the -- out specified the path to the output file. This ensured a systematic approach to uncovering the functional roles of conserved sequences within bacterial species. Declarations Supplementary Information The online version contains supplementary material available at XXX. Additional file 1: Includes supplementary figures Fig S1—Fig S6. A short description of the content of these figures is provided at the first page of the file. Peer review information XXX was the primary editor of this article at Genome Biology and managed its editorial process and peer review in collaboration with the rest of the editorial team. Review history The review history is available as XXX. Authors’ contributions LF and WJ conceived the study. WJ, WC, YW, WC and XZ conducted the experiments. FM, QQ, XL, DZ, JY, HC, YH and SW analyzed the data. FM and LF developed mKmer. FM, LF and YB wrote the paper. LF, WJ, YW and YS discussed and supervised this project. All authors have revised and approved the final manuscript. Funding This study was supported by Biological Breeding-Major (2023ZD04076), Yunnan Tobacco Company (2024530000241001) and CIC-MIC. Data availability All msmRNA-seq data used by this study are available at NGDC database (https://ngdc.cncb.ac.cn/) under project number PRJCA017256 (accession number SAMC3766839, SAMC3766838, SAMC3766837, SAMC1266599), PRJCAXXX (the soil sample) and PRJCAXXX (the two gut samples) (publicly available as of the date of publication). The mKmer package (v.1.0.0) is available at https://github.com/bioinplant/mKmer. Any additional information required to reanalyze the data reported in this work paper is available from the lead contact upon request. Ethics approval and consent to participate All data collection was approved by the Ethics Committee of the First Affiliated Hospital, Zhejiang University School of Medicine, China (2021IIT A0239). Consent for publication Not applicable. Competing interests The authors declare no competing interests References Shen Y, Qian Q, Ding L, Qu W, Zhang T, Song M et al (2024) High-throughput single-microbe RNA sequencing reveals adaptive state heterogeneity and host-phage activity associations in human gut microbiome. Protein Cell. ;pwae027 Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M et al (2015) Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell 161(5):1202–1214 Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S et al (2013) STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29(1):15–21 Eren AM, Delmont TO (2024) Bioprospecting marine microbial genomes to improve biotechnology. Nature 633(8029):287–288 Chen H, Ryu J, Vinyard ME, Lerer A, Pinello L (2024) SIMBA: single-cell embedding along with features. Nat Methods 21(6):1003–1013 Tayyebi Z, Pine AR, Leslie CS (2024) Scalable and unbiased sequence-informed embedding of single-cell ATAC-seq data with CellSpace. Nat Methods 21(6):1014–1022 Kim J, Steinegger M (2024) Metabuli: sensitive and specific metagenomic classification via joint analysis of amino acid and DNA. Nat Methods 21(6):971–973 Kwok AYC, Su S-C, Reynolds RP, Bay SJ, Av-Gay Y, Dovichi NJ et al (1999) Species identification and phylogenetic relationships based on partial HSP60 gene sequences within the genus Staphylococcus. Int J Syst Evol MicroBiol 49:1181–1192 Wood DE, Salzberg SL (2014) Kraken: ultrafast metagenomic sequence classification using exact alignments. Genome Biol 15(3):R46 Kim D, Song L, Breitwieser FP, Salzberg SL (2016) Centrifuge: rapid and sensitive classification of metagenomic sequences. Genome Res 26(12):1721–1729 Menzel P, Ng KL, Krogh A (2016) Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat Commun 7:11257 Sharma M, Khurana H, Singh DN, Negi RK (2021) The genus Sphingopyxis: Systematics, ecology, and bioremediation potential - A review. J Environ Manage 280:111744 Kim Y-S, Hwang E-M, Jeong C-M, Cha C-J (2023) Flavobacterium psychrotrophum sp. nov. and Flavobacterium panacagri sp. nov., Isolated from Freshwater and Soil. J Microbiol 61(10):891–901 Männistö MK, Tiirola M, McConnell J, Häggblom MM (2010) Mucilaginibacter frigoritolerans sp. nov., Mucilaginibacter lappiensis sp. nov. and Mucilaginibacter mallensis sp. nov., isolated from soil and lichen samples. Int J Syst Evol MicroBiol 60(Pt 12):2849–2856 Kim M-J, Park Y-J, Park M-K, Park CE, Jo Y, Tagele SB et al (2020) Complete genome sequence of Labrys sp. KNU-23 isolated from ginseng soil in the Republic of Korea. Korean J Microbiol 56(4):410–412 Abdugheni R, Wang W, Wang Y, Du M, Liu F, Zhou N et al (2022) Metabolite profiling of human-originated Lachnospiraceae at the strain level. iMeta 1(4):e58 Ramirez-Farias C, Slezak K, Fuller Z, Duncan A, Holtrop G, Louis P (2008) Effect of inulin on the human gut microbiota: stimulation of Bifidobacterium adolescentis and Faecalibacterium prausnitzii. Br J Nutr 101(4):541–550 Ichimura R, Tanaka K, Nakato G, Fukuda S, Arakawa K (2024) Complete genome sequence of Mediterraneibacter gnavus strain RI1, isolated from human feces. Microbiol Resour Announc 0:e00863–e00824 Li H, Xu H, Li Y, Jiang Y, Hu Y, Liu T et al (2020) Alterations of gut microbiota contribute to the progression of unruptured intracranial aneurysms. Nat Commun 11(1):3218 Wang K, Liao M, Zhou N, Bao L, Ma K, Zheng Z et al (2019) Parabacteroides distasonis Alleviates Obesity and Metabolic Dysfunctions via Production of Succinate and Secondary Bile Acids. Cell Rep 26(1):222–235e5 1, Wu F, Guo X, Zhang J, Zhang M, Ou Z, Peng Y (2017) Phascolarctobacterium faecium abundant colonization in human gastrointestinal tract. Experimental Therapeutic Med 14(4):3122–3126 García-López M, Meier-Kolthoff JP, Tindall BJ, Gronow S, Woyke T, Kyrpides NC et al (2019) Analysis of 1,000 Type-Strain Genomes Improves Taxonomic Classification of Bacteroidetes. Front Microbiol 10:2083 Leth ML, Ejby M, Workman C, Ewald DA, Pedersen SS, Sternberg C et al (2018) Differential bacterial capture and transport preferences facilitate co-growth on dietary xylan in the human gut. Nat Microbiol 3(5):570–580 Williams BL, Hornig M, Parekh T, Lipkin WI (2012) Application of Novel PCR-Based Methods for Detection, Quantitation, and Phylogenetic Characterization of Sutterella Species in Intestinal Biopsy Samples from Children with Autism and Gastrointestinal Disturbances. mBio 3(1):e00261–e00211 Lozupone CA, Stombaugh JI, Gordon JI, Jansson JK, Knight R (2012) Diversity, stability and resilience of the human gut microbiota. Nature 489(7415):220–230 Louis P, Flint HJ (2017) Formation of propionate and butyrate by the human colonic microbiota. Environ Microbiol 19(1):29–41 Costerton JW, Stewart PS, Greenberg EP (1999) Bacterial Biofilms: A Common Cause of Persistent Infections. Science 284(5418):1318–1322 Janda JM, Abbott SL (2007) 16S rRNA Gene Sequencing for Bacterial Identification in the Diagnostic Laboratory: Pluses, Perils, and Pitfalls. J Clincal Microbiol 45(9):2761–2764 Bermudez-Santana C, Attolini CS-O, Kirsten T, Engelhardt J, Prohaska SJ, Steigele S et al (2010) Genomic organization of eukaryotic tRNAs. BMC Genomics 11:270 Nasko DJ, Koren S, Phillippy AM, Treangen TJ (2018) RefSeq database growth influences the accuracy of k-mer-based lowest common ancestor species identification. Genome Biol 19(1):165 Rosen Y, Brbić M, Roohani Y, Swanson K, Li Z, Leskovec J (2024) Toward universal cell embeddings: integrating single-cell RNA-seq datasets across species with SATURN. Nat Methods 21(8):1492–1500 Wood DE, Lu J, Langmead B (2019) Improved metagenomic analysis with Kraken 2. Genome Biol 20(1):257 Smith T, Heger A, Sudbery I (2017) UMI-tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Res 27(3):491–499 Marçais G, Kingsford C (2011) A fast, lock-free approach for efficient parallel counting of occurrences of k -mers. Bioinformatics 27(6):764–770 Lu J, Breitwieser FP, Thielen P, Salzberg SL (2017) Bracken: estimating species abundance in metagenomics data. PeerJ Comput Sci 3:e104 Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A et al (2021) Integrated analysis of multimodal single-cell data. Cell 184(13):3573–3587e29 Kopylova E, Noé L, Touzet H (2012) SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics 28(24):3211–3217 Bailey TL, Johnson J, Grant CE, Noble WS (2015) The MEME Suite. Nucleic Acids Res 43(W1):W39–49 Nystrom SL, McKay DJ, Memes (2021) A motif analysis environment in R using tools from the MEME Suite. Pertea M, editor. PLOS Computational Biology. ;17(9):e1008991 Gupta S, Stamatoyannopoulos JA, Bailey TL, Noble W (2007) Quantifying similarity between motifs. Genome Biol 8(2):R24 Buske FA, Bodén M, Bauer DC, Bailey TL (2010) Assigning roles to DNA regulatory motifs using comparative genomics. Bioinformatics 26(7):860–866 Shen W, Le S, Li Y, Hu F, SeqKit: (2016) A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation. Zou Q, editor. PLOS ONE. ;11(10):e0163962 Shen W, Sipos B, Zhao L (2024) SeqKit2: A Swiss army knife for sequence and alignment processing. iMeta 3(3):e191 Jones P, Binns D, Chang H-Y, Fraser M, Li W, McAnulla C et al (2014) InterProScan 5: genome-scale protein function classification. Bioinformatics 30(9):1236–1240 Pagès H, Carlson M, Falcon S, Nian L, AnnotationDbi (2023) Manipulation of SQLite-based annotations in Bioconductor. Bioconductor. 10.18129/B9.bioc.AnnotationDbi The Gene Ontology Consortium, Carbon S, Douglass E, Good BM, Unni DR, Harris NL et al (2021) The Gene Ontology resource: enriching a GOld mine. Nucleic Acids Res 49(D1):D325–D334 Additional Declarations The authors declare no competing interests. Supplementary Files SupplementaryFiles.zip Includes supplementary figures Fig S1—Fig S6. A short description of the content of these figures is provided at the first page of the file. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-5748035","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Method Article","associatedPublications":[],"authors":[{"id":396498238,"identity":"752e976d-61eb-4e71-93f1-94dcec680eed","order_by":0,"name":"Fangyu Mo","email":"","orcid":"https://orcid.org/0009-0003-0560-5551","institution":"Hainan Institute, Zhejiang University, Sanya 572025, China","correspondingAuthor":false,"prefix":"","firstName":"Fangyu","middleName":"","lastName":"Mo","suffix":""},{"id":396498239,"identity":"b1fc11a5-3323-4983-be49-83433de26e09","order_by":1,"name":"Qinghong Qian","email":"","orcid":"","institution":"Institute of Crop Science, Zhejiang University, Hangzhou 310058, China","correspondingAuthor":false,"prefix":"","firstName":"Qinghong","middleName":"","lastName":"Qian","suffix":""},{"id":396498240,"identity":"e09ecb2e-aebb-4897-9149-768d999a7a1d","order_by":2,"name":"Xiaolin Lu","email":"","orcid":"","institution":"Institute of Bioinformatics and James D. Watson Institute of Genome Sciences, Zhejiang University, Hangzhou 310058, China.","correspondingAuthor":false,"prefix":"","firstName":"Xiaolin","middleName":"","lastName":"Lu","suffix":""},{"id":396498241,"identity":"822cd419-1bad-497d-a591-717540dcc0fb","order_by":3,"name":"Dihuai Zheng","email":"","orcid":"","institution":"Institute of Crop Science, Zhejiang University, Hangzhou 310058, China","correspondingAuthor":false,"prefix":"","firstName":"Dihuai","middleName":"","lastName":"Zheng","suffix":""},{"id":396498242,"identity":"650cbc12-34f7-40e9-8cee-9c1ec8dea79f","order_by":4,"name":"Wenjie Cai","email":"","orcid":"","institution":"Liangzhu Laboratory, Zhejiang University Medical Center, Hangzhou, China","correspondingAuthor":false,"prefix":"","firstName":"Wenjie","middleName":"","lastName":"Cai","suffix":""},{"id":396498243,"identity":"108eee94-6ad7-46d8-bb27-155155bd2fa7","order_by":5,"name":"Jie Yao","email":"","orcid":"","institution":"Institute of Bioinformatics and James D. Watson Institute of Genome Sciences, Zhejiang University, Hangzhou 310058, China.","correspondingAuthor":false,"prefix":"","firstName":"Jie","middleName":"","lastName":"Yao","suffix":""},{"id":396498244,"identity":"7267c49f-8967-4ef9-b2b9-541f6b19da16","order_by":6,"name":"Hongyu Chen","email":"","orcid":"","institution":"Institute of Crop Science, Zhejiang University, Hangzhou 310058, China","correspondingAuthor":false,"prefix":"","firstName":"Hongyu","middleName":"","lastName":"Chen","suffix":""},{"id":396498245,"identity":"a6a19d56-c9cc-47da-bc91-47e54edf8f3b","order_by":7,"name":"Yujie Huang","email":"","orcid":"","institution":"Institute of Crop Science, Zhejiang University, Hangzhou 310058, China","correspondingAuthor":false,"prefix":"","firstName":"Yujie","middleName":"","lastName":"Huang","suffix":""},{"id":396498246,"identity":"4771612d-e80c-4d81-b695-1234e73fd639","order_by":8,"name":"Xiang Zhang","email":"","orcid":"","institution":"Department of Colorectal Surgery, the First Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou 310003, China","correspondingAuthor":false,"prefix":"","firstName":"Xiang","middleName":"","lastName":"Zhang","suffix":""},{"id":396498247,"identity":"e6647571-bf84-4434-9159-70e655ff53f1","order_by":9,"name":"Sanling Wu","email":"","orcid":"","institution":"Analysis Center of Agrobiology and Environmental Sciences, Faculty of Agriculture, Life and Environment Sciences, Zhejiang University, Hangzhou 310058, China","correspondingAuthor":false,"prefix":"","firstName":"Sanling","middleName":"","lastName":"Wu","suffix":""},{"id":396498248,"identity":"b696421f-06b6-4c7c-8330-21edab9512a1","order_by":10,"name":"Yifei Shen","email":"","orcid":"","institution":"Institute of Department of Laboratory Medicine, the First Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou 310003, China","correspondingAuthor":false,"prefix":"","firstName":"Yifei","middleName":"","lastName":"Shen","suffix":""},{"id":396498249,"identity":"d0ae4410-cd28-4182-ac5b-b58e5bc70cc8","order_by":11,"name":"Yingqi Bai","email":"","orcid":"","institution":"BGI-Sanya, Sanya 572025, China","correspondingAuthor":false,"prefix":"","firstName":"Yingqi","middleName":"","lastName":"Bai","suffix":""},{"id":396498250,"identity":"cc8e28ca-2b64-4121-9a73-8bbe010ccce7","order_by":12,"name":"Yongcheng Wang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA4klEQVRIiWNgGAWjYLACxgYgwd7AYABjE6mF5wDJWiQSEGy8wOB47+EXP3fY5MlHPn9QzMNgI7vhAPOzB3i1nDmXZtl7Jq3Y8HZCgjEPQ5rxhgNs5gZ4tdzIMTNmbDucuHF2wgGglsOJGw7wsEkQoeV/4saZBxuAWv4TpcX4MWPbgcT5EswMQC0HCGuRPHPGjLG3LTlxA08ag+Ecg2TjmYfZzPBq4TveY/zhZ5td4vz2488M3lTYyfYdb36GV4vCAQaIMwyADANwZDLjUw8E8g0MzB9gjAcEFI+CUTAKRsEIBQC6UE5QNa3W2wAAAABJRU5ErkJggg==","orcid":"","institution":"Liangzhu Laboratory, Zhejiang University, Hangzhou 311113, China","correspondingAuthor":true,"prefix":"","firstName":"Yongcheng","middleName":"","lastName":"Wang","suffix":""},{"id":396498251,"identity":"38b1fd7c-694c-4dd8-9e58-2b8d1e9bc2a5","order_by":13,"name":"Weiqin Jiang","email":"","orcid":"","institution":"Department of Colorectal Surgery, the First Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou 310003, China","correspondingAuthor":false,"prefix":"","firstName":"Weiqin","middleName":"","lastName":"Jiang","suffix":""},{"id":396498252,"identity":"305f18b4-3e9f-416d-b0ee-a09fb4aa3f68","order_by":14,"name":"Longjiang Fan","email":"","orcid":"","institution":"Zhejiang University","correspondingAuthor":false,"prefix":"","firstName":"Longjiang","middleName":"","lastName":"Fan","suffix":""}],"badges":[],"createdAt":"2025-01-02 01:33:08","currentVersionCode":1,"declarations":{"humanSubjects":false,"vertebrateSubjects":false,"conflictsOfInterestStatement":false,"humanSubjectEthicalGuidelines":false,"humanSubjectConsent":false,"humanSubjectClinicalTrial":false,"humanSubjectCaseReport":false,"vertebrateSubjectEthicalGuidelines":false},"doi":"10.21203/rs.3.rs-5748035/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-5748035/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":72871692,"identity":"92f903ee-4b72-4018-8ed2-a9d48f5a6401","added_by":"auto","created_at":"2025-01-03 07:16:44","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":6349742,"visible":true,"origin":"","legend":"\u003cp\u003eOverview of the mKmer method. There are seven modules in mKmer for msmRNA-seq data analysis, including reads selection, species annotation, selection of \u003cem\u003eK\u003c/em\u003e and HCKs, cell-by-HCK matrix construction, and functional analysis. Traditional analysis software, such as UMI-tools, Jellyfish, and Seurat, are also employed in the mKmer pipeline.\u003c/p\u003e","description":"","filename":"Figure1.png","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/5097c466dbf3183299f24545.png"},{"id":72871689,"identity":"c3f3e30b-61cf-4acc-a7f9-354ba6cb5c3c","added_by":"auto","created_at":"2025-01-03 07:16:44","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":2492551,"visible":true,"origin":"","legend":"\u003cp\u003eFrequency of \u003cem\u003eK\u003c/em\u003e-mers depths by different \u003cem\u003eK\u003c/em\u003e sizes for msmRNA-seq datasets from soybean soil (\u003cstrong\u003eA\u003c/strong\u003e) and human gut (\u003cstrong\u003eB\u003c/strong\u003e) samples.\u003c/p\u003e","description":"","filename":"Figure2.png","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/c9d60d22ebc31a00731397fc.png"},{"id":72871696,"identity":"5b27135a-56d9-42bc-b3c8-b9fee51f03db","added_by":"auto","created_at":"2025-01-03 07:16:44","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":6226027,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eK\u003c/em\u003e-mer rank plot for feature \u003cem\u003eK\u003c/em\u003e-mers (HCKs) and barcode rank plot calling of the msmRNA sample. The Y-axis coordinates represent UMI count. \u003cstrong\u003eA\u003c/strong\u003e \u003cem\u003eK\u003c/em\u003e-mer rank plot of a soybean soil sample (\u003cem\u003eK \u003c/em\u003e= 13). \u003cstrong\u003eB\u003c/strong\u003e Barcode rank plot of the soybean soil data. \u003cstrong\u003eC\u003c/strong\u003e \u003cem\u003eK\u003c/em\u003e-mer rank plot of a human gut sample (\u003cem\u003eK\u003c/em\u003e = 12). \u003cstrong\u003eD\u003c/strong\u003e Barcode rank plot of the human gut data.\u003c/p\u003e","description":"","filename":"Figure3.png","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/37dec729b6ea04b542897f5b.png"},{"id":72871699,"identity":"f769d35a-15d5-40fa-937c-1235cf6e1462","added_by":"auto","created_at":"2025-01-03 07:16:45","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":44531584,"visible":true,"origin":"","legend":"\u003cp\u003eAn example of HCKs. A conserved region coding for a motif of the bacterial gene \u003cem\u003eHSP60\u003c/em\u003e (top), and its counts of 12-mers obtained by scanning \u003cem\u003eStaphylococcus\u003c/em\u003e, displayed as a line chart to show the distribution of HCKs (bottom). The motif’s names (nuclear/protein) from the MEME database are shown at the top. In the middle, the 12-mer composition of the conserved region (\u003cem\u003eS. warned\u003c/em\u003e as reference) is shown. The numbers of each 12-mer are provided at the end of 12-mers.\u003c/p\u003e","description":"","filename":"Figure4.png","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/baa5400c840f857f24990a10.png"},{"id":72871697,"identity":"9ff3e6eb-0f4f-4e29-9e27-7b7e9fee68ce","added_by":"auto","created_at":"2025-01-03 07:16:44","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":15152670,"visible":true,"origin":"","legend":"\u003cp\u003eBenchmark results of two samples by cell-by-gene matrix (left) and mKmer (right). UMAP clustering and species annotation of the msmRNA-seq data from the soybean soil (\u003cstrong\u003eA\u003c/strong\u003e) and human gut (\u003cstrong\u003eB\u003c/strong\u003e) samples.\u003c/p\u003e","description":"","filename":"Figure5.png","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/75ff4753781f2857566ec486.png"},{"id":72871698,"identity":"0e26c041-ec06-44dd-b73a-7a1a82bd477d","added_by":"auto","created_at":"2025-01-03 07:16:45","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":14780982,"visible":true,"origin":"","legend":"\u003cp\u003eA case study by mKmer. \u003cstrong\u003eA-B\u003c/strong\u003eUMAP clustering and species annotation of the pre-treatment (\u003cstrong\u003eA\u003c/strong\u003e) and post-treatment (\u003cstrong\u003eB\u003c/strong\u003e) gut msmRNA-seq data of a cancer patient using \u003cem\u003eK\u003c/em\u003e-mers. \u003cstrong\u003eC\u003c/strong\u003e Integrated UMAP clustering of bacterial species \u003cem\u003eP. dorei\u003c/em\u003e in the patient’s gut before (blue) and after treatment (red) with mKmer. \u003cstrong\u003eD\u003c/strong\u003eFunctional annotation of \u003cem\u003eP. dorei\u003c/em\u003e in the gut of treated patients by \u003cem\u003eKmerGOp\u003c/em\u003e.\u003c/p\u003e","description":"","filename":"Figure6.png","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/5645f424f15cabefcb87fa0c.png"},{"id":72873564,"identity":"3cd6b09d-b5a5-41df-acfa-82265f65ca1a","added_by":"auto","created_at":"2025-01-03 07:41:44","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":86114993,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/f2fd928a-dc01-4e81-829c-43e041799ccb.pdf"},{"id":72871693,"identity":"99d4640a-1961-4ec9-8a47-26fa5a9fe4a5","added_by":"auto","created_at":"2025-01-03 07:16:44","extension":"zip","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":28704547,"visible":true,"origin":"","legend":"\u003cp\u003eIncludes supplementary figures Fig S1—Fig S6. A short description of the content of these figures is provided at the first page of the file.\u003c/p\u003e","description":"","filename":"SupplementaryFiles.zip","url":"https://assets-eu.researchsquare.com/files/rs-5748035/v1/9a52fd0cc17f00b4984f47b0.zip"}],"financialInterests":"The authors declare no competing interests.","formattedTitle":"\u003cp\u003e\u003cstrong\u003emKmer: An unbiased \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eK\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e-mer embedding of microbiomic single-microbe RNA\u003c/strong\u003e \u003cstrong\u003esequencing data\u003c/strong\u003e\u003c/p\u003e","fulltext":[{"header":"Background","content":"\u003cp\u003eRecently, we developed a high-throughput single-microbe RNA sequencing (smRNA-seq) technique for microbiome samples [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e], which can generate RNAs from over 5,000 single microbes of a microbial community. This technique can effectively solve the problem of significant cell heterogeneity among bacterial populations, and thus achieve complete functional characterization of host-related microorganisms. However, there is a key challenge in the microbiomics smRNA-seq (msmRNA-seq) data analysis: constructing a high-quality gene expression matrix for downstream analysis [\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. Typically, the construction of a gene expression matrix requires a reference genome [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. While the human reference genome is largely complete and accurate, this is not the case for other organisms. In particular, many species in a microbiome sample lack reference genomes or do not have a high-quality annotated gene set [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]. These issues are especially pronounced when sequencing data encompasses multiple unclassified species (e.g., microbiome data). Any attempt to predefine which species should be included in the reference genome(s) inevitably introduces bias.\u003c/p\u003e \u003cp\u003e \u003cem\u003eK\u003c/em\u003e-mer refers to short sequences of a specific length (\u003cem\u003eK\u003c/em\u003e) and any genomic or RNA sequences are composed of different \u003cem\u003eK\u003c/em\u003e-mers. \u003cem\u003eK\u003c/em\u003e-mers have been widely used in bioinformatics analysis, including genome survey and assembly, and also single-cell omics data [e.g. 5, 6, 7]. Theoretically, RNAs from a microbiome sample can be characterized by a specific length of short sequence (i.e., \u003cem\u003eK\u003c/em\u003e-mers). In this study, we developed a new frame of msmRNA-seq analysis (named mKmer) based on high-frequency conserved \u003cem\u003eK\u003c/em\u003e-mer (HCK) rather than a gene expression matrix. Benchmark tests on seven datasets from soil and gut demonstrate that mKmer significantly improved species identification compared to the cell-by-gene matrix. To demonstrate applications of our method, we used a clinical msmRNA-seq data from the gut microbiome of colorectal cancer patients before and after treatment.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eOverview of mKmer method\u003c/h2\u003e \u003cp\u003emKmer is a tool for msmRNA data analysis by constructing a cell-by-HCK matrix to achieve efficient cell classification, species annotation, and functional analysis. The tool is currently divided into seven analysis modules, including \u003cem\u003eKmerRank\u003c/em\u003e and \u003cem\u003eKmerCell\u003c/em\u003e, aiming to provide personalized services for analysts (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). mKmer extracts biological information from raw sequencing data by testing different \u003cem\u003eK\u003c/em\u003e-mer lengths to determine the optimal \u003cem\u003eK\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). This optimal \u003cem\u003eK\u003c/em\u003e can identify key conserved sequences while effectively distinguishing noise from non-critical gene sequences. In our analysis framework, we successfully identified species and performed functional analysis for soybean soil samples and human fecal samples by selecting these HCKs before the inflection point. This demonstrates that the HCKs already densely contain a large amount of hierarchical information in microbial taxonomy, enabling precise and efficient classification at various levels, including domain, phylum, class, order, family, genus, and species. On the other hand, the low-frequency \u003cem\u003eK\u003c/em\u003e-mers contain sparse taxonomic information, and discarding them does not affect the classification results. Therefore, the HCKs means that may be translated into amino acid sequences (i.e. protein motifs) accordingly. Additionally, HCKs address issues of computational time and high matrix dimensionality caused by the excessive variety of \u003cem\u003eK\u003c/em\u003e-mers, which significantly enhanced mKmer's performance and practicality. Because it does not require any reference genome, the reproducibility of sequencing reads reached 100%.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eK-mer scanning and rank plot\u003c/h3\u003e\n\u003cp\u003eFor a raw msmRNA sequencing dataset, cells were selected as usual (such as with UMI-tools), and reads from the selected cells were filtered to remove duplicates before the downstream \u003cem\u003eK\u003c/em\u003e-mers scanning (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). We tested frequency distributions of different \u003cem\u003eK\u003c/em\u003e-mers sizes for the seven msmRNA datasets and found a change of distribution curves happening between 12-mer and 13-mer in all samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e, other five samples see Fig. S1). Due to the high-throughput smRNA-seq technique using random primers to capture RNA of individual cells, the combined amplification and release yield an average RNA sequencing depth of 1\u0026times; coverage. Therefore, low-frequency \u003cem\u003eK\u003c/em\u003e-mers with a frequency of 1, where the \u003cem\u003eK\u003c/em\u003e value is at its peak, are considered optimal. The 12-mer was therefore used as the default size for \u003cem\u003eK\u003c/em\u003e-mers scanning. We further ranked all 12-mers or 13-mers scanned from the msmRNA data by count per cell (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, other five samples see Fig. S2). The \u003cem\u003eK\u003c/em\u003e-mer rank plot (from highest to lowest \u003cem\u003eK\u003c/em\u003e-mer depth) is an interactive plot that shows all \u003cem\u003eK\u003c/em\u003e-mers detected in a microbiome sample or a msmRNA-seq dataset.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e\n\u003ch3\u003eIdentification of HCKs\u003c/h3\u003e\n\u003cp\u003eThe overall shape of the \u003cem\u003eK\u003c/em\u003e-mer rank plot (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, left) is similar to the barcode rank plot (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, right). Typically, a \u0026ldquo;cliff-and-knee\u0026rdquo; shape can be observed in the \u003cem\u003eK\u003c/em\u003e-mer rank plot of a microbiome sample. In this case, the steep cliff, followed by the plateaued knee, demonstrates that the \u003cem\u003eK\u003c/em\u003e-mer calling algorithm was able to distinguish feature \u003cem\u003eK\u003c/em\u003e-mers from others. HCKs mainly come from the evolutionary conserved regions (e.g., motifs in protein domains and DNA-binding sites) of bacterial species in a microbiome sample. As an example, the region of the genus \u003cem\u003eStaphylococcus HSP60\u003c/em\u003e gene encodes the conserved NdhRMIQE motif [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e], and a high number of \u003cem\u003eK\u003c/em\u003e-mers could be counted within this region (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). The conserved \u003cem\u003eK\u003c/em\u003e-mers exist in a wide variety of microbe species in a microbiome. From a microbial taxonomy perspective, HCK corresponds to the lowest common ancestor (LCA) sequence at the taxonomic level. The use of \u003cem\u003eK\u003c/em\u003e-mers for microbial classification has been demonstrated by many classical \u003cem\u003eK\u003c/em\u003e-mer-based metagenomic taxonomy annotation software, such as Kraken [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e], Centrifuge [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e], and Kaiju [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]. HCK was first discovered and successfully applied in the identification of taxa in single-cell data.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e\n\u003ch3\u003eIdentificaiton of maker K-mers\u003c/h3\u003e\n\u003cp\u003eBased on the cell-by-HCK matrix and routine cell clustering and dimension reduction (as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e), a visualization result by uniform manifold approximation and projection (UMAP) of a msmRNA sample and species annotation can be obtained (examples shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e, right, and the cell-by-gene matrix clustering results of this sample are shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e, left). A good clustering of the same species/cells was observed in the UMAP plot. Further, marker \u003cem\u003eK\u003c/em\u003e-mers can be identified among the different clusters (species or subspecies) using routine approaches (same as those for marker genes) such as the Seurat function (\u003cem\u003eFindAllMarkers\u003c/em\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e\n\u003ch3\u003eFunction annotation by K-motifs\u003c/h3\u003e\n\u003cp\u003eThe marker \u003cem\u003eK\u003c/em\u003e-mer can be used for functional annotation based on their encoding motifs as mentioned above. Based on motif and domain databases (MEME and Pfam), the DNA motifs and protein motifs in domains (termed the \u003cem\u003eK\u003c/em\u003e-mer-contained motifs as \u003cem\u003eK\u003c/em\u003e-motifs) can be identified for gene ontology (GO) annotation by mKmer functions (\u003cem\u003eKmerGOn and KmerGOp\u003c/em\u003e), respectively. At the protein level, the longest translated amino acids (AA) (4-mer AA) for marker \u003cem\u003eK\u003c/em\u003e-mers (12-mer nt) in a microbiome sample should have significant sequence similarity to the protein domain\u0026rsquo;s motifs in Pfam. The motif-contained \u003cem\u003eK\u003c/em\u003e-mers are those highly conserved \u003cem\u003eK\u003c/em\u003e-mers which transcript from the motif-contained regions of orthologous genes of different microbe species in a microbiome sample. Using the \u003cem\u003eK\u003c/em\u003e-motifs and their GO IDs identified, function analysis such as GO and pathway enrichment can be done as usual (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eD).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eBenchmark test with cell-by-gene matrix\u003c/h2\u003e \u003cp\u003eTo compare the performance of mKmer with the traditional gene matrix-based method, we generated a msmRNA-seq dataset from soybean (\u003cem\u003eGlycine max\u003c/em\u003e) soil and collected four publicly available msmRNA-seq datasets from human guts. Firstly, when comparing the dimensionality reduction and clustering results using the traditional cell-by-gene matrix (two examples are shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e, left), regardless of whether the msmRNA-seq data came from soil (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA) or human gut (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB), the cell-by-HCK matrix (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e, right) was more distinct and accurate (the other three examples of human gut shown in Fig. S4.). Secondly, for the same msmRNA-seq dataset, the number of species identified using the cell-by-HCK matrix was significantly higher than that identified using the cell-by-gene matrix, increasing by more than 4.5 times. At the same time, the number of each species in both the soil and the gut also increased significantly. This is particularly evident for five species with low or moderate abundance (number\u0026thinsp;\u0026lt;\u0026thinsp;300), including \u003cem\u003eBordetella pertussis\u003c/em\u003e, \u003cem\u003eFlavobacterium sp. CJ75\u003c/em\u003e, and \u003cem\u003eLabrys sp. KNU-23\u003c/em\u003e in the soybean soil sample. The newly identified strains accounted for 5/21 of the original strains. In the human gut samples, nine types of microorganisms were newly identified, including \u003cem\u003eAgathobacter rectalis\u003c/em\u003e, \u003cem\u003eFaecalibacterium prausnitzii\u003c/em\u003e, \u003cem\u003eParabacteroides distasonis\u003c/em\u003e, Phascolarctobacterium faecium, and \u003cem\u003eSutterella wadsworthensis\u003c/em\u003e. The number of newly identified bacterial species in the human gut (9 species) even exceeds the original number of species (7 species). Therefore, mKmer shows significant advantages in species identification in both complex soil environments and gut environments. This advantage is particularly pronounced in datasets where the original results were not very good and the species diversity was relatively low. Numerous studies have shown that the five bacterial species only identified by mKmer are commonly found in soil. \u003cem\u003eSphingopyxis terrae\u003c/em\u003e and \u003cem\u003eF. sp. CJ75\u003c/em\u003e are capable of degrading complex organic compounds [\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]. \u003cem\u003eMucilaginibacter mallensis\u003c/em\u003e can produce mucilaginous polysaccharides, thereby improving soil structure and fertility [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]. \u003cem\u003eL. sp. KNU-23\u003c/em\u003e is widely present in organic matter-rich soils [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. Surprisingly, \u003cem\u003eB. pertussis\u003c/em\u003e, primarily known as a human pathogen transmitted through the air, is not commonly found in soil environments, and its survival in soil has been seldom studied (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/nrmicro886\u003c/span\u003e\u003cspan address=\"10.1038/nrmicro886\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). The microbiome of the human gut is being studied more thoroughly. The species identified by mKmer, such as \u003cem\u003eA. rectalis\u003c/em\u003e, \u003cem\u003eF. prausnitzii\u003c/em\u003e, \u003cem\u003eMediterraneibacter gnavus\u003c/em\u003e, \u003cem\u003eOdoribacter splanchnicus\u003c/em\u003e, \u003cem\u003eP. distasonis\u003c/em\u003e, \u003cem\u003eP. faecium\u003c/em\u003e, \u003cem\u003ePhocaeicola coprophilus\u003c/em\u003e, \u003cem\u003eRoseburia intestinalis\u003c/em\u003e, and \u003cem\u003eS. wadsworthensis\u003c/em\u003e, are common human gut bacteria based on the literature [\u003cspan additionalcitationids=\"CR17 CR18 CR19 CR20 CR21 CR22 CR23\" citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. To further validate the reliability of our results, we applied the mKmer functions \u003cem\u003eKmerGOn\u003c/em\u003e and \u003cem\u003eKmerGOp\u003c/em\u003e to in-depth explore the functions of the mKmer identified species. The \u003cem\u003eB. pertussis\u003c/em\u003e in soybean soil exhibits a unique ability to bind iron ions in soybean soil (Fig. S3A, left). Iron is an essential micronutrient for plant growth, affecting soybean health and yield. Microorganisms can inhibit the growth of pathogens by competitively adsorbing iron in the soil, thereby reducing the occurrence of diseases. Compared to other microorganisms, acid phosphatase activity is significantly enriched in \u003cem\u003eBordetella pertussis\u003c/em\u003e (Fig. S3A, right), which can help decompose organic phosphorus compounds in the soil, releasing inorganic phosphorus that plants can absorb, thereby promoting phosphorus uptake by soybeans and enhancing crop growth. For human gut msmRNA-seq data, the \u003cem\u003eA. rectalis\u003c/em\u003e identified by mKmer revealed processes related to the metabolism of acids, including aconitate hydratase activity and the dicarboxylic acid metabolic process using \u003cem\u003eKmerGOn\u003c/em\u003e (Fig. S3B, left), consistent with findings by Abdugheni et al. [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. In addition, several processes related to the biosynthesis of amines have been found by \u003cem\u003eKmerGOp\u003c/em\u003e, such as 6-pyruvoyltetrahydropterin and tetrahydrobiopterin biosynthesis (Fig. S3B, right). Taken together, mKmer doesn\u0026rsquo;t depend on reference genomes, and can unbiasedly and effectively parse the complex biological information in msmRNA-seq data.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eA case study using mKmer\u003c/h3\u003e\n\u003cp\u003eTo demonstrate the practical applications of mKmer, we collected fecal samples from a colorectal cancer patient before and after immunotherapy for single-microbe sequencing. Using mKmer to analyze the two msmRNA-seq datasets, we identified 19 microbial species in the pre-treatment fecal sample (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA) and 27 species in the post-treatment sample (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eB, and the cell-by-gene matrix clustering results of these sample data are shown in Fig. S5). The microbial richness increased by over one-third, which is conducive to the restoration of a healthy gut microbiota [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. Through the analysis of these upregulated marker \u003cem\u003eK\u003c/em\u003e-mers with \u003cem\u003eK\u003c/em\u003e-motifs, we found that hormonal regulatory activity and carbohydrate metabolic processes were enriched in the post-treatment sample (Fig. S6). Further investigation into the shared species between the two samples, such as \u003cem\u003ePhocaeicola dorei\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eC), revealed that its populations in the pre- and post-treatment samples did not cluster together completely, indicating significant differences in gene expression. Consequently, we performed an in-depth analysis on functional changes in \u003cem\u003eP. dorei\u003c/em\u003e between these two samples. GO enrichment results for marker \u003cem\u003eK\u003c/em\u003e-mers in \u003cem\u003eP. dorei\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eD) showed that functions related to polysaccharide metabolism and outer membrane-associated defense responses were significantly upregulated in the post-treatment sample. Studies have shown that short-chain fatty acids produced by polysaccharide metabolism suggest efficient therapy and good prognosis for colorectal cancer [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e]. In addition, outer membrane binding and periplasmic space can promote biofilm formation, enhancing the host\u0026rsquo;s immune defense [\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. Therefore, the innovative mKmer method can help researchers more comprehensively analyze the dynamic changes of intestinal microecology during immunotherapy, and identify potentially beneficial microorganisms or their metabolites, which is expected to improve the immunotherapy efficacy of solid tumors.\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eOur study presents mKmer, a novel reference genome-free approach for analyzing msmRNA-seq data. The use of \u003cem\u003eK\u003c/em\u003e-mer for species taxonomic identification has been demonstrated by many classical \u003cem\u003eK\u003c/em\u003e-mer-based metagenomic taxonomy annotation software [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e, \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]. Although these tools classify species based on DNA, many conserved gene sequences remain highly consistent within species when DNA is transcripted into RNA. This implies that RNA sequences also contain species-specific conserved regions, and are useful in \u003cem\u003eK\u003c/em\u003e-mer analysis for species identification. For instance, rRNA and tRNA are widely used in taxonomic studies due to their significant conservation and variation among species [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e]. mRNA, on the other hand, reflects gene expression, which varies significantly between species. By analyzing high-frequency \u003cem\u003eK\u003c/em\u003e-mers in mRNA, species-specific expression characteristics can be captured. In this study, we discovered the presence of HCKs in every single-cell sequencing dataset, and further used them as genic sequences for downstream msmRNA-seq analysis. By leveraging the strong correlation between marker \u003cem\u003eK\u003c/em\u003e-mers and \u003cem\u003eK\u003c/em\u003e-motifs, we further explored and obtained reliable results on the functions of the microorganisms.\u003c/p\u003e \u003cp\u003eCompared to well-known tools such as Cell Ranger and STAR [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e], mKmer overcomes the limitations of incomplete or poorly annotated reference genomes by utilizing a cell-by-HCK matrix instead of the traditional cell-by-gene matrix. Benchmark tests with soybean soil and human gut msmRNA-seq data demonstrated that mKmer captures more data than those available tools and achieves clearer species clustering. This is understandable, as both Cell Ranger and STAR align the obtained msmRNA-seq data to the available reference genomes which have been sequenced. This process inevitably introduces biases. Considering the rapid evolution and variation of microorganisms, using single reference genomes per species does not align with established facts; alignment failures due to genetic variations would result in discarding a significant amount of valuable biological information obtained from the msmRNA-seq.\u0026nbsp;In contrast, mKmer, which operates without the need for reference genomes, provides a more comprehensive and unbiased analysis of microbiome samples.\u003c/p\u003e \u003cp\u003eIt is well known that scRNA-seq data contain numerous empty droplets and some doublets, as well as other impurity-containing droplets, which can severely impact the quality of sequencing results. Therefore, removing impurity information is crucial for single-cell analysis techniques. By examining the distribution of UMIs and barcodes, high-quality cells can be effectively filtered. Additionally, aligning to a reference genome to create a cell-by-gene matrix is an effective approach to exclude impurities from sequencing results. We know that the longer the \u003cem\u003eK\u003c/em\u003e-mer, the higher its specificity; thus, longer \u003cem\u003eK\u003c/em\u003e-mers are more efficient in detecting impurities. Conversely, shorter \u003cem\u003eK\u003c/em\u003e-mers have higher conservation, which increases information utilization when identifying the same species. mKmer extracts biological information from raw sequencing data by testing different \u003cem\u003eK\u003c/em\u003e-mer lengths to determine the optimal K. This optimal K can identify key conserved sequences while effectively distinguishing noise from non-critical gene sequences. As \u003cem\u003eK\u003c/em\u003e-mer gets longer,, the sequencing results from each sample showed a consistent trend (e.g. Figure\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). Specifically, at a certain \u003cem\u003eK\u003c/em\u003e, the number of \u003cem\u003eK\u003c/em\u003e-mers with a frequency of 1 was the highest among all \u003cem\u003eK\u003c/em\u003e-mer frequencies. We interpret \u003cem\u003eK\u003c/em\u003e-mers with a frequency of 1 as sequences that are useless for species clustering and may even be impurities. To effectively filter out useless sequence, we ranked each \u003cem\u003eK\u003c/em\u003e-mer in descending order of frequency and observed a distinct inflection point (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA and Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC). Based on this, we consider the high-frequency \u003cem\u003eK\u003c/em\u003e-mers before the inflection point as conservative \u003cem\u003eK\u003c/em\u003e-mers (i.e., HCKs), and cells containing these \u003cem\u003eK\u003c/em\u003e-mers are likely derived from a common ancestor (i.e., LCA) [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]. Therefore, HCKs exhibit high species recognition. The low-frequency \u003cem\u003eK\u003c/em\u003e-mers after the inflection point may indicate that the taxonomic information contained in these \u003cem\u003eK\u003c/em\u003e-mers is sparse, and they may even interfere with species identification and differentiation. It is not recommended to use low-frequency \u003cem\u003eK\u003c/em\u003e-mers in the downstream dimensionality reduction and clustering process. It is undeniable that there may be errors near the inflection point, where some highly specific conserved \u003cem\u003eK\u003c/em\u003e-mers may be misclassified as impurities. However, this has minimal impact on species identification across the entire cell set. In our analytical framework, we select these high-frequency conserved \u003cem\u003eK\u003c/em\u003e-mers at inflection points for species identification and functional analysis. At the same time, selection of HCKs address issues of computational time and high matrix dimensionality caused by the excessive variety of \u003cem\u003eK\u003c/em\u003e-mers, significantly enhancing mKmer's performance and practicality. Additionally, the high frequency of conserved \u003cem\u003eK\u003c/em\u003e-mers means that they are more likely to come from regions that can translate conserved protein motifs. These protein motifs have been used for cross-species functional annotation of single-cell RNA sequencing (scRNA-seq) [\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]. We further hypothesize that these HCKs may correspond to motif fragments of certain gene families within microorganisms. Therefore, by functionally annotating these gene motifs, mapping them onto the microbial communities, and performing enrichment analysis, we can infer the specific functions of the species in the sample.The key innovations of mKmer, including HCKs, marker \u003cem\u003eK\u003c/em\u003e-mers, and \u003cem\u003eK\u003c/em\u003e-motifs, enhance species identification and distinction. This method allows for a holistic view of microbial communities, advancing our understanding of microbial ecology and functional roles. Despite these strengths, mKmer still has room for improvement. For example, species annotation of Kraken 2 [\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e] can be corrected based on the results of clustering. Future research and development will further refine mKmer, for example, by expanding its functional modules, and optimizing its performance to provide more powerful support for microbiology research.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003emKmer is a reference genome-free approach for msmRNA-seq analysis and allows studies on cellular heterogeneity, marker motif discovery, and efficiency functional annotation. In this study, we discovered and defined HCK. More accurate clustering results confirmed that HCKs densely encapsulate a large amount of hierarchical information from microbial taxonomy. In benchmark tests on soybean soil and human gut msmRNA-seq datasets, mKmer can use more msmRNA-seq data than the traditional annotated gene-based methods, achieving more and clearer species clustering for subsequent comprehensive functional analysis. Our method therefore provides an unbiased way to analyze all species in msmRNA-seq samples and allows diverse microbiomic single-cell problems to be formulated in a unified way.\u003c/p\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003cdiv id=\"Sec13\" class=\"Section3\"\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003cdiv id=\"Sec16\" class=\"Section3\"\u003e \u003c/div\u003e \u003c/div\u003e "},{"header":"Methods","content":"\u003ch2\u003emsmRNA-seq data collection and generation\u003c/h2\u003e\n\u003cp\u003eA total of seven msmRNA-seq datasets, four from healthy donors by our previous study [\u003cspan class=\"CitationRef\"\u003e1\u003c/span\u003e] (PRJCA017256), two from a patient and one from soybean soil generated by this study, were used for performance and benchmarking of mKmer. These two fecal samples were collected from the same patient with colorectal cancer both before and after immunotherapy. The study protocol was approved by the Ethics Committee of the First Affiliated Hospital, Zhejiang University School of Medicine, China (2021IIT A0239). The protocols for sample treatment for high-throughput msmRNA-seq followed our previous study [\u003cspan class=\"CitationRef\"\u003e1\u003c/span\u003e] and the msmRNA-seq data generated by M20 Genomics were deposited at the NGDC database (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://ngdc.cncb.ac.cn/\u003c/span\u003e\u003c/span\u003e) under accession code PRJCAXXX (publicly available as of the date of publication).\u003c/p\u003e\n\u003ch2\u003eQuality control of raw data\u003c/h2\u003e\n\u003cp\u003eUMI-tools (v1.1.4) [\u003cspan class=\"CitationRef\"\u003e33\u003c/span\u003e] were used to process the unique molecular identifiers (UMIs) of our msmRNA-seq data. We utilized the \u003cem\u003eumi-tools whitelist\u003c/em\u003e for quality control of the raw data to estimate the number of cells accurately. Given that the raw data file R1 is approximately 1GB, we specified an expected cell number before determining the actual count. Therefore, the \u003cem\u003e--expect-cells\u003c/em\u003e parameter was set to 10,000. The raw data had a barcode length of 20bp and a UMI length of 8bp. To obtain the whitelist, the \u003cem\u003e--bc-pattern\u003c/em\u003e was set to \u003cem\u003eCCCCCCCCCCCCCCCCCCCC NNNNNNNN\u003c/em\u003e, and the \u003cem\u003e--set-cell-number\u003c/em\u003e parameter was set to 7,000, which corresponded to the cell number at the inflection point in the barcode rank plot. We then used the \u003cem\u003eumi-tools extract\u003c/em\u003e to filter the raw data files R1 and R2 based on the obtained whitelist, resulting in cleaned raw data.\u003c/p\u003e\n\u003cp\u003eDuring the PCR process of smRNA-seq, some molecules may have been disproportionately amplified due to sequence characteristics (e.g., GC content) or random factors, resulting in multiple reads with the same barcode and UMI. To address this, we employed the \u003cem\u003eRemoveDuplicates\u003c/em\u003e within UMI-tools to retain only the highest quality read, as determined by the Phred quality scoring system, among those with the same barcode and UMI. The sequencing data distinguishes each read in the form of 20bp_8bp, so the \u003cem\u003eRemoveDuplicates\u003c/em\u003e defaults to the last 29 characters of the sequence information line in the FASTQ file as a unique identifier.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eK\u003c/strong\u003e \u003cstrong\u003e-mers scanning and counting\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eJellyfish (v2.2.10) [\u003cspan class=\"CitationRef\"\u003e34\u003c/span\u003e] was used for fast, memory-efficient counting of \u003cem\u003eK\u003c/em\u003e-mers in DNA sequences. In our experiment, the \u003cem\u003e--m\u003c/em\u003e parameter of the \u003cem\u003ejellyfish count\u003c/em\u003e was set to the default value of 10, to determine the \u003cem\u003eK\u003c/em\u003e-mer length (\u003cem\u003eK\u003c/em\u003e). This \u003cem\u003eK\u003c/em\u003e corresponds to the smallest \u003cem\u003eK\u003c/em\u003e where the peak in \u003cem\u003eK\u003c/em\u003e-mer frequency distribution occurs at x\u0026thinsp;=\u0026thinsp;1. The \u003cem\u003e--s\u003c/em\u003e parameter was set to the default value of 10M. To ensure detection of all high-frequency \u003cem\u003eK\u003c/em\u003e-mers, the \u003cem\u003e--h\u003c/em\u003e parameter was set to 100,000,000 based on experimental testing.\u003c/p\u003e\n\u003cp\u003eTo confirm the selected \u003cem\u003eK\u003c/em\u003e value was reasonable, we observed the distribution of peaks with different \u003cem\u003eK\u003c/em\u003e values by drawing KmerFrequency plots. The \u003cem\u003e--put\u003c/em\u003e of the \u003cem\u003eKmerFrequency\u003c/em\u003e was three histo files with different \u003cem\u003eK\u003c/em\u003e values specified, and the \u003cem\u003e--out\u003c/em\u003e argument specified the path to the output KmerFrequency plot file. The histo file generated with the selected \u003cem\u003eK\u003c/em\u003e was used as the input for the \u003cem\u003eKmerRank\u003c/em\u003e to create a \u003cem\u003eK\u003c/em\u003e-mer rank plot, where the x value at the inflection point indicates the number of HCKs. The \u003cem\u003ejellyfish dump\u003c/em\u003e was then used to convert the \u003cem\u003ejf\u003c/em\u003e file into a readable format for extracting the top-counts \u003cem\u003eK\u003c/em\u003e-mers.\u003c/p\u003e\n\u003cp\u003eHCK\u0026thinsp;=\u0026thinsp;select_top {sort [count (K-mer, C)]}\u003c/p\u003e\n\u003cp\u003eM\u003csub\u003eij\u003c/sub\u003e = count (HCK\u003csub\u003ei\u003c/sub\u003e, C\u003csub\u003ej\u003c/sub\u003e)\u003c/p\u003e\n\u003cp\u003eIn the frequency matrix \u003cem\u003eM\u003c/em\u003e, a specific element \u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003eij\u003c/em\u003e\u003c/sub\u003e represents the occurrence count of the \u003cem\u003ei\u003c/em\u003e-th HCK in the \u003cem\u003ej\u003c/em\u003e-th cell among the selected HCKs.\u003c/p\u003e\n\u003ch2\u003eGeneration of cell-by-HCKs matrix\u003c/h2\u003e\n\u003cp\u003eAfter detecting each \u003cem\u003eK\u003c/em\u003e-mer, they were sorted by detection depth in descending order. HCKs were then selected based on this sorted list. Using the generated HCK list, \u003cem\u003eK\u003c/em\u003e-mers were read sequentially from the cleaned R2 reads. Each \u003cem\u003eK\u003c/em\u003e-mer in every cell was counted to generate the cell-by-HCKs matrix. To minimize memory usage during execution, the program generated cell-by-HCKs matrix for every 1,000 cells read and then merged these matrices. For ease of use, we integrated the entire process of cell-by-HCKs matrix into a single command named \u003cem\u003eKmerCell\u003c/em\u003e. The \u003cem\u003e--kmercount\u003c/em\u003e argument required the \u003cem\u003ekmer_counts_dumps.fa\u003c/em\u003e file output from jellyfish dump (the file suffix must be \u003cem\u003e_counts_dumps.fa\u003c/em\u003e); \u003cem\u003e--fastq\u003c/em\u003e required the clean R2 FASTQ file; \u003cem\u003e--topkmer\u003c/em\u003e specified the number of HCKs; and \u003cem\u003e--k\u003c/em\u003e specified the selected \u003cem\u003eK\u003c/em\u003e.\u003c/p\u003e\n\u003ch2\u003eIdentification of microbial species\u003c/h2\u003e\n\u003cp\u003eWe employed a \u003cem\u003eK\u003c/em\u003e-mer-based root-to-leaf classification strategy for microbial species annotation, which is integrated in our software under the name \u003cem\u003esmAnnotation\u003c/em\u003e. We first applied Kraken 2 (v 2.0.7-beta) [\u003cspan class=\"CitationRef\"\u003e32\u003c/span\u003e], a \u003cem\u003eK\u003c/em\u003e-mer-based read classification method, on every read in each barcode based on standard refseq of Kraken 2(Refseq archaea, bacteria, viral, plasmid, human1, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://benlangmead.github.io/aws-indexes/k2\u003c/span\u003e\u003c/span\u003e). After all the reads were assigned into each node of different taxonomic levels (e.g., order, family, genus, species), we calculate the sum of reads in each node from the leaf to the root. Then we performed the taxonomic classification from the root to leaf taxonomic levels. In the root taxonomic level, we ranked all the nodes from the highest to lowest based on the number of reads of the nodes, and selected the node with the most reads as a potential annotation candidate. Based on the annotation results, then we performed the same annotation process in the next lower taxonomic level, until the leaf nodes (species level). Then Bracken (v 2.5) [\u003cspan class=\"CitationRef\"\u003e35\u003c/span\u003e] was used to count the \u003cem\u003efraction_total_reads\u003c/em\u003e of the species classified into each cell, and the species with the largest value was selected as the final annotation result. For Kraken 2, the comparison database was the NCBI standard database by default (Archive size: 60GB) and resolution parameter \u003cem\u003e--r\u003c/em\u003e was set to 100. For \u003cem\u003esmAnnotation\u003c/em\u003e, clean R2 as the specified file for \u003cem\u003e--input\u003c/em\u003e, and the output file named \u003cem\u003esmAnnotation.report\u003c/em\u003e was placed in the current working directory by default.\u003c/p\u003e\n\u003ch2\u003eVisualization and clustering\u003c/h2\u003e\n\u003cp\u003eTo visualize the data, we further reduced the dimensionality of all filtered cells using Seurat (v4) [\u003cspan class=\"CitationRef\"\u003e36\u003c/span\u003e] and used UMAP to project the cells into 2D space. The annotation results of Kraken 2 and Bracken were mapped to the Seurat object; only the annotation results with \u003cem\u003efraction_total_reads\u003c/em\u003e values greater than 0.5 were retained, and strains with abundance less than 0.1% were filtered out. The steps include: (i) Using the LogNormalize method of the \u003cem\u003eNormalizeData\u003c/em\u003e of Seurat to calculate the expression values of \u003cem\u003eK\u003c/em\u003e-mers. The scale.factor argument is set to the default 10000, nfeatures to 6000, and the \u003cem\u003eScaleData\u003c/em\u003e object to all genes; (ii) PCA was performed using the normalized expression value; among all the principal components, the top 30 principal components were used to do clustering and UMAP analysis; (iii) To find clusters, a weighted graph-based clustering method, Shared Nearest Neighbour (SNN), was selected, and the resolution is set to 0.5. Marker genes for each cluster were identified with the MAESTRO test with default parameters via the \u003cem\u003eFindAllMarkers\u003c/em\u003e in Seurat and the min.pct parameter was set to 0.25.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGO annotation by marker\u003c/strong\u003e \u003cstrong\u003eK\u003c/strong\u003e\u003cstrong\u003e-mers and\u003c/strong\u003e \u003cstrong\u003eK\u003c/strong\u003e\u003cstrong\u003e-motifs\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eBefore performing functional enrichment analysis on sequencing data from different samples, we first need to filter out rRNA from the raw transcriptome data. In this study, we used SortMeRNA (v4.2.0) [\u003cspan class=\"CitationRef\"\u003e37\u003c/span\u003e] for filtering, with the reference dataset including multiple species, such as bacteria, eukaryotes, archaea, and various sequencing databases, for 5S, 16S, 23S, and other rRNA types. After using the \u003cem\u003eFindAllMarkers\u003c/em\u003e, we obtained a list of marker \u003cem\u003eK\u003c/em\u003e-mers, which were used for functional enrichment analysis of clusters or species of interest. Marker \u003cem\u003eK\u003c/em\u003e-mers are considered to be identified from highly conserved sequences, which are likely to represent individual motifs. The functional analysis of bacterial species based on their specific motifs is reliable. The MEME (Multiple Em for Motif Elicitation) suite (v5.0.5) [\u003cspan class=\"CitationRef\"\u003e38\u003c/span\u003e] is a comprehensive resource for discovering and analyzing sequence motifs in DNA, RNA, and protein sequences. Memes [\u003cspan class=\"CitationRef\"\u003e39\u003c/span\u003e] is an R package that provides a seamless R interface to a selection of popular MEME Suite tools. By analyzing the conserved sequences of each strain, we aimed to elucidate the specific functions of the strains. To obtain the specific functions of bacterial species, we designed two partitioning workflows to conduct Gene Ontology (GO) enrichment analysis on \u003cem\u003eK\u003c/em\u003e-motifs.\u003c/p\u003e\n\u003ch2\u003eNucleotide motif analysis (KmerGOn)\u003c/h2\u003e\n\u003cp\u003eThe first analysis workflow involved converting each marker \u003cem\u003eK\u003c/em\u003e-mer into a motif file in MEME format. These motifs were then compared against motifs in the microbial nucleotide motif database using the \u003cem\u003etomtom\u003c/em\u003e [\u003cspan class=\"CitationRef\"\u003e40\u003c/span\u003e] integrated into the MEME suite. This step identified the best-matching known motifs. Subsequently, \u003cem\u003eama\u003c/em\u003e and \u003cem\u003egomo\u003c/em\u003e [\u003cspan class=\"CitationRef\"\u003e41\u003c/span\u003e] in the MEME suite were used to compare the identified known motifs against the Escherichia coli database, obtaining GO functions for each motif. Finally, GO functional enrichment and visualization were performed on the clusters or species of interest.\u003c/p\u003e\n\u003ch2\u003eProtein motif analysis (KmerGOp)\u003c/h2\u003e\n\u003cp\u003eThe second analysis workflow utilized SeqKit (v2.8.2) [\u003cspan class=\"CitationRef\"\u003e42\u003c/span\u003e, \u003cspan class=\"CitationRef\"\u003e43\u003c/span\u003e] to translate each marker \u003cem\u003eK\u003c/em\u003e-mer into amino acids using six reading frames (since \u003cem\u003eK\u003c/em\u003e\u0026thinsp;=\u0026thinsp;12, only sequences with 4 AA were retained). These motifs were then compared against motifs in the all-species motif database using the \u003cem\u003etomtom\u003c/em\u003e integrated into the MEME suite. Each protein motif was further analyzed using InterProScan (v5.47-82.0) [\u003cspan class=\"CitationRef\"\u003e44\u003c/span\u003e] to search domain databases (e.g., Pfam, PROSITE, PRINTS, etc.) and obtain GO IDs. Finally, used \u003cem\u003eselect\u003c/em\u003e of the AnnotationDbi (v1.64.1) [\u003cspan class=\"CitationRef\"\u003e45\u003c/span\u003e] package to match the corresponding term and GO ID from the GO.db (v3.18.0) [\u003cspan class=\"CitationRef\"\u003e46\u003c/span\u003e] database.\u003c/p\u003e\n\u003cp\u003eFor both workflows, the output list of marker \u003cem\u003eK\u003c/em\u003e-mers from the \u003cem\u003eFindAllMarkers\u003c/em\u003e served as the input file. The --\u003cem\u003ecluster\u003c/em\u003e was set to the target cluster, and the --\u003cem\u003eout\u003c/em\u003e specified the path to the output file. This ensured a systematic approach to uncovering the functional roles of conserved sequences within bacterial species.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eSupplementary Information\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe online version contains supplementary material available at XXX.\u003c/p\u003e\n\u003cp\u003eAdditional file 1: Includes supplementary figures Fig S1\u0026mdash;Fig S6. A short description of the content of these figures is provided at the first page of the file.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePeer review information\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eXXX was the primary editor of this article at Genome Biology and managed its editorial process and peer review in collaboration with the rest of the editorial team.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eReview history\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe review history is available as XXX.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026rsquo; contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eLF and WJ conceived the study. WJ, WC, YW, WC and XZ conducted the experiments. FM, QQ, XL, DZ, JY, HC, YH and SW analyzed the data. FM and LF developed mKmer. FM, LF and YB wrote the paper. LF, WJ, YW and YS discussed and supervised this project. All authors have revised and approved the final manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis study was supported by Biological Breeding-Major (2023ZD04076), Yunnan Tobacco Company (2024530000241001) and CIC-MIC.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll msmRNA-seq data used by this study are available at NGDC database (https://ngdc.cncb.ac.cn/) under project number PRJCA017256 (accession number SAMC3766839, SAMC3766838, SAMC3766837, SAMC1266599), PRJCAXXX (the soil sample) and PRJCAXXX (the two gut samples) (publicly available as of the date of publication). The mKmer package (v.1.0.0) is available at https://github.com/bioinplant/mKmer. Any additional information required to reanalyze the data reported in this work paper is available from the lead contact upon request.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll data collection was approved by the Ethics Committee of the First Affiliated Hospital, Zhejiang University School of Medicine, China (2021IIT A0239).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare no competing interests\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eShen Y, Qian Q, Ding L, Qu W, Zhang T, Song M et al (2024) High-throughput single-microbe RNA sequencing reveals adaptive state heterogeneity and host-phage activity associations in human gut microbiome. Protein Cell. ;pwae027\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMacosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M et al (2015) Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell 161(5):1202\u0026ndash;1214\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S et al (2013) STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29(1):15\u0026ndash;21\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eEren AM, Delmont TO (2024) Bioprospecting marine microbial genomes to improve biotechnology. Nature 633(8029):287\u0026ndash;288\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen H, Ryu J, Vinyard ME, Lerer A, Pinello L (2024) SIMBA: single-cell embedding along with features. Nat Methods 21(6):1003\u0026ndash;1013\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTayyebi Z, Pine AR, Leslie CS (2024) Scalable and unbiased sequence-informed embedding of single-cell ATAC-seq data with CellSpace. Nat Methods 21(6):1014\u0026ndash;1022\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim J, Steinegger M (2024) Metabuli: sensitive and specific metagenomic classification via joint analysis of amino acid and DNA. Nat Methods 21(6):971\u0026ndash;973\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKwok AYC, Su S-C, Reynolds RP, Bay SJ, Av-Gay Y, Dovichi NJ et al (1999) Species identification and phylogenetic relationships based on partial HSP60 gene sequences within the genus Staphylococcus. Int J Syst Evol MicroBiol 49:1181\u0026ndash;1192\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWood DE, Salzberg SL (2014) Kraken: ultrafast metagenomic sequence classification using exact alignments. Genome Biol 15(3):R46\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim D, Song L, Breitwieser FP, Salzberg SL (2016) Centrifuge: rapid and sensitive classification of metagenomic sequences. Genome Res 26(12):1721\u0026ndash;1729\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMenzel P, Ng KL, Krogh A (2016) Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat Commun 7:11257\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSharma M, Khurana H, Singh DN, Negi RK (2021) The genus Sphingopyxis: Systematics, ecology, and bioremediation potential - A review. J Environ Manage 280:111744\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim Y-S, Hwang E-M, Jeong C-M, Cha C-J (2023) Flavobacterium psychrotrophum sp. nov. and Flavobacterium panacagri sp. nov., Isolated from Freshwater and Soil. J Microbiol 61(10):891\u0026ndash;901\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eM\u0026auml;nnist\u0026ouml; MK, Tiirola M, McConnell J, H\u0026auml;ggblom MM (2010) Mucilaginibacter frigoritolerans sp. nov., Mucilaginibacter lappiensis sp. nov. and Mucilaginibacter mallensis sp. nov., isolated from soil and lichen samples. Int J Syst Evol MicroBiol 60(Pt 12):2849\u0026ndash;2856\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim M-J, Park Y-J, Park M-K, Park CE, Jo Y, Tagele SB et al (2020) Complete genome sequence of Labrys sp. KNU-23 isolated from ginseng soil in the Republic of Korea. Korean J Microbiol 56(4):410\u0026ndash;412\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAbdugheni R, Wang W, Wang Y, Du M, Liu F, Zhou N et al (2022) Metabolite profiling of human-originated Lachnospiraceae at the strain level. iMeta 1(4):e58\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRamirez-Farias C, Slezak K, Fuller Z, Duncan A, Holtrop G, Louis P (2008) Effect of inulin on the human gut microbiota: stimulation of Bifidobacterium adolescentis and Faecalibacterium prausnitzii. Br J Nutr 101(4):541\u0026ndash;550\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eIchimura R, Tanaka K, Nakato G, Fukuda S, Arakawa K (2024) Complete genome sequence of Mediterraneibacter gnavus strain RI1, isolated from human feces. Microbiol Resour Announc 0:e00863\u0026ndash;e00824\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi H, Xu H, Li Y, Jiang Y, Hu Y, Liu T et al (2020) Alterations of gut microbiota contribute to the progression of unruptured intracranial aneurysms. Nat Commun 11(1):3218\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang K, Liao M, Zhou N, Bao L, Ma K, Zheng Z et al (2019) Parabacteroides distasonis Alleviates Obesity and Metabolic Dysfunctions via Production of Succinate and Secondary Bile Acids. Cell Rep 26(1):222\u0026ndash;235e5\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003e1, Wu F, Guo X, Zhang J, Zhang M, Ou Z, Peng Y (2017) Phascolarctobacterium faecium abundant colonization in human gastrointestinal tract. Experimental Therapeutic Med 14(4):3122\u0026ndash;3126\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGarc\u0026iacute;a-L\u0026oacute;pez M, Meier-Kolthoff JP, Tindall BJ, Gronow S, Woyke T, Kyrpides NC et al (2019) Analysis of 1,000 Type-Strain Genomes Improves Taxonomic Classification of Bacteroidetes. Front Microbiol 10:2083\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLeth ML, Ejby M, Workman C, Ewald DA, Pedersen SS, Sternberg C et al (2018) Differential bacterial capture and transport preferences facilitate co-growth on dietary xylan in the human gut. Nat Microbiol 3(5):570\u0026ndash;580\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWilliams BL, Hornig M, Parekh T, Lipkin WI (2012) Application of Novel PCR-Based Methods for Detection, Quantitation, and Phylogenetic Characterization of Sutterella Species in Intestinal Biopsy Samples from Children with Autism and Gastrointestinal Disturbances. mBio 3(1):e00261\u0026ndash;e00211\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLozupone CA, Stombaugh JI, Gordon JI, Jansson JK, Knight R (2012) Diversity, stability and resilience of the human gut microbiota. Nature 489(7415):220\u0026ndash;230\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLouis P, Flint HJ (2017) Formation of propionate and butyrate by the human colonic microbiota. Environ Microbiol 19(1):29\u0026ndash;41\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCosterton JW, Stewart PS, Greenberg EP (1999) Bacterial Biofilms: A Common Cause of Persistent Infections. Science 284(5418):1318\u0026ndash;1322\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJanda JM, Abbott SL (2007) 16S rRNA Gene Sequencing for Bacterial Identification in the Diagnostic Laboratory: Pluses, Perils, and Pitfalls. J Clincal Microbiol 45(9):2761\u0026ndash;2764\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBermudez-Santana C, Attolini CS-O, Kirsten T, Engelhardt J, Prohaska SJ, Steigele S et al (2010) Genomic organization of eukaryotic tRNAs. BMC Genomics 11:270\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNasko DJ, Koren S, Phillippy AM, Treangen TJ (2018) RefSeq database growth influences the accuracy of k-mer-based lowest common ancestor species identification. Genome Biol 19(1):165\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRosen Y, Brbić M, Roohani Y, Swanson K, Li Z, Leskovec J (2024) Toward universal cell embeddings: integrating single-cell RNA-seq datasets across species with SATURN. Nat Methods 21(8):1492\u0026ndash;1500\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWood DE, Lu J, Langmead B (2019) Improved metagenomic analysis with Kraken 2. Genome Biol 20(1):257\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSmith T, Heger A, Sudbery I (2017) UMI-tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Res 27(3):491\u0026ndash;499\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMar\u0026ccedil;ais G, Kingsford C (2011) A fast, lock-free approach for efficient parallel counting of occurrences of \u003cem\u003ek\u003c/em\u003e -mers. Bioinformatics 27(6):764\u0026ndash;770\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLu J, Breitwieser FP, Thielen P, Salzberg SL (2017) Bracken: estimating species abundance in metagenomics data. PeerJ Comput Sci 3:e104\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A et al (2021) Integrated analysis of multimodal single-cell data. Cell 184(13):3573\u0026ndash;3587e29\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKopylova E, No\u0026eacute; L, Touzet H (2012) SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics 28(24):3211\u0026ndash;3217\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBailey TL, Johnson J, Grant CE, Noble WS (2015) The MEME Suite. Nucleic Acids Res 43(W1):W39\u0026ndash;49\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNystrom SL, McKay DJ, Memes (2021) A motif analysis environment in R using tools from the MEME Suite. Pertea M, editor. PLOS Computational Biology. ;17(9):e1008991\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGupta S, Stamatoyannopoulos JA, Bailey TL, Noble W (2007) Quantifying similarity between motifs. Genome Biol 8(2):R24\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBuske FA, Bod\u0026eacute;n M, Bauer DC, Bailey TL (2010) Assigning roles to DNA regulatory motifs using comparative genomics. Bioinformatics 26(7):860\u0026ndash;866\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShen W, Le S, Li Y, Hu F, SeqKit: (2016) A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation. Zou Q, editor. PLOS ONE. ;11(10):e0163962\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShen W, Sipos B, Zhao L (2024) SeqKit2: A Swiss army knife for sequence and alignment processing. iMeta 3(3):e191\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJones P, Binns D, Chang H-Y, Fraser M, Li W, McAnulla C et al (2014) InterProScan 5: genome-scale protein function classification. Bioinformatics 30(9):1236\u0026ndash;1240\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePag\u0026egrave;s H, Carlson M, Falcon S, Nian L, AnnotationDbi (2023) Manipulation of SQLite-based annotations in Bioconductor. Bioconductor. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.18129/B9.bioc.AnnotationDbi\u003c/span\u003e\u003cspan address=\"10.18129/B9.bioc.AnnotationDbi\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eThe Gene Ontology Consortium, Carbon S, Douglass E, Good BM, Unni DR, Harris NL et al (2021) The Gene Ontology resource: enriching a GOld mine. Nucleic Acids Res 49(D1):D325\u0026ndash;D334\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[{"identity":"4298afbf-b796-4aaa-bb86-24e5603b6f71","identifier":"10.13039/100006206","name":"Biological and Environmental Research","awardNumber":"2023ZD04076","order_by":0},{"identity":"18bedcfa-d8a4-47a7-a1fd-a4f703e5d06b","identifier":"10.13039/501100010825","name":"China Tobacco Yunnan Industrial Corp","awardNumber":"2024530000241001","order_by":1}],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"Hainan Institute of Zhejiang University","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"K-mer, msm RNA-seq, HCK, reference genome-free, K-motif","lastPublishedDoi":"10.21203/rs.3.rs-5748035/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-5748035/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe advanced single-microbe RNA sequencing (smRNA-seq) technique addresses the pressing need to understand the complexity and diversity of microbial communities, as well as the distinct microbial states defined by different gene expression profiles. Current analyses of smRNA-seq data heavily rely on the integrity of reference genomes within the queried microbiota. However, establishing a comprehensive collection of microbial reference genomes or gene sets remains a significant challenge for most real-world microbial ecosystems. Here, we developed an unbiased embedding algorithm utilizing \u003cem\u003eK\u003c/em\u003e-mer signatures, named mKmer, which bypasses gene or genome alignment to enable species identification for individual microbes and downstream functional enrichment analysis. By substituting gene features in the canonical cell-by-gene matrix with highly conserved \u003cem\u003eK\u003c/em\u003e-mers, we demonstrate that mKmer outperforms gene-based methods in clustering and motif inference tasks using benchmark datasets from crop soil and human gut microbiomes. Our method provides a reference genome-free analytical framework for advancing smRNA-seq studies.\u003c/p\u003e","manuscriptTitle":"mKmer: An unbiased K-mer embedding of microbiomic single-microbe RNA sequencing data","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-01-03 07:16:39","doi":"10.21203/rs.3.rs-5748035/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"429d26a3-e268-44b5-8a46-80a52b4c6181","owner":[],"postedDate":"January 3rd, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":42242116,"name":"Bioinformatics"}],"tags":[],"updatedAt":"2025-01-03T07:16:39+00:00","versionOfRecord":[],"versionCreatedAt":"2025-01-03 07:16:39","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-5748035","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-5748035","identity":"rs-5748035","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}