key: cord-0917822-0osna8qb
authors: Roy, Sanjiban Sekhar; Taguchi, Y.-H.
title: Identification of genes associated with altered gene expression and m6A profiles during hypoxia using tensor decomposition based unsupervised feature extraction
date: 2021-04-26
journal: Sci Rep
DOI: 10.1038/s41598-021-87779-7
sha: ee5c8f816b2c14acf29caf1f052974a11484a180
doc_id: 917822
cord_uid: 0osna8qb

Although hypoxia is a critical factor that can drive the progression of various diseases, the mechanism underlying hypoxia itself remains unclear. Recently, m6A has been proposed as an important factor driving hypoxia. Despite successful analyses, potential genes were not selected with statistical significance but were selected based solely on fold changes. Because the number of genes is large while the number of samples is small, it was impossible to select genes using conventional feature selection methods with statistical significance. In this study, we applied the recently proposed principal component analysis (PCA), tensor decomposition (TD), and kernel tensor decomposition (KTD)-based unsupervised feature extraction (FE) to a hypoxia data set. We found that PCA, TD, and KTD-based unsupervised FE could successfully identify a limited number of genes associated with altered gene expression and m6A profiles, as well as the enrichment of hypoxia-related biological terms, with improved statistical significance.

Hypoxia 1 , also known as tissue hypoxia, is a serious symptom with various causes. For example, hypoxia could result in death, such as in the case of COVID-19, a serious pandemic 2 . Hypoxia also plays a critical role in cancer 3 . Both brain hypoxia 4 and lung cell hypoxia 5 can be fatal. Despite the significance of hypoxia, the critical factors of hypoxia are not yet fully understood 6 . Recently, m6A was reported to be a newly discovered regulator of hypoxia 7 . Wang et al. 8 found that many genes are simultaneously associated with altered m6A and gene expression profiles in hypoxia. Although the investigations were successful, there was one methodological issue with their study; they selected genes associated with altered m6A and gene expression in hypoxia without determining statistical significance. They selected genes based on fold change (FC). Usually, only using FC to select altered expression or any other measurements might be erroneous because a sufficiently large FC might be observed simply by chance when a large number of candidates are considered. In their analysis, all human genes (as many as a few tens of thousands) and whole genome m6A were considered. In this case, if the FC was not validated statistically, a sufficiently large FC might have been observed simply by chance. The genes associated with altered m6A and gene expression based on statistical significance could not be identified because of the small number of samples; there were only four time points (including the control) measured without any replicates. If we consider the large number of genes as well as m6A peaks in the genome, it is unlikely that four samples are enough to achieve statistical significance; small samples result in larger P-values, whereas a large number of genes and m6A peaks result in relatively larger P-values. In this study, we applied principal component analysis (PCA) and tensor decomposition (TD)-based unsupervised feature extraction (FE) to select genes associated with altered m6A as well as gene expression in hypoxia to determine statistical significance. Enrichment analyses of selected genes are reasonable and consistent with previous findings 8 and can now be supported with statistical significance. Thus, not only were the critical roles of m6A in hypoxia validated but also the usefulness of PCA-and TD-based unsupervised FE in the case where there are very few samples with a large number of variables.

There are a limited number of genomic studies using TD 9, 10 . Fang proposed tightly integrated genomic and epigenomic data mining using TD 11 applied TD to DNA copy-number alterations 18 (a few hundred samples). All methods other than that used by Li et al. included as few as 53 genes and required as many as several hundred samples, whereas our methods generally require only a few samples (in this study as few as eight samples). To our knowledge, ours is the only method applicable to a data set that includes a few samples with as many as 10 4 genes. Figure 1 shows the flow chart of analyses performed in this study.

PCA based unsupervised FE applied to gene expression profiles. Because gene expression profiles are measured solely over four time points, it is formatted as a matrix x it ∈ R N×4 , and PCA-based unsupervised FE was applied to x it . Figure 2 shows the PC loading attributed to time points. Because the second PC loading is mostly correlated with time, we decided to employ the second PC score, u 2i , in order to attribute P-values to genes i. Using Eq. (7) with ℓ = 2 , P i s are attributed to gene i. Then 52 is (genes) with associated corrected P-values less than 0.01, were selected. Table 1 shows the enrichment terms in "KEGG 2019 Human" categories in Enrichr for 52 selected gene symbols. Not all terms are related to hypoxia, whereas a term such as "Oxidative phosphorylation" is known to be related to hypoxia 19 . "Cardiac muscle contraction" is also known to be related to hypoxia 20 . Retrograde endocannabinoid signaling is known to be related to hypoxia 21 . Although three representative neurodegenerative diseases, "Parkinson's disease, " "Alzheimer's disease, " and "Huntington's disease", are listed, hypoxia is known to be related to neurodegenerative diseases 22 . Glycolysis is also related to hypoxia 23 . Most importantly, HIF-1, a hypoxia-inducible factor, is listed. There are additional identified enrichments that can support the success of PCA-based unsupervised FE. Although they are not always top-ranked, the 52 identified genes are also known to be up/downregulated in independent hypoxia experiments ( Table 2 ). In the "GO Biological Process 2018" category in Enrichr, various glucose/glucogenesis related terms are enriched (Table 3) . These results suggested that the analyses performed by PCA-based unsupervised FE were successful.

