The transcriptome data used in this study were obtained from the Protein Atlas database (Uhlén et al., 2015), including 1,206 human cell lines (quantile normalized). For the calculation of the four indexes, the SPM and PEM scores were specified using formula 4, 5 without modification. The index τ and H_g score only estimated the overall sample-specific expression per gene. The EE score was the same as the inverse logarithm form of PEM. These three indexes were not included in this study. The expression differential (ED) score was calculated using formula 6. The star coordinate transformation method is as follows: ED = x i − ∑ i = 1 n x i ​ * ​ s i / ∑ i = 1 n s i (6)

The star coordinate (START) method is based on parallel coordinates (Heinrich and Weiskopf, 2009) and star coordinates (Kandogan, 2000). This study adopted a similar strategy: first, gene expression data were mapped to a two-dimensional space, and then the degree of differential gene expression was calculated (Supplementary Figure S1). The arrangement order of each coordinate axis was determined by hierarchical clustering method. The START method defines a two-dimensional vector A = [A1, A2, A3, … , An] that starts at the origin (Ox, Oy) and extends along the coordinate axis. Let the initial value of the data points be D = (D1, D2, D3, … , Dn) and set to the lengths of the two-dimensional vectors Ai. The Cartesian coordinates of the two-dimensional vector Ai corresponding to Di are (Di * cos Ɵ, Di * sin Ɵ), where Ɵ is the angle of coordinates. The final Cartesian coordinate value of data point D can be obtained from the sum of Ai, as specified in formula 7: P x , y = O x + ∑ i = 1 n D i ∙ cos ⁡ θ , O y + ∑ i = 1 n D i ∙ sin ⁡ θ (7)

The distance from point V, which is the vertical intersection point for point P to the coordinate axis of the target sample to the origin site O is a perfect measure of the degree of gene expression difference (Supplemental document Figure 1A). Let δ be the angle of OP, the length of VO can be calculated by formula 8: VO = P x 2 ​ * ​ P y 2 2 ​ * ​ ⁡ cos θ − δ (8)

In our previous study, we developed a genome-wide transcription factor binding site prediction software, Grit (Huang et al., 2022), based on comparative genomics and the mixed Student’s t-test. In this study, we revised the calculation method for the raw binding site score (RS), specified in formula 9, 10. This revised score was denoted as “Jindex”. The reason for this revision is that, according to Markov Chain theory, the relationship between the probabilities of individual parts should be multiplicative rather than additive. RS = ln 1 M s ∑ L = 1 M s ∏ k = 1 w q k , L k p L k (9) Jindex = Max s ln ∏ k = 1 w q k , L k p L k w ;   1 ≤ s ≤ l − w + 1 (10)

The Jindex calculation represents the maximum of repeated averaging of log likelihood ratios (LLRs). The averaged LLR indicated the possibility for a motif being present at one particular location in a sequence, where w was the width of the motif, L denoted the location being considered, L _k was the nucleotide at position k within this location, p(L _k ) is the background probability of observing nucleotide L _k estimated from the frequency of L _k in that sequence, and q(k, L _k ) is the probability of observing nucleotide L _k estimated from the frequency of the Kth location in the motif. The Jindex for a motif present in a sequence with length l was the maximum of the average of LLRs taken over all locations of s, where M _s was the number of locations in the sequence calculated as l–w + 1. Jindex value >1.0 indicated that the probability of observing the motif in the sequence is higher than the background probability.

A total of 1,575 position weight matrixs (PWMs) for transcription factors were obtained from public sources (Vorontsov et al., 2024; Rauluseviciute et al., 2024). The newest promoter sequence set (550-bp set) was obtained from Jinzer’s website (Yao et al., 2024). Jinzer software was run with PWMs and the promoter sequence set identified candidate transcription factor binding site (TFBS) with Jindex ≥1.0. A gene set was created by assigning target genes to TFs if the gene had at least one TFBS. Each TF–target gene pair was assigned a rank value, which was the Jindex of Jinzer’s output. However, it should be aware that the analysis presented in the manuscript focused on the −1 ∼ −500bp upstream sequence which will lose distal transcription factor binding site located out of this region.

