Differences in the coding region of the canonical transcripts

Compared to general profile features (e.g., number of ORFs), the sequences themselves provide more direct mapping to the protein function and structure [38]. Here, we encoded 344 discrete features and 7,026 categorical features from complementary DNA (cDNA) of the canonical transcript to explore features specific to ISGs. We divided the discrete features into 4 categories (nucleotide composition/dinucleotide composition/codon usage/nucleotide 4-mer composition) and compared their representations among 3 different groups of human genes, including recompiled ISGs from dataset S2, recompiled non-ISGs from dataset S2, and the background human genes from dataset S1 (Fig. 5).

Open in a separate window

Figure 5:

Differences in the representation of discrete features encoded from coding regions (canonical). Mann–Whitney U tests are applied for hypothesis testing on the whole comparing data without sampling, and the results are provided in Supplementary Data S2. Here, the ISGs and non-ISGs are taken from dataset S2 (No. = 620 and 874) while the background human genes are from dataset S1 (No. = 10,836) (Table 5).

First, guanine and cytosine were both more depleted in ISGs than non-ISGs, leading to an underrepresentation of GC content in the ISGs (Mann–Whitney U test: M₁ = 52%, M₂ = 55%, P = 2.3E-11). This attribute is the opposite to the GC-biased gene conversion (gBGC) process and would result in ISGs being less stable with weak evolutionary conservation (Fig. 4) [39]. Additionally, the underrepresentation of GC content also influenced the representation of other dinucleotide features. Among all dinucleotide depletions in ISGs, CpG depletion was ranked first followed by GpG and GpC depletions (P = 2.9E-14, 4.9E-13, and 1.2E-10, respectively). In turn, adenine- and thymine-related dinucleotide composition, exemplified by ApT and TpA, were more enriched in ISGs than non-ISGs (P = 8.0E-10 and 8.5E-10, respectively).

We compared the usage of 64 different codons in the third category as their frequencies influence transcription efficiency [40]. Differences between the ISGs and background human genes were observed in codons for 11 amino acids, including leucine (L), isoleucine (I), valine (V), serine (S), threonine (T), alanine (A), glutamine (Q), lysine (K), glutamic acid (E), arginine (R), and glycine (G). The most significant difference was observed in the usage of codon “AGA.” Among all arginine-targeted alternative codons, codon “AGA” was usually favoured, and its presence reached an estimated 25% in the ISGs but reduced to 22% in the background human genes (P = 1.4E-05). It was significantly lower in the non-ISGs, at 18% (P = 1.9E-13). On the other hand, compared to the background human genes, the codon “CAG” coding for amino acid “Q” was the most underrepresented in the ISGs. It was less favoured by the ISGs than non-ISGs (M₁ = 72%, M₂ = 78%, P =7.3E-13), although it dominated in coding patterns. As for the 3 stop codons, compared with the background human genes, the usage of the TAA stop codon was overrepresented in the ISGs (M₁ = 28%, M₂ = 33%, P =9.7E-03). In this category of codon usage, the features with different frequencies between the ISGs and background human genes became more discriminating when comparing the ISGs with non-ISGs. Significant differences in codon usages between the ISGs and non-ISGs were widely observed except for methionine (M) and tryptophan (W). Hence, despite the limited differences of codon usages between the ISGs and background human genes, these features were useful for discriminating the ISGs from non-ISGs.

In the last category, we calculated the occurrence frequency of 256 nucleotide 4-mers to add some positional resolution for finding and comparing interesting organisational structures [41]. Among the 256 4-mers, 46 of them were differentially represented between the ISGs and background human genes (Supplementary Data S2). Most of these 4-mers were overrepresented by the ISGs except 2 with the pattern “TAAA” and “CGCG.” Interestingly, the feature of “TAAA” composition became a positive factor when comparing ISGs and non-ISGs (M₁ = 4.1%, M₂ = 3.7%, P =4.1E-06), suggesting it might be a suitable feature to discern potential or incorrectly labelled ISGs. We found that 6 nucleotide 4-mers (“ACCC,” “AGTC,” “AGTG,” “TGCT,” “GACC,” and “GTGC”) were overrepresented in the ISGs when compared to the background human genes. However, they were not differentially represented when comparing the ISGs with non-ISGs. These 6 features might be inherently biased for some reason and were not powerful enough to contribute to distinguishing the ISGs from non-ISGs. In addition to the aforementioned 40 features (except 4-mer “ACCC,” “AGTC,” “AGTG,” “TGCT,” “GACC,” and “GTGC”) that were differentially represented in ISGs compared to background human genes, we found a further 39 features that nucleotide 4-mers differentially represented between ISGs and non-ISGs (Supplementary Data S2).

To check the effect of these aforementioned 343 features on the level of stimulation in the IFN-α system (log₂(fold change) >0), we calculated the PCC for the normalised features (Eq. 2) and found 106 features were positively related to the increase of fold change, and 34 features were suppressed when human genes were more upregulated after IFN-α treatments (Student t-test: P < 0.05) (Supplementary Data S3). ApA composition showed the most obvious positive correlation with stimulation level (PCC = 0.464, P =8.8E-06), while a negative association between the representation of 4-mer “CGCG” and IFN-α–induced upregulation was the most significant (PCC = −0.593, P =3.2E-09). Human genes with higher upregulation in the presence of IFN-α contained more codons “CAA,” rather than “CAG” for coding amino acid “Q.” The depletion of GC content, especially cytosine content, promotes the suppression of many nucleotide compositions in the cDNA (e.g., CpG composition).