TD-based unsupervised FE applied to m6A profiles. Although we successfully applied PCA-based unsupervised FE to some gene expression profiles of hypoxia, the identification of the relationship between altered gene expression and hypoxia was not the primary purpose of this study. Instead, its purpose was to identify the relationship between m6A profiles and hypoxia. To identify the relationship between hypoxia and m6A profiles, HOSVD was applied to x ktj , as described in "Materials and methods". The left panel in Fig. 3 shows the singular value vectors attributed to time points; the second singular value vector is most significantly correlated with time. Remarkably, u 2t is almost identical to v 2t , which is the 2nd PC loading attributed to the time points when PCA was applied to x it (right panel in Fig. 3 ). Considering that gene expression and m6A profiles are distinct from each other, this coincidence between u 2t , which is attributed to m6A profiles, and v 2t , which is attributed to gene expression profiles, suggests that our analysis correctly detects the regulatory relationship between m6A and the gene expression profile. In addition to the fact that u 2t is most significantly correlated with time points, u 2j s have opposite signs between j = 1 and j = 2 (not shown here), which means that ℓ 3 = 2 is associated with the distinction between the input and m6A. We then determined which G(ℓ 1 , 2, 2) has the largest absolute value to determine which u ℓ 1 k is used to select genomic regions, k (Table 4 ). Because it is obvious that G(2, 2, 2) has the largest absolute value, we decided to employ u 2k to select genomic regions. P k s are attributed to k using Eq. (9) with ℓ 1 = 2 . Then 106 ks (genomic regions 25,000 nucleotides in length, see "Materials and methods") associated with corrected P-values less than 0.01, were selected. These 106 genomic regions included 196 unique gene symbols that were uploaded to Enrichr to evaluate enrichment. In contrast to the 52 genes identified by PCA-based unsupervised FE applied to gene expression, no KEGG pathway terms or GO BP terms were enriched in these 196 gene symbols. Nevertheless, there are some hypoxia experiments in which genes with altered expression are enriched in 196 gene symbols (Table 5) . Therefore, even though 196 genes are less biologically significant than the 52 genes identified in the gene expression analysis, they still have some potential to be related to hypoxia.

Integrated analysis of gene expression and m6A profiles using KTD-based unsupervised FE. Since TD-based unsupervised FE applied to m6A profiles was not fully successful, we needed to employ more advanced methodology: Kernel TD-based (KTD) unsupervised FE. HOSVD was applied to x tjt ′ j ′ , as described in "Materials and methods". Figure 4 shows the results that are consistent with the results obtained by non-integrated analysis (Figs. 2, 3) . The second singular value vectors, u 2t and u 2t ′ , are most consistent with time points; it is coincident with the second PC loading attributed to gene expression, and the second singular value vectors attributed to m6A are coincident with time points. u 2t and u 2t ′ are also identical; it is coincident with the second PC loading attributed to gene expression, and the second singular value vectors attributed to m6A are identical. In addition, u 2j has the opposite sign between j = 1 (control) and j = 2 (m6A). Then u ℓ 1 i was computed using Eq. (16), with ℓ 1 = 2 , and u ℓ 2 ℓ 3 k was computed using Eq. (17), with ℓ 2 = ℓ 3 = 2 . P i and P k are attributed to i and k, respectively, with Eqs. (18) and (19) , respectively. The 53 is (genes) and 128 ks (genome regions) associated with adjusted P-values less than 0.01 were selected. Two hundred gene symbols were retrieved from 128 genomic regions, as previously described.

