The mRNA-seq raw sequence data in BAM format of liver cancer patients’ normal and tumor tissue biopsies for this study were downloaded from the Seven Bridges Cancer Genomics Cloud (CGC https://www.cancergenomicscloud.org/). The data set included 98 samples (49 pairs of RNA-Seq sequencing samples of primary liver cancer tumor tissue and adjacent normal tissue). The read sequences that did not map to the reference genome GRCH38 were extracted from the BAM files. These foreign reads were then mapped to the RefSeq (NCBI Reference Sequence Database) (https://www.ncbi.nlm.nih.gov/refseq/) - a large collection of bacterial genomes. The mapping results were then input to the GRAMMy pipeline for the relative abundance estimation of infiltrating bacteria. The obtained expression and microbial profiles were used for downstream analysis. The data processing flow is shown in Figure 1.

We first filtered the reads that were not matched to the human genome reference GRCH38 from the RNA-Seq samples in BAM format with SAMtools, which included reads that were not mapped at both ends and reads that were not mapped at one end (mapped at one end but not at the other). Apply Fastp to remove low-quality reads and partial reads, and excise poor quality bases. Use the BWA tool to remap the unmapped reads to a comprehensive microbial reference library RefSeq, containing the full genome sequence of 23,790 bacteria and archaea. We then applied GRAMMy algorithm to estimate the bacteria relative abundance.

The mRNA analysis pipeline begins with the Alignment Workflow, which is performed using a two-pass method with STAR. STAR aligns each read group separately and then merges the resulting alignments into one. Following alignment, BAM files are processed through the RNA Expression Workflow to determine RNA expression levels. The reads mapped to each gene are enumerated using HT-Seq-Count. Expression values are provided in a tab-delimited format.

Three real data sets were used to validate the PSO-SVR algorithm as compared to the CIBERSORTX, EPIC, and MCP-counter algorithms. The first validation data set from the CIBERSORT article consists of the gene expression profiles of 20 peripheral blood mono nuclear cells (PBMC) samples and the ratio of immune cells determined through flow cytometry for the same samples. The second data set comes from the National Center for Biotechnology Information Search database (NCBI database) under the GEO Datasets GSE64385. It consists of the gene expression profile of 10 colorectal cancer (CRC) samples and the composition of immune cells as determined by immunohistochemistry of the same samples. The third data set from CIBERSORTX consists of the gene expression profiles of 19 melanoma samples and the composition of immune cells obtained through single-cell sequencing technique of the same samples. The full description of these data sets is shown in Supplementary Table 1.

We first filtered the reads that did not matched to the human genome reference GRCH38 from the RNA-Seq samples in BAM format using SAMtools (27). The filtered reads included those that were not mapped at both ends and those that were not mapped at one end (mapped at one end but not the other). The process was as follows:

Extracting unmapped single reads: we used the command “samtools view -u -f 4 -F 264” to extract single-end unmapped reads in the format of FASTQ (28).

Extracting unmapped paired-end reads: we used the command “samtools view -u -f 12 -F 256” to extract paired-end unmapped reads with the format FASTQ.

Sequencing quality issues like low-confidence bases and sequence-specific bias complicate mRNA-seq analyses (29). We applied comprehensive quality control (QC) before analyses (30). We applied Fastp (31) to remove low-quality reads and partial reads, and to excise poor quality bases.

For the extracted unmapped paired-end reads, we used the command “fastp -q 0 -u 100 -n 10 -l 36 -A -G -M 0 -i” to delete sequences with base quality lower than 40% of Q15, sequences with N greater than 5, sequences with length less than 36, and broke up sequences, respectively. For the extracted unmapped single reads, we used the command “fastp -q 0 -u 100 -n 10 -l 36 -A -G -M 0 -i” to perform the same QC process on single-ended sequenced sequences. To merge the extracted single-end and paired-end unmapped reads, we used the SMS2 (32) software to add a reverse complementary sequence to the unmapped single reads. We then converted the resulting FASTQ formatted reads into the FASTA (33) format through SeqKit (34).

We then used the BWA tool (35) to remap the unmapped reads to a comprehensive microbial reference library RefSeq, containing the full genome sequence of 23,790 bacteria and archaea. We then applied the GRAMMy algorithm to estimate the relative abundance of bacteria to mitigate the problem of ambiguous mapping of short reads to relative reference sequences.

The Mann Whitney U test, also known as “Mann Whitney rank sum test”, was proposed by H.B.mann and D.Rwhitney in 1947. It assumes that any two samples are from two populations that are exactly the same except the population mean, in order to test for significant difference between the mean of the two populations.

Assuming that the mean values of two populations exist, they are recorded as μ₁,μ₂ respectively. With only one translation difference between f₁ and f₂ at most, we get: μ₁ = μ₂ – α. The assumptions to be tested are as follows:

The steps of Mann Whitney U test are:

In the test, because the critical value table of Mann Whitney U test only gives a smaller critical value, the smaller U value in U_A and U_B is used as the test statistic.

4. Select the smaller U value to compare with the critical value of U. if U > U_α (α = 0.05, accept the original assumption H₀. If U < U_α (α = 0.05, Then reject H₀ and accept H₁. The acceptance domain is the same as Wilcoxon test. U test can also be divided into small samples and large samples. In case of small samples, the critical values of U have been compiled into a table. In large samples, the distribution of U tends to be normal, so it can be treated by normal approximation.

The loss function is an index to measure the performance of the prediction model in predicting the expected results. The commonly used loss functions are mean squared error (MSE) and root mean squared error (RMSE). Since MSE squares the error (y – y^predicted = e), if e > 1, the value of the error will increase a lot. If there is an outlier in our data, the value of e will be very high and will be much greater than │e│. This will make the model with MSE loss give a higher weight to outliers. In order to minimize this outlier data point, we use the RMSE value, namely root mean square error, but at the expense of the prediction effect of other normal data points, which will eventually reduce the overall performance of the model.

Correlation analysis refers to the analysis of two or more variable elements with correlation, so as to measure the correlation degree of two variable factors Correlation analysis can be carried out only when there is a certain connection or probability between the elements of correlation.

(1) Pearson correlation coefficient

Given two continuous variables x and y, the Pearson correlation coefficient is defined as:

Where x ¯ and y ¯ are the mean values of the variables x and y respectively.

ρ close to 0 indicates that there is no correlation between the two variables; whereas close 1 or - 1 indicates that the two variables are strongly correlated.

(2) Spearman correlation coefficient

Spearman correlation coefficient is defined as Pearson correlation coefficient ρ between hierarchical variables. For samples with a sample size of N, N original data are converted into hierarchical data. Compared with Pearson correlation coefficient, Spearman correlation coefficient is insensitive to data errors and extreme values, which is defined as.

Where R and S are the grades of observed values i respectively, R ¯ and S ¯ are the average grades of variables x and y respectively, and N is the total number of observed values.

Logistic regression, also known as log probability regression, is a machine learning method used to solve the binary classification problem, which is used to estimate the possibility of something. It does not need to scale the input features, and the interpretability of the model is very good. The influence of different features on the final result can be seen from the weight of features. We fit the following regularization model to binary features as:

Where β is the regression coefficient. Parameter α is the balanced lasso (L1) and ridge (L2) regularization, and λ determines their weights.

Therefore, we use the logistic regression classifier to classify 98 samples. Their corresponding bacterial relative abundance data and immune cell proportion data were used as input characteristics, and the sample status (normal/tumor) was used as the prediction variable 0 or 1. 75% of the data were used for training and 25% for testing. The association between hepatocarcinogenesis and tumor bacteria and invasive immune cells was analyzed through the classification results under different input characteristics.

Support Vector Regression (SVR) is a commonly used technique for non-linear unsupervised learning. The generalization ability and prediction accuracy of SVR depend on the choice of kernel function coefficients and penalty factors. Due to the large feature size, small sample size and unknown sample distribution, in order to avoid over fitting as much as possible, we used a simple and effective linear kernel function in order to avoid overfitting as much as possible. We introduced the particle swarm optimization (PSO) technique, which uses a powerful algorithm to invoke an iterative approach to solve parameter optimization. It can improve the accuracy of SVR results, and increase the robustness of the results through multiple iterations. The overall flowchart is shown in Figure 2.

We briefly introduce the SVR model for immune cell infiltration estimation. Based on a large amount of tissue gene expression profile data combined with a priori knowledge of purified leukocyte subpopulation expression profile data (i.e., a “signature matrix” representing the expression profile data of each gene in different immune cells), the proportion of immune cells in tumor biopsies can be accurately estimated. The idea was first used in CIBERSORT. The gene expression profile data were first transformed into a linear combination of marker genes for immune cells to solve for f (representing the proportion of immune cells) in a linear combination equation. The linear combination equation was expressed as:

Where X_n×L denotes the expression profile data of n genes from L samples for deconvolution. f_m×L is the proportion of immune cells to be sought, specifically indicating the proportion of each immune cell in L samples among the m immune cell types. S_n×m is the signature matrix, indicating the expression profile data of n genes in each immune cell for m immune cells. The SVR algorithm solves the linear combinatorial equation using a deconvolution scheme.

The accuracy of the immune cell ratio calculated through the above SVR model is influenced by the model parameters such as the penalty factor C, sensitivity ϵ, and kernel function coefficient ϕ. These parameters are crucial to an SVR model’s accuracy and generalizability. However, in practice, these parameters are largely manually chosen without justification. Here we propose to optimize (ϵ, C, ϕ) parameters through iterations using the particle swarm algorithm. In each iteration, the swarming particles update their velocity and position vectors by tracking two “optimal values” (p_ibest, g_best), where p_ibest denotes the individual optimal value of the particle and g_best is the global optimal value. In the computation, the expression data of each gene (i.e., each row vector in the X_{n ×L} matrix) was treated as a particle in the swarming iteration.

Where v_i and x_i represent the velocity vector and position vector of the i^th particle respectively, each gene is regarded as a particle respectively, and n represents the size of the population, specifically the number of genes; q is a non-negative inertia factor. The larger the value of q is, the stronger the global optimization ability is and the weaker the local optimization ability is; c₁ and c₂ are learning factors, general c₁ = c₁ = 2; r₁ and r₁ both represent random coefficients belonging to [0,1]; p_ibest represents the individual optimal value of the ith particle, and g_best is the global optimal value. We provide more detailed parameter explanation and algorithm flow in the Supplementary Materials. Table 1 provides an explanation of key parameters.

We apply this iterative optimization process of particle swarm to solve the SVR model, so that its optimal input parameters are optimized according to the swarm fitness metric. Upon initialization, we specify the approximate ranges of the parameters ϵ, C, and ϕ as: ϵ = (0, 0.2), C = (1, 100), ϕ = (0.01, 2.0), and q_max and q_min values are 0.9 and 0.4 for the model. Upon convergence, we find the optimal parameter solution (ϵ, C, ϕ). SVR models using optimal parameters are validated using benchmark data sets. We apply model with these parameters to the expression data of liver cancer patients’ tissues to solve for immune cell infiltration ratios.

The major innovation of PSO-SVR over previous SVR deconvolution methods is that it applies the machine learning technique – particle swarm algorithm, to do SVR iterative optimization. In the process, the swarming algorithm searches for SVR hyperplane formed by support vectors that capture as many data points as possible while satisfying the given constraints and avoiding overfitting, using a linear “ω-insensitive” loss function. The function penalizes only those data points outside a specific error radius. The support vectors were selected from the signature matrix of genes, and the standard reference expression profile (signature matrix) was selected from a composition of 22 immune cells including CD4 T, and CD8 T immune cells in CIBERSORTX, but other signature gene sets can be applied as well.

We gathered three data sets (see Supplementary Table 1), in which the immune cell proportions were determined through orthogonal flow cytometry (PBMC-FC), immunohistochemistry (CRC-IC), and single-cell RNA sequencing technologies (Melanoma-scRNA). Flow cytometry was an established technology to identify and determine different cell types in heterogeneous cell populations (36). immunohistochemistry is also well-established, using chemical reaction to label antibodies and to identify and quantify antigens within tissue cells. The genetic heterogeneity of cells of the same tissue can be also analyzed by scRNA-Seq at the level of individual cells to cluster and compute the composition immune cells (37). In order to verify the accuracy of PSO-SVR algorithm, we used data sets based on these technologies as the orthogonal control, and studied the deviation between PSO-SVR results and the control.

As shown in Figure 3A, for the PBMC-FC data, that the PSO-SVR estimated B cell, CD4 T cells, CD8 T cells and Monocytes levels were very close to the flow cytometry results. Overall, the four types of immune cells present a good correlation. The CD8 T cells showed relatively less consistency as compared to the other two types, maybe because CD8 T cells accounted for a relatively small number of T lymphocytes in PBMC (5%–20%) thus is subject to less accurate estimation. As shown in Figure 3B, the PSO-SVR estimates of B cells, monocytes, NK cells, and T cells were also in good concordance with the immunohistochemistry results. As shown in Figure 3C, PSO-SVR estimated immune cell fractions, such as CD4 T and CD8 T cells, also showed a good correlation with the scRNA-seq levels.

To further evaluate the PSO-SVR estimates using the PBMC-FC, CRC-IC and Melanoma-scRNA data sets, we used three metrics, namely root mean square error (RMSE), Pearson correlation coefficient, and Spearman’s correlation coefficient. The individual indicators are shown in Table 2.

We observe in Table 2 that the error was generally larger for B cells and that both Pearson and Spearman’s correlation coefficients were relatively low. For CD4 T and CD8 T cells, the difference between the results presented in the PBMC-FC and Melanoma-scRNA is very small, and the RMSE and correlation coefficients for both are close, proving the reliability of the algorithm in this study.

Finally, the results calculated through PSO-SVR were compared with other algorithms such as CIBERSORTX, EPIC, MCP-counter, for calculating the proportion of immune cells in infiltrating tumors. From the results of PBMC-FC and Melanoma-scRNA (refer to Supplementary Tables 2 and 5), it can be observed that PSO-SVR outperforms CIBERSORTX, EPIC, MCP-counter, in terms of the three indicators, RMSE, Pearson, and Spearman correlation coefficients. However, the results of CRC-IC (refer to Supplementary Table 3) are not as good as the other three methods, which could be attributed to the small amount of CRC-IC data. Overall, our algorithm was highly accurate.

In order to further analyze the impact of tumor infiltrating immune cells and bacteria on the occurrence and development of liver cancer, we used the downloaded gene expression data and transcriptome data of 98 liver cancer samples, and calculated the corresponding bacterial abundance data by using the estimation method of microbial relative abundance in this paper. The PSO-SVR algorithm described in Section 3.1 was used to calculate the relative proportion of immune cells in liver cancer samples.

Then the logistic regression classification method was applied to three cases: case 1, only bacterial abundance data were used as the input feature; In case 2, only the immune cell proportion data were used as the input feature; and in case 3, both bacterial abundance and immune cell ratio were used as input characteristics. From the classification results of the three cases (see Table 4), the classification accuracy of using only a single input feature is lower - only bacteria: 0.70, only immune cell feature: 0.74, and both features: 0.84. No matter which data feature is added, the classification accuracy can be improved. After adding the immune cell proportion data, the classification accuracy is improved by 20% compared with case 1; After adding bacterial abundance data, the accuracy of classifier is improved by 13.51% compared with case 2. This shows that the bacteria and infiltrating immune cells in liver cancer tissue contain key characteristics for cancer diagnosis, which can be used as a reference index for clinical cancer diagnosis.