Low compositions of human toll-like receptor 7/8-stimulating RNA motifs in the MERS-CoV, SARS-CoV and SARS-CoV-2 genomes imply a substantial ability to evade human innate immunity

Data collection

The complete (full-length) genome sequences of 1,002 coronaviruses infecting humans (including 22 human coronavirus 229E (HCoV-229E), 61 human coronavirus NL63 (HCoV-NL63), 26 human coronavirus HKU1 (HCoV-HKU1), 139 human coronavirus OC43 (HCoV-OC43), 63 severe acute respiratory syndrome coronavirus (SARS-CoV), 121 Middle East respiratory syndrome coronavirus (MERS-CoV), 432 severe acute respiratory SARS-CoV-2) and 138 unclassified coronaviruses (Other) were retrieved from the Virus Pathogen Resource (ViPR, https://www.viprbrc.org/) (Pickett et al., 2012) and analyzed in this study (Table S1). In addition, 1,762 complete (full-length) genomic sequences of coronaviruses from six nonhuman mammalians (bat, camel, cat, civet, dog and pig) and avian hosts were retrieved from the Virus Pathogen Resource and used for analysis (Table S2).

Ninety-five oligoribonucleotides (ORNs) and 39 ribonucleotide tetramers with experimentally validated human TLR 7/8-stimulating activity were identified from 17 research reports (Table S3). The sequences of TLR 7 proteins from different organisms exhibit high variations. For example, the sequence identity of TLR 7 proteins from humans and mice is 81%. The preferences for ligand nucleotide compositions of TLRs from different organisms might be different. Since the experiments validating the TLR 7/8-stimulating activity of these ORN sequences were conducted using human cells, the TLR-stimulating triribonucleotide composition and TLR-stimulating scores described in this study should be considered to be specific for human TLR 7/8.

Weighted triribonucleotide compositions of single-stranded RNA virus genomes

A method to evaluate the human TLR 7/8-stimulating ability of single-stranded RNA virus genomes based on their human TLR 7/8-stimulating triribonucleotide compositions was developed by Yang & Chen (2012). The 4³ = 64 possible trimers are labeled as X₁, X₂, …, X₆₄. Each trimer frequency f_Xi is defined as (1) $$f_{\mathrm{X}\mathrm{i}}=c_{\mathrm{X}\mathrm{i}}/l,i=1,2,\dots ,64$$here c_Xi is the number of trimer Xi and l is the total number of trimers. The trimer weights w_Xi were computed using the following formula: (2) $$w_{Xi}={\log }_{10}\left({\displaystyle \frac{f_{Xi}}{1/64}}\right)=\{\begin{array}{c}{w}_{Xi}^{+}\mathrm{i}\mathrm{f}{f}_{Xi}>{\displaystyle \frac{1}{64}}\\ w_{Xi}^{-}\mathrm{i}\mathrm{f}{f}_{Xi}<{\displaystyle \frac{1}{64}}\end{array}$$where $w_{Xi}^{+}and{w}_{Xi}^{-}$ are the weights of overrepresented and underrepresented human TLR 7/8-stimulating trimers, respectively. If the relative frequency of a trimer in the human TLR 7/8-stimulating ORN sequences is greater than 1/64 (the expected value of a random distribution), the trimer is considered to be human TLR 7/8 stimulatory. Otherwise, the trimer is considered to be nonhuman TLR 7/8-stimulatory. Each trimer is assigned a weight based on the logarithm of its relative frequency in the human TLR 7/8-stimulating ORN sequences (Fig. 1).

For any individual RNA virus genome, the positive and negative weighted trimer compositions were calculated and are referred to as Score S and Score N, respectively, these scores are collectively referred to as the human TLR 7/8-stimulating scores. Score S for stimulating trimers was calculated as (3) $$\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}S=\Sigma (c_{Xi}{w}_{Xi}^{+})/l$$and Score N for nonstimulating trimers was calculated as (4) $$\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}N=\Sigma (c_{Xi}{w}_{Xi}^{-})/l$$where c_Xi is the number of trimer Xi that appear in the viral genomic RNA (with i = 1, … 64). l is the total number of trimers in the viral genomic RNA. Higher Score S and lower Score N values indicates greater numbers of human TLR 7/8-stimulating triribonucleotides in the viral RNA genome and implies that stronger human TLR 7/8-mediated innate immunity may be induced by this viral RNA. Conversely, lower Score S and higher Score N values indicate greater numbers of human TLR 7/8 nonstimulating triribonucleotides in the viral RNA genome and, as a consequence, a higher potential for evasion of the human TLR 7/8-mediated innate immune responses (Fig. 1).

Data analysis

Data manipulation was performed with Perl scripts written by the author. The heatmap, stripchart of Logit P values and scatter plots of Score S and Score N values were plotted using the ggplot2 package of R (the R package for statistical computing). Logistic regression was conducted with the glm function of R. Linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA) were performed with the lda and qda functions, respectively, in the MASS package of R. Naive Bayes and support vector machine (svmLinear, svmPoly and svmRadial) classifiers in the caret package of R were used.

Tenfold cross-validation