The 53 and 200 gene symbols were uploaded to Enrichr. www.nature.com/scientificreports/ enrichment terms were identified, whereas no terms were identified when TD-based unsupervised FE was applied to the m6A profile. Table 7 shows the results of the "Disease Perturbations from GEO down/up" category in Enrichr. Compared with Tables 2 and 5, although enrichment in the "Disease Perturbations from GEO down" for m6A is missing, it still has enrichment for both gene expression and m6A profiles. Table 8 shows the enrichment in the "GO Biological Process 2018" category of Enrichr. Compared with Table 3 , enrichment for gene expression does not change, and enrichment for m6A is identified, whereas it was not identified when TD-based unsupervised FE was applied to the m6A profile. Thus, KTD-based unsupervised FE improved the enrichment of m6A profiles without affecting the enrichment for gene expression. www.nature.com/scientificreports/ Additional improvement from PCA-and TD-based unsupervised FE to integrated analysis using KTD-based unsupervised FE identified significant associations between gene expression and m6A (Table 9 ). When PCA-and TD-based unsupervised FE were separately applied to gene expression and m6A profiles, 52 and 196 genes were identified, respectively. The number of common genes between them was seven. However, integrated analysis of gene expression and m6A profiles with TKD-based unsupervised FE identified 53 genes for gene expression Table 3 . "GO Biological Process 2018" category of Enrichr for 52 genes selected by applying PCA-based unsupervised FE to gene expression profiles. Six glucose-related terms with adjusted P-values less than 0.05 are listed.

Overlap P-value Adjusted P-value www.nature.com/scientificreports/ and 200 genes for m6A profiles. The number of common genes increased to 12. Although we cannot estimate their significance very accurately, if we can tentatively assume that there are 20,000 human genes in total, both coincidences are significant, and coincidence in KTD-based unsupervised FE is more significant.

In conclusion, integrated analysis of gene expression and m6A profiles using KTD-based unsupervised FE substantially increased the results compared with applying PCA-and TD-based unsupervised FE separately to gene expression and m6A profiles, respectively. KTD-based unsupervised FE could identify the relationship of gene expression and m6A with hypoxia simultaneously for the first time in a statistically significant manner.

Although we have shown that integrated analysis of gene expression and m6A profiles simultaneously identified the relationship between hypoxia and gene expression as well as hypoxia and m6A profiles, if other simpler conventional methods can achieve similar performances, it is useless to employ complicated methods such as KTD-based unsupervised FE. To confirm that other conventional methods cannot identify similar relationships, we applied a few conventional feature selection methods. As can be seen in the following text, no conventional feature selections were found to be useful.

Linear regression analysis. Linear regressions were applied to the gene expression profiles, x it , and m6A profiles, x ktj . No genes or genomic regions was associated with adjusted P-values less than 0.05, respectively; thus, no genes or genomic regions were associated with adjusted P-values less than 0.01 either.

Although we tried to apply SAM to gene expression profiles, x it , and m6A profiles, x ktj , we found that SAM requires at least two replicates for each class. In this study, there are no replicated classes in four classes in x it or eight classes in x ktj ; therefore, we could not apply SAM to these data sets. www.nature.com/scientificreports/ Limma. Limma was applied to the gene expression profiles, x it , and m6A profiles, x ktj . No genes but 72252 genomic regions were associated with adjusted P-values of less than 0.01, respectively. Therefore, limma was not useful.

Random forest. Random forest was applied to gene expression profiles, x it , and m6A profiles, x ktj . Four hundred and eighty genes and 722 genomic regions had non-zero importance, respectively. Thus, random forest successfully selected a reasonable number of genes and genomic regions. Nevertheless, no hypoxia-related biological terms were enriched in the 480 genes or gene symbols included in the 722 genomic regions. Thus, random forest was not a useful method. www.nature.com/scientificreports/

One might wonder why linear regression, SAM, limma, and random forest failed to select genes associated with altered gene expression, genomic regions associated with m6A profiles in hypoxia, or genes biologically related to hypoxia. This is because it is a very difficult problem. There are more than 17140 genes as well as 123817 genomic regions, whereas the number of samples measured was four and eight, respectively, which were too small to obtain sufficiently significant P-values. These numbers of genes and genomic regions were too large to obtain significant P-values; although random forest is free from P-values, too small sample numbers often prevent random forest from obtaining results that are not obtainable by chance.

To demonstrate how the KTD-based unsupervised FE outperforms the other four methods, we applied them to two synthetic data sets with N = 1000 and N 1 = 10 . When linear regression was applied to the 1st synthetic data set, there were no is associated with adjusted P-values less than 0.01. When limma and random forest were applied to the 1st synthetic data set, there were as many as 500 is associated with adjusted P-values less than 0.01. Thus, neither linear regression nor limma was useful.

The problem is when time points are regarded as a quantitative property (i.e., in linear regression); in this case, eight samples were too small to give significant P-values because there are 1000 features on which P-values must be corrected by considering multiple comparison criteria. However, if they are classified into eight classes with one replicate, too many is were regarded as distinct between eight classes because those values that are constant between any pairs of eight classes are unlikely to be fulfilled when using a null hypothesis.

