Single-cell simultaneous metabolome and transcriptome profiling revealing metabolite-gene correlation networ

Metabolic studies at single cell level could directly define the cellular phenotype closest to physiological or disease states. However, the current single cell metabolome (SCM) study using mass spectroscopy has the difficulty to give a complete view of the metabolic activity in the cell, while the prediction of metabolism-phenotype relationship is limited by the potential inconsistency between transcriptomic and metabolic levels. Here, single-cell simultaneous metabolome and transcriptome profiling method (scMeT-seq) is developed, based on sub-picoliter sampling from the cell for the initial metabolome profiling followed by single cell transcriptome sequencing. This design does not only provide sufficient cytoplasm for SCM, but also nicely keeps the cellular viability for the accurate transcriptomic analysis in the same cell. Diverse relationships between the two omics are revealed in macrophages with the stimulation of lactate and lactate transporter (MCT1) antagonist. Moreover, we mapped metabolite states to the single-cell differentiation trajectory and gene correlation network of macrophages, which allows the unsupervised functional interpretation of metabolome. Thus, the established scMeT-seq should lead to a new perspective in metabolic research by transforming metabolomics from metabolite snapshot to functional approach. The scRNA-seq library was constructed following a modified Smart-seq2 protocol(21). In detail, 4ul lysed sample was incubated at 72C for 3 min. The first-strand cDNA was reverse-synthesized using oligo(dT) and template-switching oligonucleotides. The synthesized first-strand cDNAs were amplified by 20 cycles. After purification, 0.1 ng of cDNA was used for Nextera tagmentation (Illumina) and library construction. Sequencing was then performed on Illumina HiSeq 4000.Single-cell samples were consolidated and subjected to sequencing on the Illumina Novaseq 6000 platform, targeting an average of 3 million reads per library in a paired-end 150-bp format. The raw sequencing data underwent demultiplexing and were subsequently processed with the fastp software (version 0.20.0) for quality trimming. The quality-trimmed reads were then aligned to the mouse reference genome mm10 using the HISAT2 aligner (version 2.1.0) with standard settings. The counts of reads per gene were quantified utilizing the StringTie package. The processed data matrix was import into Seurat (version 4.0). We took the transcriptome matrix as input and then scaled it before using dimensionality reduction techniques. Next, we performed PCA on the scaled data and selected the top principal components to preserve biological variation and remove technical noise. We added the quality-controlled metabolic matrices to the new assay and performed data standardization and normalization. We used a graph-based clustering method to cluster macrophages. We selected gene sets in GSEA (https://www.gsea-msigdb.org/gsea/index.jsp) and assessed the pathway activity of the cells using the GSVA method. We combined the pathway matrix and the normalized metabolome matrix, calculated the correlation between pathways and metabolites using the spearman method, and visualized the results by scatter plot in ggplot2. Based on the grouping information, Wilcoxon rank sum test was performed using FindMarkers function. Genes with |LogFC|>0.25 and P<0.05 were considered as signature genes, and the results were visualized by volcano plots of ggplot2.The Monocle2 package was used to analyze macrophage trajectories. We imported the transcriptome and metabolite matrices into Monocle and used the feature genes and metabolites to sort the pseudotemporal trajectories of the cells. The "DDRTree" was applied to reduce the space to two dimensions and the trajectories in the reduced dimensional space were visualized using the function "plot_cell_trajectory". Cells were divided into three states based on the trajectories, and Gene Ontology (GO) enrichment analysis was performed for the differentially expressed genes of each state using the ClusterProfiler package. The GENIE3 (https://bioconductor.org/packages/release/bioc/html/GENIE3.html) package was used to analyze the metabolite-gene covariance network. We imported the merged matrix of bi-omics and performed functional enrichment analysis for the genes in the three selected modules using the above software and visualized them by bar plot of ggplot2. The schematic was drawn by the BioRender Graphics Tool (http://BioRender.com).