To find conserved sequence patterns relating to gene regulation [42], we checked the existence of 2,940, 44,100, and 661,500 short linear nucleotide patterns (SLNPs) consisting of 3 to 5 consecutive nucleobases in the group of the ISGs and non-ISGs. By using a positive 5% difference in the occurrence frequency as the cutoff threshold, we found 7,884 SLNPs with a maximum difference in representation of around 15%. After using Pearson's chi-squared tests and Benjamini–Hochberg correction to avoid type I error in multiple hypotheses [43], 7,025 SLNPs remained with an adjusted P-value lower than 0.01 (Supplementary Data S4), hereon referred to as “flagged” SLNPs. The differentially represented 7,025 SLNPs were ranked according to the adjusted P-value. As shown in Fig. 6A, dinucleotide “TpA” dominates in the top 10, top 100, top 1,000, and all differentially represented SLNPs even if TpA representation is suppressed in the cDNA of genes’ canonical transcripts compared to other dinucleotides. Dinucleotide “ApT” and “ApA” are also frequently observed in the flagged SLNPs, but their occurrences do not show significant differences in the top 100 SLNPs (Pearson's chi-squared test: P > 0.05). GC-related dinucleotides (e.g., “CpC,” “GpC,” and “GpG”) are rarely observed in the flagged SLNPs, especially in the top 10 or top 100. In view of this, we hypothesize that the differential representation of nucleotide compositions influences and reflects on the pattern of SLNPs in the ISGs. By checking the co-occurrence status of the flagged SLNPs, we found that these sequence patterns had a cumulative effect in distinguishing the ISGs from non-ISGs, especially when the number of co-occurring SLNPs reached around 5,320 (Pearson's chi-squared test: P =7.9E-13, Fig. 6B). There were 8 (~1.3%) ISGs in dataset S2 containing all the flagged 7,025 SLNPs. Their upregulations after IFN-α treatment were generally low with a fold change fluctuating around 2.2. However, some of these 8 genes, such as desmoplakin (DSP), were clearly highly upregulated in endothelial cells isolated from human umbilical cord veins after not only IFN-α treatments (fold change = 11.1) but also IFN-β treatments (fold change = 13.7). We also found some non-ISGs (e.g., hemicentin 1 [HMCN1]) and human genes with limited expression in the IFN-α experiments (ELGs) (e.g., tudor domain containing 6 [TDRD6]) containing the flagged SLNPs, but their frequencies were lower than that in the ISGs.

Open in a separate window

Figure 6:

Short linear nucleotide patterns (SLNPs) in the coding regions (canonical). (A) Influence of dinucleotide composition on the flagged SLNPs. (B) The co-occurrence status of SLNPs in different human genes. Ranks in (A) are generated based on the adjusted P-value given by Pearson's chi-squared tests after the Benjamini–Hochberg correction procedure. Detailed results of the hypothesis tests are provided in Supplementary Data S4. Here, the ISGs and non-ISGs are taken from dataset S2 while the background human genes are from dataset S1 (Table 5).

Differences in the protein amino acid sequence

We used the amino acid sequences generated by the canonical transcript to extract features at the proteomic level. In addition to the basic composition of 20 standard amino acids, we considered 17 additional features related to physicochemical (e.g., hydropathy and polarity) or geometric properties (e.g., volume) [44, 45]. We found several amino acids that were either enriched or depleted in the ISG products compared to the background human proteins, which were produced by genes in dataset S1 (Fig. 7). The differences were even more marked between protein products of the ISGs and non-ISGs, highlighting some differences that were not observed when comparing the ISG products to the background human proteins (e.g., isoleucine composition). The differences observed in the amino acid composition were at least in part associated with the patterns previously observed in features encoded from genetic coding regions. For example, asparagine (N) showed significant overrepresentation in the ISG products compared to the non-ISG products or background human proteins (Mann–Whitney U test: P =2.8E-12 and 1.2E-03, respectively). This was expected as there are only 2 codons (i.e., “AAT” and “AAC” coding for amino acid “N”), and dinucleotide “ApA” showed a remarkable enrichment in the coding region of ISGs. A similar explanation could be given for the relationship between the deficiency of GpG content and amino acid “G.” The translation of amino acid “K” was also influenced by ApA composition but was not significant due to the mild representation of dinucleotide “ApG” in the genetic coding region. Additionally, as previously mentioned, the ISGs showed a significant depletion in the CpG content, and consequently, the amino acids “A” and “R” in the ISG products were significantly underrepresented. Cysteine (C) was not frequently observed in human proteins but still showed a relatively significant enrichment in the ISG products (M₁ = 2.3%, M₂ = 2.5%, P =1.8E-03).

Open in a separate window

Figure 7:

Differences in the representation of discrete features encoded from protein sequences. Mann–Whitney U tests are applied for hypothesis testing on the whole data without sampling and the results are provided in Supplementary Data S2. Here, the ISGs and non-ISGs are taken from dataset S2 (No. = 620 and 874) while the background human genes are from dataset S1 (No. = 10,836) (Table 5). Aliphatic group: amino acids “A,” “G,” “I,” “L,” “P,” and “V”; aromatic/huge group: amino acids “F,” “W,” and “Y” (volume >180 cubic angstroms); sulphur group: amino acids “C” and “M”; hydroxyl group: amino acids “S” and “T”; acidic/negative_charged group: amino acids “D” and “E”; amide group: amino acids “N” and “Q”; positive_charged group: amino acids “R,” “H,” and “K”; hydrophobic group: amino acids “A,” “C,” “I,” “L,” “M,” “F,” “V,” and “W” that participate in the hydrophobic core of the structural domains [46]; neutral group: amino acids “G,” “H,” “P,” “S,” “T,” and “Y”; hydrophilic group: amino acids “R,” “N,” “D,” “Q,” “E,” and “K”; tiny group: amino acids “G,” “A,” and “S” (volume <90 cubic angstroms); small group: amino acids “N,” “D,” “C,” “P,” and “T” (volume ranged from 109 to 116 cubic angstroms); medium group: amino acids “Q,” “E,” “H,” and “V” (volume ranged within 138 to 153 cubic angstroms); large group: amino acids “R,” “I,” “L,” “K,” and “M” (volume ranged within 163 to 173 cubic angstroms); uncharged group: the remaining 15 amino acids except electrically charged ones; polar group: amino acids “R,” “H,” “K,” “D,” “E,” “N,” “Q,” “S,” “T,” and “Y”; nonpolar group: the remaining 10 amino acids except polar ones.

When focusing on the composition of amino acid sequences grouped by physicochemical or geometric properties, we found some features differentially represented between the ISG products and background human proteins. The result showed that hydroxyl (amino acids “S” and “T”), amide (amino acids “N” and “Q”), or sulphur amino acids (amino acids “C” and “M”) were more abundant in the ISG products compared to the background human proteins (Mann–Whitney U test: P =0.04, 1.0E-03, and 0.02, respectively). Small amino acids (amino acids “N,” “C,” and “T”; aspartic acid [D]; and proline [P]; the volume ranging from 108.5 to 116.1 cubic angstroms) were more frequently observed in the ISG products than in background human proteins (M₁ = 22.1%, M₂ = 21.7%, P =0.02). These differences became more marked when comparing the representation of these features between the ISG and non-ISG products. For example, features relating to chemical properties of the side chain (e.g., aliphatic), charge status, and geometric volume showed differences between proteins produced by the ISGs and non-ISGs. Some features such as neutral amino acids that include amino acids “G,” “P,” “S,” and “T”; histidine (H); and tyrosine (Y) were not differentially represented between the ISG and non-ISG products, but they indicated an obvious association with the change of IFN-α–triggered stimulations (PCC = −0.556, P =4.1E-08) (Supplementary Data S3).

Next, we searched the sequence of the ISG products against that of the non-ISG products to find conserved short linear amino acid patterns (SLAAPs), which might be constrained by strong purifying selection [47]. As opposed to the analysis of the genetic sequence, we obtained only 19 enriched sequence patterns with a Pearson's chi-squared P-value ranging from 1.5E-04 to 0.02 (Table 1), hereon referred to as flagged SLAAPs. They were greatly influenced by 4 polar amino acids, “K,” “N,” “E,” and “S,” and 1 nonpolar amino acid: “L.” Some of these flagged SLAAPs (e.g., SLAAP “NVT” and “S-N-E”) were clearly overrepresented in the ISG products compared to the background human proteins and could be used as features to differentiate the ISGs from background human genes. The third column in Table 1 indicates a number of patterns that are lacking in the non-ISG products and hence may be the reason for the lack of upregulation in the presence of IFN-α. Particularly, we noticed that SLAAP “KEN” was a destruction motif that could be recognised or targeted by anaphase promoting complex (APC) for polyubiquitination and proteasome-mediated degradation [48, 49]. Results shown in Fig. 8A illustrate that the co-occurrence of differentially represented SLAAPs (flagged) has a cumulative effect in distinguishing the ISGs from non-ISGs. This cumulative effect can even be achieved with only 2 random SLAAPs (Pearson's chi-squared test: P =4.6E-10). The bias in the co-occurring SLAAPs (flagged) in the background human proteins towards a pattern similar to the non-ISG products further proves the importance of these 19 SLAAPs. However, their co-occurrence is not associated with the level of IFN-triggered stimulations (PCC = 0.015, P > 0.05) (Fig. 8B).