However, when TD-based unsupervised FE was applied to the first synthetic data set, the situation was very different. Figures 5 and 6 show the two combinations of the u ℓ 2 ,j , u ℓ 3 k , and G(ℓ 1 ℓ 2 ℓ 3 ) . Figure 5 Fig. 7 ). However, Fig. 6 ( ℓ 2 = ℓ 3 = 2 ) corresponds to x ijk , N 1 < i ≤ 2N 1 because x ij1 = −x ij2 , which coincides with u ℓ 3 1 = −u ℓ 3 2 (see Fig. 7 ). Thus, in contrast to other supervised methods, TD-based unsupervised FE can detect the distinction between x ij1 = x ij2 and x ij1 = −x ij2 . To determine whether TD-based unsupervised FE can correctly identify x ijk s, we need to attribute P-values to is. Because |G(ℓ 1 , 2, 1)| and |G(ℓ 1 , 2, 2)| have the largest values when ℓ 1 = 3 and ℓ 1 = 2 , respectively, we decided to attribute P-values to is using Eq. (13) by assigning ℓ 1 = 3 and ℓ 1 = 2 , respectively. Computed P-values were corrected, and is associated with adjusted P-values less than 0.01 were selected. Table 10 shows the performance of TD-based unsupervised FE applied to the first synthetic data set. It perfectly selects i coincident with Eq. (1). As a result, the reason why only PCA and TD-based unsupervised FE could select a reasonable number of genes, while other methods failed, is because PCA and TD-based Table 8 . "GO Biological Process 2018" category of Enrichr for 53 genes (based on gene expression) and 200 genes (based on m6A) selected by applying KTD-based unsupervised FE to integration of gene expression and m6A profiles. Nine glucose-related terms with adjusted P-values less than 0.05 are listed.

Overlap P-value Adjusted P-value

Glycolytic process through glucose- Table 9 . Confusion matrices of selected genes between gene expression and m6A profiles. PCA-and TD-based unsupervised FE were separately applied to gene expression and m6A profiles, or KTD-based unsupervised FE was applied to the integration of gene expression and m6A profiles Next, we attempted to demonstrate how KTD-based unsupervised FE integrates two data sets in order to identify common features between the two. Figures 8 and 9 show singular value vectors obtained when KTDbased unsupervised FE was applied to x jk 1 k 2 j ′ k ′ 1 k ′ 2 . Figure 8 shows u 1j , u 1k 1 , u 1j ′ , and u 1k ′ 1 , which is coincident with x ijk , i ≤ N 1 , and 2N 1 < i ≤ 3N 1 and x ′ ijk , i ≤ N 1 , and 4N 1 < i ≤ 5N 1 , since both singular value vectors and x ijk and x ′ ijk for these is increase as j increases and do not depend on k (see Fig. 10 ). To determine whether we can select these x ijk and x ′ ijk in these is, we reproduced singular value vectors attributed to is using Eq. (22) with ℓ 1 = ℓ 2 = 1 or Eq. (23) with ℓ 4 = ℓ 5 = 1 . Then P-values were attributed to is using Eq. (24) or Eq. (25) . Computed P-values were corrected, and is associated with adjusted P-values less than 0.01 were selected. Table 11 shows a perfect performance. Figure 9 shows u 1j , u 2k 1 , u 1j ′ , and u 2k ′ 1 , which is coincident with x ijk , N 1 < i ≤ 2N 1 , and 3N 1 < i ≤ 4N 1 , and x ′ ijk , N 1 < i ≤ 2N 1 , and 5N 1 < i ≤ 6N 1 , since u 1j and u 1j ′ increase as j increases, and u 2k 1 and u 2k ′ 1 have opposite www.nature.com/scientificreports/ signs between k 1 = k ′ 1 = 1 and k 2 = k ′ 2 = 2 , while x ij1 and x ′ ijk1 for these is increase as j increases, and x ij2 and x ′ ij2 for these is decrease as j increases (see Fig. 10 ). To determine whether we can select these x ijk and x ′ ijk in these is, we reproduced singular value vectors attributed to is using Eq. (22) with ℓ 1 = 1, ℓ 2 = 2 or Eq. (23) with ℓ 4 = 1, ℓ 5 = 2 . Then, P-values were attributed to is using Eq. (24) or Eq. (25) . Computed P-values were corrected, and is associated with adjusted P-values less than 0.01 were selected. Table 11 shows a perfect performance. Table 12 shows the confusion matrix between is selected for x ijk and those selected for x ′ ijk . This corresponds to Table 9 , where genes were selected based on gene expression and m6A profiles. This might be the reason why KTD-based unsupervised FE could identify a significantly overlapping set of genes between gene expression and m6A profiles.