Flaver software (Yao et al., 2024), which implements weighted Kendall’s tau statistics, was used to identify the key transcription factors within each cell line. This method allows for testing the significance of the correlation between the rank orders of genes in the gene set and their corresponding rank orders within the gene list. Let Si and Li, i = 1, … , n be the ranks of the gene set or list, respectively. Furthermore, let (i, Ri), i = 1, … , n be paired ranks, where Ri is a rank entity of L, whose corresponding S has rank i among Sj, j = 1, … , n. Kendall’s tau has the form of formula 11 and the limiting distribution (LD), following the U-statistics of formula 12 approximated to N (0, 1). This method was implemented in Faver software as an STD (-w 0 option): τ = 2 / n n − 1 ∙ ∑ i > j n sgn i − j sgn R i − R j (11) LD = 3 n n − 1 2 2 n + 5 ∙ τ (12)

Shieh (1998) discovered a weighted version of the rank correlation, the weighted Kendall’s tau, which has the form of formula 13: τ w = 2 / ∑ i n v i 2 − ∑ i n v i 2 ∙ ∑ i > j n v i v j sgn i − j sgn R i − R j (13)

The sgn(x) = −1, 0 or 1, if x <, = or >0, and v_i represents the weighting function bounded to [1, n] and ranges from (0, 1). The limiting distribution (LD) can be derived using formula 14: LD w = n τ w 3 lim n → ∞ n − 1 ∑ i n v i 2 lim n → ∞ n − 1 ∑ i n v i 2 (14)

Where when n→∞, LD approximated to N(0, 1).