The cross-validation test is well-known and has been used in many computation-based studies (Le, Ho & Ou, 2017; Le et al., 2019; Le, Yapp & Yeh, 2019). The weighted triribonucleotide compositions of coronavirus genomes were randomly split into 10 subsets. Ten pairs of training (90% data for model building) and test (10% data for model evaluation) data were combined from the 10 subsets and used to perform 10-fold cross validations. For each test run, a training set was used to train a model, and the model was then tested using the test set. All methods (logistic regression, LDA, QDA, naive Bayes, svmLinear, svmPoly and svmRadial) were used to perform 10-fold cross validations. Sensitivity was computed with the following formula: sensitivity = true positive/(true positive + false negative). Specificity was computed with the following formula: specificity = true negative/(true negative + false positive). Accuracy was computed with the following formula: accuracy = (true positive + true negative)/total number of samples.

Analysis of coronaviruses from nonhuman animals

Seven models derived from the seven methods (logistic regression, LDA, QDA, naive Bayes, svmLinear, svmPoly and svmRadial) were used to analyze data of coronaviruses from nonhuman animals. All seven methods were selected for binary classification to distinguish human coronaviruses causing common colds and severe acute respiratory syndromes. Using the results of logistic regression as a standard, the overall agreements between the results of the logistic regression model and those of the other 6 models were computed. The overall agreements were computed by the following formula: (true positive + true negative)/total number of samples.

Compositions of triribonucleotides in genomes of coronaviruses infecting humans

The similarity of triribonucleotide compositions was not consistent with the similarity of genome sequences (phylogenetic analysis in Fig. 2A). The overall triribonucleotide compositions of genomic (plus) strand and complementary (minus) strand RNAs of coronaviruses infecting humans are shown in Fig. 2B. These results indicate that triribonucleotide compositions provide novel information that cannot be revealed by phylogenetic analysis.

Human TLR 7/8-stimulating potential of human coronavirus genomes

The human TLR 7/8-stimulating scores of coronavirus genomic (positive strand) RNAs are shown in Fig. 3. The human TLR 7/8-stimulating potential of coronavirus genomic (positive strand) RNAs followed the order of NL63-CoV > HKU1-CoV >229E-CoV ≅ OC63-CoV > SARS-CoV-2 > MERS-CoV > SARS-CoV. The human TLR 7/8-stimulating scores of the complementary (negative) strand of coronavirus genomic RNAs (replication intermediates) are shown in Fig. 4. The human TLR 7/8-stimulating potential of the complementary (negative) strand of coronavirus genomic RNAs followed the order of SARS-CoV-2 > SARS-CoV >229E-CoV > MERS-CoV > NL63-CoV ≅ OC63-CoV > HKU1-CoV. The highest Score S value of the complementary (negative) strand of SARS-CoV-2 genomic RNA wass 0.0813 which was less than the minimum Score S value of SARS-CoV genomic (positive strand) RNA. Greater numbers of human TLR 7/8-stimulating triribonucleotides were found in the genomic (positive strand) RNAs than in the complementary strand RNAs of all coronaviruses analyzed. These results suggest a greater potential for human TLR 7/8-stimulating activity in the early stage of viral infection than the complementary (negative) strand RNAs produced during coronavirus replication stage.

A logistic regression model for predicting the human TLR 7/8-stimulating potential

To predict the human TLR 7/8-stimulating potential of coronaviruses, a logistic regression model was constructed as follows: (5) $$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}(p)=\log (p/1-p)={\mathrm{\beta }}_0+{\mathrm{\beta }}_1\cdot S+{\mathrm{\beta }}_2\cdot N+{\mathrm{\beta }}_3\cdot \mathrm{S}\mathrm{r}+{\mathrm{\beta }}_4\cdot \mathrm{N}\mathrm{r}$$S and N are the Score S and Score N values for the positive strand of a viral genome sequence. Sr and Nr are the Score S and Score N values for the negative strand of a viral genome sequence. 229E-CoV, OC43-CoV, NL63-CoV and HKU1-CoV were used as the highly stimulating group. MERS-CoV, SARS-CoV and SARS-CoV-2 were used as the poorly stimulating group. After the model selection procedure, the following model was selected as the final model: (6) $$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}(p)=\log (p/1-p)={\mathrm{\beta }}_0+{\mathrm{\beta }}_1\cdot S$$The results of 10-fold cross-validation are shown in Table 1. The averages of the intercepts and coefficients of the 10 models were used to construct the logistic regression model as follows: (7) $$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}(p)=\log (p/1-p)=-71.01+730.49\cdot S$$The logistic regression model using all data is: (8) $$\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{i}\mathrm{t}(p)=\log (p/1-p)=-69.05+710.36\cdot S$$

Comparison with other methods

Six methods were used to validate the results of logistic regression. The results of 10-fold cross-validations using linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), naive Bayes and three support vector machine classifiers (linear, polynomial and radial) were the same as those using the logistic regression model (Table 1). The high values (≅1) of sensitivity, specificity and accuracy may be due to the almost complete separation of the Score S values of the two groups. These results are consistent with the data shown in Fig. 3.