Although TD-and KTD-based unsupervised FE can outperform conventional supervised methods when applied to a small number of samples with a large number of features, TD-and KTD-based unsupervised FE have yet another advantage: j dependence is not monotonic (see the open red triangles in Figs. 2, 3, and 4) . Such a non-linear dependence on j cannot be assumed by supervised methods in advance. Wrongly assumed j dependence results in decreased feature selection performance. This is another reason why PCA-, TD-, and KTD-based unsupervised FE can outperform other conventional supervised feature selection methods.

In order to see if our findings are robust, we tried to find alternative data sets in which gene expression and m6A were simultaneously measured for hypoxia, but we could not find any such data sets. Thus, we employed GSE120860, in which only m6A was measured. TD-based unsupervised FE applied to these data sets gave us 54 genes associated with adjusted P-values less than 0.01 ( ℓ 2 = ℓ 3 = 1 and ℓ 4 = 2 were selected, and u 4i was used Figure 7 . Schematic figure that explains the correspondence between x ijk , i ≤ N 1 (red solid arrows) and u 2j , u 1k and that between x ijk , N 1 < i ≤ 2N 1 (red broken arrows) and u 2j and u 2k , respectively. www.nature.com/scientificreports/ to attribute P-values to gene i with Eq. (11) because G(4, 1, 1, 2) has the largest absolute value given ℓ 2 = ℓ 3 = 1 and ℓ 4 = 2 ). Uploading these genes to Enrichr did not identify any terms associated with both hypoxia and significant P-values. As shown in Table 5 , when genes were selected by m6A only, there were fewer significant terms. Therefore, we may need to have alternative data sets associated with both gene expression and m6A simultaneously in order to validate the robustness of our results. One might wonder why we did not compare the proposed methods with conventional unsupervised methods using PCA and TD but only with supervised methods. The reasons for this are as follows. Although there are many papers whose titles include "Feature selection using principal component analysis", feature selections in these papers mean selecting limited numbers of latent vectors generated by PCA or TD. Thus, they are not applicable to the present study, which needed to select not generated features but original ones (i.e., genomic regions). Although there are a few studies that aim to select original features, and not generated latent vectors, they did not attribute P-values to the features, which would have allowed us to evaluate the significance of the feature selections. For example, Song et al. 24 selected a limited number of original features associated with relatively larger absolute values of eigenvectors, and no P-values were attributed to the individual original features. The purpose of the present study is not simply to select features but to evaluate the significance of selected features; as denoted above, Song et al's study could not help us to evaluate the significance of the feature selections. This is why we did not compare our method with other unsupervised methods using PCA or TD but compared ours with the supervised methods that could give us P-values, by which we could evaluate the significance of the feature selections.

In this study, we applied PCA-, TD-, and KTD-based unsupervised FE to gene expression and m6A profiles in hypoxia. Although these methods identified a limited number of genes significantly related to hypoxia, other conventional methods failed. To understand why PCA-, TD-, and KTD-based unsupervised FE could outperform other conventional methods, we applied these methods to synthetic data sets with small numbers of samples and large numbers of features. As a result, we successfully reproduced the superior performance of TD-and KTD-based unsupervised FE over other conventional methods. Thus, the superiority of PCA-, TD-, and KTDbased unsupervised FE is possibly due to having a small number of samples with a large number of features. In conclusion, despite the limitations of previous studies, we validated a set of genes associated with altered gene expression and m6A profiles in hypoxia in a statistically significant manner. 

ijk , i ≤ N 1 , and 4N 1 < i ≤ 5N 1 (red solid arrows) and u 1j and u 1k and that between www.nature.com/scientificreports/ length were obtained. For gene expression, four profiles included in GSE141941_normoxiaVShypoxia6h.12h.24h_ RNA-seq.PROCESSED.DATA.xlsx, which is also available as a part of the Supplementary Information, were employed. Eight files for m6A were composed of four time points (including the control), time input, or treated files. Four profiles for gene expression were composed of four times points as well.

As an alternative data set, we employed GSE120860, in which only the m6A profile was measured. We downloaded 16 bed files provided in the Supplementary Information section in GEO, which correspond to four healthy controls and four patients, of which the tumor and paratumor were measured.