Four weighting functions were developed in this study, ranging from 0 to 1. The first two weighting functions were either based on the geometric mean of the gene ranks in gene set i_s and gene list (i_l) or separately (formulas 15, 16) to determine how the genes in the top-ranked TF–target gene pairs in the gene set correlated with the top-ranked differentially expressed gene in the gene list. The remaining two weighting functions were either based on the geometric mean of the density of genes in gene set d_s and gene list (d_l) or separately (formulas 17, 18). The source code was deposited in GitHub and is available under a free academic license (https://github.com/thua45/flaver). The v₁ to v₄ methods were implemented in Flaver software as MIX-LINEAR, LINEAR, MIXED-DENSITY-CURVE, and DENSITY-CURVE, and can be specified by options w 1, 2, 7, and 8, respectively. v 1 = 1 − i s − 0.5 n ∙ 1 − i l − 0.5 n (15) v 2 = 1 − i x − 0.5 / n , x = s  or  l (16) v 3 = 1 − d s max d s ∙ 1 − d l max d l (17) v 4 = 1 − d x / max d x , x = s  or  l (18)

HL-60, MOLM-13, OCI-AML-3, and U-937 cell lines were purchased from Procell (Wuhan, China). HL-60, MOLM-13 and U937 cells were cultured in an RPMI 1640 medium supplemented with 10% FBS, 1% glutamine and 1% penicillin–streptomycin (ThermoFisher Scientific, Dreieich, Germany). OCI-AML-3 cells were cultured in a complete medium consisting of 84% IMDM, 15% FBS, and 1% penicillin streptomycin (ThermoFisher Scientific, Dreieich, Germany). All these suspension cell lines were maintained in culture at a density below 1 × 10⁶ cells/mL and were used for seeding in 6-well plates with 1 × 10⁶ cells per well. All cell lines were grown in a humidified air incubator at 37°C containing 5% CO₂. The cells were treated with shRNA lentiviral vectors targeting ZNF460, SPI1, SPIB, ZNF384, ZNF784, and BATF3 (detailed information available in Supplemental document) following manufacturer’s instructions, and cell proliferation was measured using CellTrace™ CFSE Cell Proliferation Kit (ThermoFisher Scientific, Cat NO. C34554).

The Jinzer run took 2 h on a 8-core Dell desktop computer and identified a total of 5.91 million significant TFBS (Jindex ≥1.0). The number of target genes found in at least one TFBS for the PWMs varied from 1,548 to 7,150. TGIF1 and OVOL2 had the highest and lowest number of target genes, respectively. Compared to the 11,539 human H3K27ac Chip-Seq datasets collected from the Chip-Atlas database (Zou et al., 2024), 60.4% of the TFBS candidates were supported by experimental evidence. For example, NFKB1 transcription has 3,879 candidate target genes with at least one TFBS identified by Jinzer. Among these, 68.3% of the targets had a Chip-Seq pick covering the TFBS, in line with our prediction. The Jinzer gene set file was created by assigning a RANK score to each TFtarget gene pair, which was assigned from the Jindex value.

The distribution of Jindex for the transcription factor binding sites identified by Jinzer is shown in Figure 1. The results indicated that the Jindex values of the binding sites for the TF-target gene pairs were approximately normally distributed (Figure 1A). Because binding sites with higher Jindex values are more likely to be true binding sites, we set a threshold value for Jindex in practice; that is, binding sites with Jindex scores greater than the threshold value were considered candidate sites (Figure 1B). These candidate loci can be converted into a ranked gene set and used as an input file for Flaver, denoted as “Jinzer set” in this study.

Transcriptome data from 1,206 human cell lines were transformed into ranked cell-specific gene lists using the four indexes described in the Materials and Methods section. The ranks of the ED (Figure 2A) and PEM (Figure 2B) indexes were similar and approximately normally distributed. The major difference between the distributions of the ranks of the ED and PEM indexes was that the standard deviation of the ED index was significantly higher than that of the PEM index. The results of the START (Figure 2C) were unique, informative, and interesting. The gene list created using the SPM index was severely skewed and distributed (Figure 2D). We performed a gene search according to the gene’s Gene Ontology (GO) annotation and found an average of 22.3% overlap between the GO gene sets annotated with sample-specific functions and the cell-specific gene lists identified by the four indexes. This suggested that many genes with known sample-related functions were also dominantly or recessively expressed in the corresponding transcriptome data and vice versa. The gene lists created by the four indexes, named the CELLINE list, were used as inputs for the Flaver software.

R (version 4.2) was used to generate four synthetic gene sets and gene list data to simulate the gene expression and transcription factor binding site data to the real data distribution properties. The gene set and gene list in simulated dataset 1 were normally distributed random data with covariance (R) range of [-1, 1] (Figure 2A shows the simulated data with R = 0.8). For simulated dataset 2, a threshold condition of RANK >0 was set for both the gene set and gene list, in addition to simulated dataset 1, which was used to simulate filtering high-scoring binding sites and upregulated genes when analyzing transcriptome data (Figure 2B shows the simulated data with R = 0.8). The gene set and gene list in simulation dataset 3 were also normally distributed random data with a covariance range in [-1, 1]. However, in contrast to simulation dataset 2, only a filter condition of RANK gene set greater than zero was set, and the gene list had a full-spectrum distribution (Figure 2C shows the simulation data with R = 0.8). Simulation dataset 4 introduces a covariance gradient in addition to simulation dataset 3. The gene correlation coefficient for genes with RANK >1.5 was set to R, the gene correlation coefficient for genes with 1.0 < RANK ≤1.5 was set to R/2, and the gene correlation coefficient for genes with RANK ≤1.0 was set to zero. Simulation dataset 4 was used to simulate situations in which genes with higher RANK had a higher degree of correlation, and genes with lower RANK had lower or no correlation (Figure 2D shows simulated data with R = 0.8).

The four weighting methods described in the Materials and Methods section were used to test the four sets of simulated data. For simulation data 1, the P-value and the covariance of the gene set and gene list showed a good response relationship, with a “U-shaped” line graph. As shown in Figure 4A, on both sides of covariance = 0, the P-value decreased rapidly with an increase in |covariance|, and no significant difference was observed among the four weighting methods. For simulation data 2, the overall shape of the line graph was different from that of simulation data 1; namely, it was “J-shaped.” Specifically, on the covariance (0, 1) side, the P-value decreased rapidly with an increase in |covariance|, which is consistent with simulation test data 1. By contrast, in the covariance (−1, 0) side, no marked decrease was observed in the P-value with an increase in |covariance|, and the results of the four weighting methods were highly consistent. For simulated data 3, the trend between the obtained P-value and the covariance of the gene set and gene list was similar to that of simulation data 1; namely, it was also “U-shaped,” that is, the P-value decreased rapidly with an increase in |covariance| on both sides of covariance = 0. The steepness of the curve was significantly lower than that of simulation data 1, and slight differences were observed between the four weighting methods. For simulated data 4, the curve between the P-value and the covariance of the gene set and gene list was also “U-shaped.” However, the steepness of the curve was flatter than that of simulated data 3, and the P-values obtained by the four weighting methods were significantly different. Specifically, the W7 method was the steepest, followed by the W1, W2, W8, and W0 methods.

The Flaver software was used to identify key transcription factors in the sample-specific gene lists. The real-world gene list, CELLINE list, described in the Materials and Methods section, and the Jinzer set created in the Jinzer software, were used as input files for Flaver. The MIXED-DENSITY-CURVE function, which has been proven to be the most sensitive method (W7), was used to run the Flaver. A total of 644 transcription factors were identified within a running time of 23 h (FDR <0.05). Figure 5 shows the clustering results according to the sign (Dir) * -log (P-value) value of the transcription factors. Based on the Spearman’s rank correlation coefficient distance, the 1,206 cell lines were divided into nine categories (Figures 5A–I), and the transcription factors were divided into 12 categories (Figures 1–5, 5–12).

The proportions of cell line types in each category are plotted in Figure 6. Brain cancer samples were predominantly located in cluster D, while leukemia, lymphoma, and myeloma samples were mainly located in cluster H. Skin cancer samples were located in cluster A, breast cancer samples in class F, and colorectal cancer samples in clusters G and I. Head and neck cancer samples were mainly located in class F, while kidney cancer samples were distributed in cluster G. Uterine and ovarian cancer samples were distributed in cluster C. Lung cancer samples were found in clusters B, C, E, I, and D. Non-cancerous samples were mainly located in cluster D.

From the clustering analysis results (Figure 5), we identified a set of transcription factors (row 9) that showed the same regulatory behavior in leukemia samples (column H). A significant positive correlation was observed between the RANK values of the gene set and the gene list for these transcription factors. This set of transcription factors could be used to distinguish leukemia cell lines from other cell types (Figure 7A). Based on the comparison of the differences in -Sign (P-value) * Log10(P-value) for these transcription factors between the leukemia cell lines and other types of cell lines, we selected six transcription factors with the highest differences between the leukemia and non-cancer groups for further experimental verification, namely ZNF460, SPI1, SPIB, ZNF384, ZNF784, and BATF3 (Figure 7B).

Four cell lines (HL-60, MOLM-13, OCI-AML-3, and U-937) were used for experimental validation. The results which were presented in Figure 8 showed that the ZNF460 shRNA lentiviral vector significantly inhibited the proliferation of HL-60, MOLM-13, OCI-AML-3, and U-937 cells, and that the ZNF460, SPI1, ZNF784, and BATF3 shRNA lentiviral vectors exerted strong inhibitory effects. The inhibition rates of ZNF460, SPIB, ZNF784, and BATF3 in MOLM-13 cells were over 20%. Only one transcription factor (ZNF384) lentiviral vector did not reach significance in the OCI-AML-3 cell lines. The ZNF460 lentivirus vector inhibited the proliferation of 80% of the U-937 cells.