Open in a separate window

Figure 8:

Representation of co-occurring short linear amino acid patterns (SLAAPs, flagged), in our main dataset. (A) The co-occurrence status of SLAAPs in different classes. (B) Relationship between co-occurrence of the marked SLAAPs and log₂(fold change) after IFN-α treatments. Here, the ISGs and non-ISGs are taken from dataset S2 while the background human genes are from dataset S1 (Table 5). Points in (B) are located based on the average feature representation of genes with similar expression performance in IFN-α experiments.

Regions that lacked stable structures under normal physiological conditions within proteins are termed intrinsically disordered regions (IDRs). They play an important role in cell signalling [50]. Compared with ordered regions, IDRs are usually more accessible and have multiple binding motifs, which can potentially bind to multiple partners [51]. According to the results calculated by IUPred [52], we identified 6,721, 10,510, and 119,071 IDRs (IUpred score no less than 0.5) in proteins produced by the ISGs, non-ISGs, and background human genes, respectively. We hypothesize that enriched SLAAPs widely detected in the IDRs may be important for human protein–protein interactions or potentially virus mimicry [53]. For instance, in the ISG products, about 40.8% of SLAAP “SxNxT” were observed in the IDRs, 14.9% higher than that in non-ISG products (Table 1). This difference reflected the importance of SLAAP “SxNxT” for target specificity of IFN-α–induced protein–protein interactions (PPIs) [9], even if it was not statistically significant. By contrast, the conditional frequency of SLAAP “SxNxE” in the IDRs of the ISG and non-ISG products was almost the same, indicating that SLAAP “SxNxE” might have an association with some inherent attributes of the ISGs but was less likely to be involved in the IFN-α–induced PPIs. SLAAP “KEN” in the IDRs also showed some interesting differences: in the non-ISG products, 41.9% of SLAAP “KEN” were observed in the IDRs, 14.6% higher than that in the ISG products, which provided an effective approach to distinguish the ISGs from non-ISGs. When SLAAP “KEN” is discovered in the ordered globular region of a protein sequence, statistically, the protein is more likely to be produced by an ISG, but this assumption is reversed if the SLAAP is located in an IDR (Pearson's chi-squared tests: P =0.03). Despite the relatively low conditional frequency of SLAAP “KEN” in the IDRs of the ISG products, these SLAAPs in the IDR are more likely to be functionally active than those falling within ordered globular regions [54].

Differences in network profiles

We constructed a network with 332,698 experimentally verified interactions among 17,603 human proteins (confidence score >0.63) from the Human Integrated Protein–Protein Interaction rEference (HIPPIE) database [55] to investigate if the connectivity among human proteins has an association with genes’ expression in the IFN-α experiments. In total, 10,169 out of 10,836 human proteins produced by genes in our background dataset S1 were included in the network. Based on this network, we calculated 8 features as defined in the methods, including the average shortest path, closeness, betweenness, stress, degree, neighbourhood connectivity, clustering coefficient, and topological coefficient.

As illustrated in Fig. 9B and G, ISG products tend to have higher values of betweenness and stress than background human proteins (Mann–Whitney U test: P =0.01 and 0.03, respectively), which means they are more likely to locate at key paths connecting different nodes of the PPI network. Some ISG products with high values of betweenness and stress (e.g., tripartite motif containing 25 [TRIM25]) can be considered the shortcut or bottleneck of the network and play important roles in many PPIs, including those related to the IFN-α–triggered immune activities [56, 57]. However, such differential representation of betweenness does not mean ISG products are more likely to be or even be close to bottlenecks of the network compared to the background human proteins. Some examples shown in Table 2 indicate that ISG products are less connected by top-ranked bottlenecks and hubs of the network than non-ISG products or the background human proteins. This conclusion is not influenced by the hub/bottleneck protein's performance in the IFN-α experiments. Comparing proteins produced by the ISGs and non-ISGs, we found the former tends to have lower values of clustering coefficient and neighbourhood connectivity (Mann–Whitney U test: P =0.04 and 7.9E-03, Fig. 9D and F). This discovery indicates that the ISG products and some of their interacting proteins are less likely to be targeted by lots of proteins. It also supports the finding that the ISG products are involved in many shortest paths for nodes but are away from hubs or bottlenecks in the network. To some extent, this location also increases the length of the average shortest paths through ISG products in the network (Fig. 9A).

Open in a separate window

Figure 9:

Differences in network properties. The included features are (A) average shortest path, (B) betweenness, (C) closeness, (D) clustering coefficient, (E) degree, (F) neighbourhood connectivity, (G) stress, and (H) topological coefficient. Mann–Whitney U tests are applied for hypothesis testing on the whole comparing data without sampling, and the results are provided in Supplementary Data S2. Here, the ISGs and non-ISGs are taken from dataset S2 (No. = 620 and 874) while the background human genes use samples from dataset S1 (No. = 10,836) (Table 5).

When investigating the association between IFN-α–induced gene stimulation and network attributes of gene products, we only found the feature of neighbourhood connectivity was underrepresented as the level of differential expression in the presence of IFN increases (PCC = −0.392, P =2.2E-04). This suggests that proteins produced by genes that are highly upregulated in response to IFN-α are further away from hubs in the PPI networks.

Features highly associated with the level of IFN stimulations

In this study, we encoded a total of 397 discrete and 7,046 categorical features covering the aspects of evolutionary conservation, nucleotide composition, transcription, amino acid composition, and network profiles. In order to find out some key factors that may enhance or suppress the stimulation of human genes in the IFN-α system, we compared the representation of discrete features of human genes with different but positive log₂(fold change). Two features on the co-occurrence of SLNPs and SLAAPs were not taken into consideration here as they were more subjective than the other discrete features and were greatly influenced by the number of sequence patterns. Upon the calculation of PCC and the result of hypothesis tests, we found 168 features highly associated with the level of IFN-α–triggered stimulations (Student t-tests: P < 0.05) (Supplementary Data S3). Among them, 118 features showed a positive correlation (Fig. 10) while the remaining 50 features showed a negative correlation (Fig. 11) with the change of upregulation in IFN-α experiments. Among these 168 features, the number of ORFs, alternative splicing results, and counts of exons used for coding were encoded from characteristics of the gene. Average dN/dS and average dS within human paralogues were encoded based on the sequence alignment results from Ensembl [58]. In total, 140 and 22 features were encoded from the genetic sequence and proteomic sequence, respectively. The last one, neighbourhood connectivity, was obtained from the network profile of a human interactome constructed based on experimentally verified data in the HIPPIE database [55].

Open in a separate window

Figure 10:

A total of 118 features positively associated with higher upregulation after IFN-α treatments. Features here are screened based on Pearson’s correlation coefficient (PCC) and results of Student t-tests (P < 0.05). Features with a higher PCC indicate a stronger positive correlation. Detailed results about PCC and hypothesis tests are provided in Supplementary Data S3.

Open in a separate window

Figure 11:

Fifty features negatively associated with higher upregulation after IFN-α treatments. Features here are screened based on Pearson’s correlation coefficient (PCC) and results of Student t-tests (P < 0.05). Features with a lower PCC indicate stronger negative correlation. Detailed results about PCC and hypothesis tests are provided in Supplementary Data S3.

In the positive group, “large” amino acid composition that includes the composition of 5 amino acids with geometric volume ranging from 163 to 173 cubic angstroms was ranked first for having the highest PCC at 0.593 (Student t-test: P =2.8E-09). This feature was not highlighted previously as it did not have a strong signal for discriminating the ISGs from non-ISGs (Mann–Whitney U test: P > 0.05). Similar phenomena were found on 87 features (64 positive correlations and 23 negative correlations) such as AG content, ApG content, and previously mentioned neutral amino acid composition. The strongest negative correlation between feature representation and IFN-α–triggered stimulations was found on the feature of 4-mer “CGCG” (PCC = −0.593, P =3.2E-09). This feature also showed a differential distribution between the ISGs and non-ISGs, providing useful information to distinguish the ISGs from non-ISGs. Similar phenomena were found on 81 features (54 positive correlations and 27 negative correlations) such as previously mentioned GC content, CpG content, and the usage of codon “GCG” coding for amino acid “A.”

Collectively, the biased effect on the basic composition of nucleotide sequences influences the correlation between the representation of sequence-based features and IFN-α–triggered stimulations. Human genes that show overrepresentation in more features listed in Fig. 10 are expected to be more upregulated after IFN-α treatments at least in the human fibroblast cells. Meanwhile, the underrepresentation of features listed in Fig. 11 also contributes to the level of upregulation in the IFN-α experiments.

Difference in feature representation of interferon-repressed genes and genes with low levels of expression

We grouped human genes into 2 classes based on their response to IFN-α in the human fibroblast cells. Genes significantly upregulated in IFN-α experiments were included in the ISG class, while those that did not were put into the non-ISG class. However, there is also another group of human genes downregulated in the presence of IFN-α (i.e., the IRGs). They were labelled as the non-ISGs but contain unique patterns that constitute an important aspect of the IFN response [3]. Some of these IRGs were not upregulated in any known type I IFN systems, and thus they have been placed in a refined non-ISG class for analyses and predictions. Additionally, a number of genes have insufficient levels of expression in the experiments to determine a fold change (i.e., ELGs). Here, we used the previously defined features to compare the ISGs from dataset S2 with the IRGs and ELGs divided from the background dataset S1 (Table 5).

As shown in Fig. 12, the IRGs are differentially represented to a lower extent in the majority of nucleotide 4-mer composition features than the ISGs, indicating the deficiency of some nucleotide sequence patterns in the coding region of IRGs. Note that many nucleotide 4-mer composition features are more suppressed in the ISGs than non-ISGs, although the differences are small. The biased representation of these features in the IRGs suggests that the IRGs have characteristics similar to the ISGs rather than non-ISGs. Additionally, there are a very limited number of features relating to evolutionary conservation, nucleotide sequence composition, or codon usage showing obvious differences between the ISGs and IRGs, but many of them are differentially represented when comparing the ISGs with non-ISGs. Therefore, involving the IRGs in the class of the non-ISGs will increase the risk for machine learning models to produce more false positives. However, there are some informative features differentiating the IRGs from ISGs. For example, compared to the ISGs, the IRGs are more enriched in CpGs (Mann–Whitney U test: P =5.6E-03), which is also mentioned in [59]. The IRGs tend to have higher closeness centrality and neighbourhood connectivity than the ISGs (Mann–Whitney U test: P =0.04 and 6.4E-06, respectively), suggesting that the IRGs are more central in the human PPI network and connected to key proteins with many interaction partners. Differences in some amino acid composition features between the ISGs and IRGs can also be observed in Fig. 12. Therefore, accurate predictability is still expected when using features extracted from protein sequences.

Open in a separate window

Figure 12:

Differential expressions of discrete features between different genes and their coded proteins. Mann–Whitney U tests are applied for hypothesis testing on the whole comparing data without sampling, and the results are provided in Supplementary Data S2. Here, the ISGs and non-ISGs are taken from dataset S2 (No. = 620 and 874); the IRGs and ELGs are taken from dataset S4 (No. = 1,006) and dataset S8 (No. = 2,217); the background human genes are from dataset S1 (No. = 10,836) (Table 5).

Figure 12 illustrates 161 features showing significant differences (Mann–Whitney U tests: P < 0.05) in the representation of the ISGs and ELGs. An estimated 82% of these features were also differentially represented between the ISGs and non-ISGs. Seventy-nine percent of these significant features displayed similar overrepresentation or underrepresentation in 2 comparisons: ISGs versus ELGs and ISGs versus non-ISGs. These ratios indicate that the majority of the ELGs are less likely to be ISGs based on their feature profile as well as their low expression levels in cells induced with IFN-α. Network analyses showed that the ELG products tended to have lower values of all calculated network features than ISG products, with the exception of topological coefficient. This means that the ELG products are less connected to other human proteins in the human PPI network. Particularly, their abnormal representation on the feature of average shortest paths indicates that some ELGs (e.g., vascular cell adhesion molecule 1 [VCAM1] and ubiquitin D [UBD]) may still have high connectivity in the human PPI network.

Implementation with machine learning framework

In this study, we encoded 397 discrete and 7,046 categorical features for the analyses. As excess of features will greatly increase the dimension of feature spaces and complicate the classification task for the classifier, we limited the number of SLNPs to the top 100 based on the adjusted P-value, and we expected these to be sufficient to provide a picture of short linear sequence patterns in the coding region of the canonical transcript. Accordingly, features measuring the co-occurrence status of multiple SLNPs were recalculated based on the selected 100 SLNPs. As a result, we prepared 518 features (Supplementary Data S5) for our machine learning framework. To reduce the impact of noisy data on classifications, we used only the refined ISGs and non-ISGs from dataset S2 for training and modelling.

Table 3 first shows the comparisons of prediction performance among different machine learning methods. The threshold is determined by maximising the value of the Matthews correlation coefficient (MCC). As the random forest (RF) classifier was built based on randomly selected samples and features [60], we repeated its modelling procedure 10 times. These initial comparisons showed that the support vector machine (SVM) [61] is superior to k-nearest neighbours (KNN) and RF [60]. The poor prediction performance of the best base classifier (SVM, area under the curve [AUC] = 0.6509) indicates that there are a number of poorly performing features hidden in the set. As some genes respond to IFNs in a cell-specific manner [2], it is hard to produce predictions unless we detect key discriminating features, which are robust to the change of biological environment.

Here, we considered 2 iterative strategies for selecting predictive features. The first one is the forward feature selection (FFS) [63] wherein features are added sequentially based on their individual performance. This strategy did not work well (Table 3) as the prediction performances were all poor when the feature was used individually (Supplementary Data S5). The second strategy is developed based on the backward feature elimination scheme but uses fewer iterations to achieve the end result, namely, AUC-driven subtractive iteration algorithm (ASI) (Fig. 15).

Open in a separate window

Figure 15:

Expression of the ELGs in different tissues. Expression data for 10 ELGs are collected from the Genotype-Tissue Expression project (https://gtexportal.org/) [65]. The tissues in red are not included in the Interferome database [24]. White boxes in the heatmap indicate that there are no data available for genes in the corresponding tissues. The overall expression level of these 10 ELGs is reflected via human perspective photo retrieved from Expression Atlas (https://www.ebi.ac.uk/gxa) [66].

Preprocessing using the ASI algorithm showed that there were at least 28% of features disrupting the prediction model. The loss of some of the individual nucleotide 4-mer feature seemed not to influence the performance of the classifier at this stage, but the similarities between IRGs and ISGs (Fig. 12), particularly in the 4-mer features, were a cause for concern when the model was used to predict new data, especially unknown IRGs.

When using the ASI algorithm, the number of disrupting features did not stabilise until the algorithm reached the 11th iteration. The remaining 74 features constituted our optimum feature set for predicting the ISGs (Table 4). Among them, 14 and 9 features displayed positive and negative correlations with the level of upregulation in IFN-α experiments (P < 0.05). During the procedure, the AUC kept increasing steadily and reached 0.7479 at the end (Table 3). The MCC also showed an overall improvement, although it fluctuated slightly during the last few iterations. By ranking the scores calculated by the prediction model, we found 68.1% of the 496 genes (equal to the number of ISGs in the training dataset) were successfully predicted as the ISGs. Fig.13B illustrates the distribution of prediction scores generated by the ASI-optimised model for human genes with different expressions in IFN-α experiments. Human genes with higher upregulation in IFN-α experiments tend to obtain higher prediction scores from our optimised machine learning model (PCC = 0.243, P =4.2E-10).

Open in a separate window

Figure 13:

The optimisation of the machine learning model with the ASI algorithm. (A) Change of the prediction models based on the one generated with all 518 features (poorly performing feature vector = 144, best MCC = 0.250, SN_496 = 0.615, and AUC = 0.6509). (B) Distribution of prediction scores generated by the ASI-optimised model for human genes with different expression levels in the IFN-α system. The ISGs and non-ISGs shown in (B) are randomly selected through an undersampling strategy [64] on dataset S2. The list of gene names can be found in Supplementary Data S1.

However, there were also some ISGs incorrectly predicted by our model, even though they were highly upregulated, for example, basic leucine zipper ATF-like transcription factor 2 (BATF2, prediction score = 0.34). The model produced 33 ISGs with a prediction score higher than 0.8, but this number for the non-ISGs reduced to 6, including 1 IRG (tripartite motif containing 59 [TRIM59]). The highest prediction score within the non-ISGs was found on ubiquitin conjugating enzyme E2 R2 (UBE2R2, prediction score = 0.88). It contains many features similar to the ISGs but was not differentially expressed in the presence of IFN-α in the human fibroblast cells [3]. The lowest prediction score within ISGs was found on cap methyltransferase 1 (CMTR1, prediction score = 0.12) due to the weak signal from its features. For instance, CMTR1 protein does not contain any ISG-favoured SLAAPs listed in Table 1. The influence of the IRGs on the prediction was reflected in the training dataset but was not significant. Compared with human genes not differentially expressed in the IFN-α experiments (non-ISGs but not IRGs), there were slightly more IRGs unsuccessfully classified when using a threshold of 0.549 (Pearson's chi-squared tests: M₁ = 27%, M₂ = 24%, P > 0.05).

Generation of discrete features

We encoded 397 discrete features from aspects of evolution, nucleotide composition, transcription, amino acid composition, and network properties. Original values of these features for our compiled 10,836 human genes are accessible at [90].

From the perspective of evolution, we used the number of transcripts, ORFs, and count of exons used for coding to quantify the alternative splicing process. Genes with more transcripts and ORFs have higher alternative splicing diversity to produce proteins with similar or different biological functions [33, 91, 92]. Frequent use of protein-coding exons indicates more complex alternative splicing products [93]. Here, duplication and mutation features were measured by the number of within-species paralogues and substitutions [34, 35]. These data were collected from BioMart (RRID:SCR_002987) [58] to assess the selection on protein sequences and mutational processes affecting the human genome [94].

From the perspective of nucleotide composition, we calculated the percentage of adenine, thymine, cytosine, guanine, and their 4-category combinations in the coding region of the canonical transcript. The first category measured the proportion of 2 different nitrogenous bases out of the implied 4 bases (e.g., GC content). The second category also focused on the combination of 2 nucleotides but added the impact of phosphodiester bonds along the 5′ to 3′ direction (e.g., CpG content) [95]. The third category calculated the occurrence frequency of 4-mers (e.g., “CGCG” composition to involve some positional resolution) [41]. The last category considered the co-occurrence of SLNPs. From the perspective of transcription, we calculated the usage of 61 coding codons and 3 stop codons in the coding region of the canonical transcripts. Codon usage biases are observed when there are multiple codons available for coding 1 specific amino acid. They can affect the dynamics of translation and thus regulate the efficiency of translation and even the folding of the proteins [40, 96].

From the perspective of amino acid composition, we calculated the percentage of 20 standard amino acids and their combinations based on their physicochemical properties [46]. Patterns in the amino acid level are considered to have a direct impact on the establishment of biological functions or to reflect the result of strong purifying selection [47]. Based on the chemical properties of the side chain, we grouped amino acids into 7 classes, including aliphatic, aromatic, sulphur, hydroxyl, acidic, amide, and basic amino acids. We also grouped amino acids based on geometric volume, hydropathy, charge status, and polarity but found some overlaps among these features. For instance, amino acids with basic side chains are all positively charged. Aromatic amino acids all have large geometric volumes (volume >180 cubic angstroms). Likewise, we also considered the co-occurrence of short linear sequence patterns at the protein level. These co-occurring SLAAPs may relate to potential mechanisms regulating the expression of the ISGs [97].

To infer network properties for the gene products, we constructed a human PPI network based on 332,698 experimentally verified interactions (confidence score >0.63) from HIPPIE (RRID:SCR_014651) [56]. Nodes and edges of this network are provided at our web server. Eight network-based features, including the average shortest path, closeness, betweenness, stress, degree, neighbourhood connectivity, clustering coefficient, and topological coefficient, were calculated from this network. Isolated nodes or proteins were not included in our network and were assigned zero values for all these 8 features. The shortest path measures the average length of the shortest path between a focused node and others in the network. Closeness of a node is defined as the reciprocal of the length of the average shortest path. Proteins with a low value of the shortest paths or closeness are close to the centre of the network. Betweenness reflects the degree of control that 1 node exerted over the interactions of other nodes in the network [98]. Stress of a node measures the number of shortest paths passing through it. Proteins with a high value of betweenness or stress are close to the bottleneck of the network. Degree of a node counts the number of edges linked to it while neighbourhood connectivity reflects the average degree of its neighbours. Proteins with high values of degree or neighbourhood connectivity are close to the hub of the network. They are considered to play an important role in the establishment of the stable structure of the human interactome [99]. Clustering and topological coefficient measure the possibility of a node to form clusters or topological structures with shared neighbours. The former coefficient can be used to identify the modular organisation of metabolic networks [100] while the latter one may be helpful to find out virus mimicry targets [53].

Generation of categorical features

In this study, categorical features were used to check the occurrence of short linear sequence patterns in the genome and proteome. SLNPs constructed in this study contained 3 to 5 random nucleotides, producing 708,540 alternative choices. SLNPs with no restrictions on their first or last position were not taken into consideration as their patterns could be expressed in a more concise way. A SLNP was picked out to encode a binary feature when its occurrence level in the coding region of the canonical ISG transcripts was significantly higher than that for the non-ISGs (Pearson's chi-squared test: P < 0.05). SLAAPs were constructed with 3 to 4 fixed amino acids separated by putative gaps. The gap could be occupied by at most 1 random amino acid, producing 1,312,000 alternative choices. Likewise, binary features were prepared for SLAAPs showing significant enrichment in the ISG products than in the non-ISG products (Pearson's chi-squared test: P < 0.05). Since there were lots of results rejecting the null hypothesis, we adopted the Benjamini–Hochberg correction procedure to avoid type I error [43]. Additionally, we also encoded 2 features to check the co-occurrence or absence of multiple SLNPs and SLAAPs. This co-occurrence status might be a better representation of functional sites composed of short stretches of adjacent nucleobases or amino acids surrounding SLNPs or SLAAPs [47].

Assessment of associations between feature representation and IFN-triggered stimulations

We obtained 8,619 human genes with expression data from the OCISG database [3]. In total, 4,111 of them were annotated with a positive log₂(fold change) ranging from 0 to 12.6, which meant they were upregulated after IFN-α treatments in the human fibroblast cells. In order to measure the average level of feature representation (AREP) for genes with similar expression during IFN stimulations, we introduced a 0.1-length sliding window to divide the data into 126 bins with different log₂(fold change). Here, PCC was introduced to test the association between the representation of discrete features and IFN-α–triggered stimulation (log₂(fold change) >0). It can be formulated as

(1)

where n is the number of divided parts that equals 126 in this study; and are the value of log₂(fold change) and AREP in the ith part; and are the mean and standard deviation of log₂(fold change), which are set as 6.4 and 3.7, respectively, in this study; and and are the mean and standard deviation of 126 AREPs that reflect the representation of the considered feature. To make fair comparisons among features with different scales, we normalised them based on the major value of their representations:

(2)

where and are the lower and upper bound representing the 5th and 95th percentile within representation values for the target feature. The representation of feature was considered to have a stronger positive/negative association with IFN-α–triggered stimulations if the PCC calculated from the normalised features was closer to 1.0/−1.0 and the P-value calculated by the Student t-test was lower than 0.05.

Machine learning and optimisation

We designed a machine learning framework for the prediction of ISGs. First, all features were encoded and normalised based on their major representations (Eq. 2). Then we used an undersampling procedure [64] to generate a balanced dataset from dataset S2 for training and modelling. The SVM with radial basis function [61] was used as the basic classifier. It maps the normalised feature space to a higher dimension to generate a space plane to better classify the majority of positive and negative samples. In order to avoid overfitting [77] and make it easier for the SVM model to generate an appropriate classification plane that involved fewer false positives and false negatives, here we propose a subtractive iteration algorithm driven by the change of AUC. This algorithm is developed based on the traditional backward feature elimination method [63] but uses fewer iterations to filter out poorly performing features (Fig. 16). In each iteration, we traversed the features and removed those that did not improve the AUC of the prediction results. In the testing procedure, we encoded the optimum features for testing samples and placed them in the optimised feature space. Samples with longer distance to the optimised classification plane indicated a stronger signal of being the ISGs or non-ISGs. They were more likely to get higher prediction scores (close to 0 or 1) from the SVM model.

Open in a separate window

Figure 16:

The pseudo-code of the AUC-driven subtractive iteration algorithm.

Performance evaluation

In this study, the prediction results were evaluated with 3 threshold-dependent criteria, including sensitivity, specificity, and MCC [101], and 2 threshold-independent criteria: SN_n and AUC. Sensitivity and specificity were used to assess the quality of the machine learning model in recognising ISGs and non-ISGs, respectively, while MCC provided a comprehensive evaluation for both positives and negatives. The number of “n” in the SN_n criterion was determined based on the number of ISGs used for testing. It was used to measure the upper limit of the prediction model as well as to check the existence of important false positives close to the class of ISGs from the perspective of data expression. Finally, AUC was a widely used criterion to evaluate the prediction ability of a binary classifier system. The group of interest was almost unpredictable in a specific binary classifier system if the AUC of the classifier was close to 0.5.

The authors thank Andrew Davison, Suzannah Rihn, and Sam Wilson for helpful discussions and recommendations and Scott Arkison for help setting up the website.