Synthetic data set. In order to demonstrate how well KTD-based unsupervised FE can work when there are a small number of samples associated with a large number of features (observations) and where other conventional supervised methods fail, we prepared a synthetic data set. It has N variables attributed to eight samples whose number is the same as that of the m6A profiles. In the first data set, we aimed to demonstrate the performance when KTD-based unsupervised FE is applied to a single data set. It is composed of x ijk ∈ R N×4×2 , where j corresponds to the first to fourth time points, and k = 1 and k = 2 correspond to two distinct experimental conditions. ǫ ijk obeys N (0, 1 2 ) , where N (µ, σ ) is a Gaussian distribution with a mean of µ and a standard deviation of σ . For i ≤ N 1 , the two conditions have the same dependence on time points, whereas for N 1 < j ≤ 2N 2 , the two conditions have opposite time point dependence.

In the second dataset, we aimed to demonstrate how KTD-based unsupervised FE can identify features that share the same time point dependence on two measurements represented as two tensors, x ijk and x ′ ijk ∈ R N×4×2 , which obey Eq. (1) for i ≤ 2N 1 , and for i > 2N 1 ,

ijk , and x ′ (i+4N 1 )jk share the same time-point dependence.

PCA-based unsupervised FE applied to gene expression. Although the details of PCA-based unsupervised FE have been described in a recently published book 25 , we briefly outline this method. Suppose we have gene expression profiles as a matrix, x it ∈ R N×4 , which represents the gene expression of the ith gene at the tth time point. The ℓ th PC score attributed to gene i, u ℓi ∈ R N , can be obtained as the ith component of the ℓ th eigenvector, u ℓ , of the gram matrix XX T ∈ R N×N , where X is an N × 4 matrix composed of x it , and X T is a transposed matrix of X as

(4) XX T u ℓ = ℓ u ℓ . www.nature.com/scientificreports/ In contrast, the PC loading attributed to the tth time point can be computed as the tth component of the ℓ th PC loading vector, v ℓ ∈ R 4 , which can be computed as This is also the ℓ th eigenvector of X T X ∈ R 4×4 because In order to select genes, we first need to determine which v ℓ is associated with time dependence. After identifying the time-dependent v ℓ , we attribute P-values to the ith gene using u ℓi by assuming that u ℓi follows a Gaussian distribution (null hypothesis)

where P χ [> x] is the cumulative χ 2 distribution in which the argument is larger than x, and σ ℓ is the standard deviation. The computed P-values were corrected by the BH criterion 25 , and genes associated with adjusted P-values less than 0.01 were selected.

TD-based unsupervised FE applied to the m6A profile. Although the details of TD-based unsupervised FE are described in a recently published book 25 , we briefly outline this method. For GSE141941, suppose that a tensor x ktj ∈ R K×4×2 represents the m6A of the kth genomic region at the tth time point of input (control, j = 1 ) or m6A ( j = 2 ) sample. Individual genomic regions are 25,000-nucleotide sequence length regions sequentially defined over the whole genome without overlaps and adjusted with each other. Higher-order singular value decomposition 25 (HOSVD) was applied to x ktj , and TD was obtained as where G ∈ R K×4×2 is the core tensor, u ℓ 1 k ∈ R K×K , u ℓ 2 t ∈ R 4×4 , and u ℓ 3 j ∈ R 2×2 are singular value matrices and orthogonal matrices.

To select genomic regions that are associated with time dependence and to distinguish between input and m6A treatment, we need to specify which u ℓ 2 t and u ℓ 3 j are associated with time dependence and the distinction between control (input) and m6a, respectively. Once u ℓ 2 t and u ℓ 3 j are fixed, we attempt to find G(ℓ 1 ℓ 2 ℓ 3 ) with the largest absolute value, given ℓ 2 and ℓ 3 . Finally, we attribute the P-value to the kth genomic region by assuming that u ℓ 1 k obeys a Gaussian distribution (null hypothesis) as where σ ℓ 1 is the standard deviation. The computed P-values were corrected by the BH criterion, and genes associated with adjusted P-values less than 0.01 were selected.

For GSE126860, m6A profiles were formatted as x ijkm ∈ R N×4×2×2 , which represents the m6A profiles of the ith gene at the jth subject of the kth group ( k = 1 : annotated as 0204 in GEO, k = 2 : annotated as patient in GEO) of the mth tissue ( m = 1:tumor, m = 2:paratumor). HOSVD was applied, and we obtained where G ∈ R N×4×2×2 is the core tensor, and u ℓ 1 i ∈ R N×N , u ℓ 2 j ∈ R 4×4 , and u ℓ 3 k , u ℓ 4 m , ∈ R 2×2 are singular value matrices and orthogonal matrices.

In order to select genes that are associated with the distinction between tumor and paratumor, but independent of subjects as well as groups, we need to specify which u ℓ 4 m are associated with the distinction between the tumor and paratumor, but u ℓ 2 j and u ℓ 3 k take constant values. Once u ℓ 2 j , u ℓ 3 k , and u ℓ 4 m are fixed, we then attempt to find that G(ℓ 1 ℓ 2 ℓ 3 ℓ 4 ) with the largest absolute value, given ℓ 2 , ℓ 3 , and ℓ 4 . Finally, we attribute the P-value to the ith gene by assuming that u ℓ 1 i obeys a Gaussian distribution (null hypothesis) as where σ ℓ 1 is the standard deviation. The computed P-values were corrected by the BH criterion, and genes associated with adjusted P-values less than 0.01 were selected.

When HOSVD was applied to the 1st synthetic data set, we obtained www.nature.com/scientificreports/ where G ∈ R N×4×2 is the core tensor and u ℓ 1 i ∈ R N×N , u ℓ 2 j ∈ R 4×4 , and u ℓ 3 k ∈ R 2×2 are singular value matrices and orthogonal matrices. After identifying u ℓ 2 j and u ℓ 3 k associated with properties of interest, we attempt to find G(ℓ 1 ℓ 2 ℓ 3 ) with the largest absolute value, given ℓ 2 , ℓ 3 . Using the selected ℓ 1 , P-values are attributed to the ith by assuming that u ℓ 1 i obeys a Gaussian distribution,

The computed P-values were corrected by the BH criterion, and genes associated with adjusted P-values less than 0.01 were selected.

Integrated analysis of gene expression and m6A profiles. To integrate gene expression and m6A profiles, we employed a recently proposed KTD-based unsupervised FE 26 . We define a tensor x tjt ′ j ′ ∈ R 4×2×4×2 as HOSVD was applied to x tjt ′ j ′ , and we obtained where G ∈ R 4×2×4×2 is the core tensor, and u ℓ 1 t ∈ R 4×4 , u ℓ 2 t ′ ∈ R 2×2 , u ℓ 3 j ∈ R 4×4 , and u ℓ 4 j ′ ∈ R 2×2 are singular value matrices and orthogonal matrices. Here, it should be noted that u ℓ 2 j , u ℓ 3 t ′ , and u ℓ 4 j ′ are attributed to m6A profiles, and only u ℓ 1 t is attributed to the gene expression profiles.

In order to identify genes whose expression profiles depend on time and genomic regions where m6A profiles depend on time associated with the distinction between m6A and control, we need to find which u ℓ 1 t and u ℓ 3 t ′ depend on time and which u ℓ 2 j and u ℓ 4 j ′ are distinct between control and m6A (since x tjt ′ j ′ does not change even if j is replaced with j ′ , u ℓ 2 j = u ℓ 4 j ′ ). Once ℓ 1 , ℓ 2 , ℓ 3 , and ℓ 4 are identified, we can compute the singular value vectors attributed to gene expression samples, u ℓ 1 i , and m6A profiles, u ℓ 2 ℓ 3 k , can be computed as

The P-values are attributed to is and ks as where σ ℓ 1 and σ ℓ 2 ℓ 3 are standard deviations. The computed P-values were corrected by the BH criterion, and genes, i, and genomic regions, k, associated with adjusted P-values less than 0.01 were selected.

Integrated analysis of the 2nd synthetic data set. To integrate x ijk and x ′ ijk in the second synthetic data set, we define a tensor,

After applying HOSVD to x jk 1 k 2 j ′ k ′ 1 k ′ 2 , we obtained (12) www.nature.com/scientificreports/ Singular value vectors attributed to is can be reproduced as After identifying u ℓ 1 j , u ℓ 2 k 1 , u ℓ 4 j , and u ℓ 5 k 1 is considered, P-values are attributed to i, as

The computed P-values were corrected by the BH criterion, and is associated with adjusted P-values less than 0.01 were selected.

After selecting genomic regions, we needed to retrieve the gene symbols included in the selected genomic regions. This could be done using the biomaRt package implemented in R by specifying the hg19 human genome to which short reads were mapped.

Ensembl gene ID to gene symbol. Since gene expression profiles are defined using Ensembl gene IDs, we needed to convert these IDs to gene symbols. This was done by uploading gene symbols selected by TD-based unsupervised FE to DAVID 27 . Uploaded Ensembl gene IDs were converted to gene symbols using the gene ID conversion tool implemented in DAVID by specifying the official gene symbol as the target of conversion. Enrichment analysis. Identified gene symbols were uploaded to Enrichr 28 , which is an enrichment server, to evaluate various enrichments within sets of identified gene symbols.

Various conventional feature selections. Linear regression-based feature selection. To select genes or genomic regions using linear regression analysis, the ls function in the base package in R was used. P-values computed by ls were corrected by the BH criterion, and genes or genomic regions associated with adjusted Pvalues less than 0.01 or 0.05 were selected.

When linear regression was applied to gene expression, x it , was assumed, where T(1) = 0, T(2) = 6, T(3) = 12, and T(4) = 24.

When linear regression was applied to m6A, x ktj , was assumed. When linear regression was applied to the 1st synthetic data, was assumed.

SAM. When SAM 29 was applied to gene expression, x it , or m6A, x ktj , are assumed to be classified into four classes based on t (for gene expression) or eight classes based on the combination of t and j (for m6A), respectively. The sam function was implemented in the siggenes package in R.

Limma. When limma 30 was applied to gene expression, x it , or m6A, x ktj , respectively, they were assumed to be classified in the same way as in SAM. Limma was applied to logarithmically converted x it or x ktj . The limma function was implemented in the limma package in R. When limma was applied to the first synthetic data set, x ijk was classified into eight classes based on the pairs of j and k. Because x ijk takes both positive and negative values, x ijk s themselves were regarded as logarithmically converted valRes. (21) x jk 1 k 2 j ′ k ′

G(ℓ 1 ℓ 2 ℓ 3 ℓ 4 ℓ 5 ℓ 6 ) ×u ℓ 1 j u ℓ 2 k 1 u ℓ 3 k 2 u ℓ 4 j ′ u ℓ 5 k ′ 1 u ℓ 6 k ′ 2 .

x ′ ijk 1 u ℓ 4 j u ℓ 5 k 1 .

(26) x it = a i + b i T(t) (27) x ktj = a k + b k T(t)j www.nature.com/scientificreports/ Random forest. When random forest 31 was applied to gene expression, x it , or m6A, x ktj , respectively, they are assumed to be classified in the same way as in SAM. When it was applied to the first synthetic data set, x ijk was classified into eight classes, as in the case of limma. The randomForest function was implemented in the randomForest package. Features included in OOB were selected by selecting features with non-zero importance given by the importance function implemented in the randomForest package in R.

The pathophysiology of 'happy' hypoxemia in COVID-19

The role of hypoxia in cancer progression, angiogenesis, metastasis, and resistance to therapy

The hypoxic brain

Lung cell hypoxia: Role of mitochondrial reactive oxygen species signaling in triggering responses

Mechanisms of hypoxemia

N6-methyladenosine is required for the hypoxic stabilization of specific mRNAs

Reprogramming of m6a epitranscriptome is crucial for shaping of transcriptome and proteome in response to hypoxia

Tensor factorization toward precision medicine

A survey of some tensor analysis techniques for biological systems

Tightly integrated genomic and epigenomic data mining using tensor decomposition

Tensor decomposition for multiple-tissue gene expression experiments

Tensor decomposition of stimulated monocyte and macrophage gene expression profiles identifies neurodegenerative disease-specific trans-eqtls

Three-way clustering of multi-tissue multi-individual gene expression data using semi-nonnegative tensor decomposition

Classification of clinical gene-sample-time microarray expression data via tensor decomposition methods

Differentially expressed genes extracted by the tensor robust principal component analysis (TRPCA) method. Complexity 1-13

Tensor decomposition for sub-typing of complex diseases based on clinical and genomic data

GSVD-and tensor GSVD-uncovered patterns of DNA copy-number alterations predict adenocarcinomas survival in general and in response to platinum

Hypoxia and mitochondrial oxidative metabolism

vascular smooth muscles and endothelium

Hypoxia-induced inhibition of the endocannabinoid system in glioblastoma cells

Hypoxia-induced signaling activation in neurodegenerative diseases: Targets for new therapeutic strategies

Transcriptional regulation of genes encoding glycolytic enzymes by hypoxiainducible factor 1

Feature selection using principal component analysis

Unsupervised Feature Extraction Applied to

Application of tensor decomposition to gene expression of infection of mouse hepatitis virus can identify critical human genes and efffective drugs for SARS-CoV-2 infection

Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources

Enrichr: A comprehensive gene set enrichment analysis web server 2016 update

Significance analysis of microarrays applied to the ionizing radiation response

limma powers differential expression analyses for RNA-sequencing and microarray studies

Classification and regression by randomforest

This study was supported by KAKENHI 20K12067, 20H04848, and 19H05270.

All authors planned the study. Y.H.T. performed the analyses. All authors validated and discussed the results. All authors reviewed the manuscript.

The authors declare no competing interests.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-87779-7.Correspondence and requests for materials should be addressed to Y.-H.T.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.