ART FounDATion-LOG (ARabidopsis Transcription regulatory Factor Domain-domain interaction Analysis Tool-Liquid liquid phase separation, Oligomerization, Go analysis), a toolkit for interaction data based domain analysis

1. Introduction

Peach (Prunus persica) has been bred for more than 4000 years [1]. Traditional breeding has selected peach cultivars with improved fruit quality and traits over last thousands of years. In the last several decades, marker-assisted breeding has been developed based on advanced next generation sequencing technologies and has gained popularity among breeding scientists [2]. Genome-wide association study (GWAS) have been employed to improve the marker-assisted breeding [2]. However, identifying important functional genetic variants in GWAS data remains challenging due to the complexity of biological systems. There are only limited resources available in peach, compared to model species- Arabidopsis thaliana (A. thaliana). Furthermore, considerable portions of regulatory mechanisms have been conserved across plant species; TF families of A. thaliana are subsets of those of peach. To effectively solve the problem, we took a strategic approach: integration of the immense reservoir of omics data from the model species into the genetic variant analysis in peach; We developed bioinformatics programs based on A. thaliana in order to predict determinative genetic variants using AI and statistics methods.

Well-organized A. thaliana databases can expedite processes of gene regulatory networks (GRN) construction; values of nodes of GRN can be derived from various data such as RNA-seq, transcriptome wide association study (TWAS), epigenome-wide association study (EWAS), phosphorylation site information (proteome), and metabolic pathway in addition to data on genome-wide positions of domain and key regulatory elements [3,4]. However, important pieces of information is missing- interaction modes between regulatory elements/domains, which causes the research bottleneck. For biological processes employ complex regulatory systems, complete annotations of individual interactions are not available. However, integration of LOG, domain-domain interaction (DDI), domain-motif (peptide) interaction (DMI), physicochemical and structural properties of domains, and post-translational modification (PTM) help elucidate mechanisms underlying protein interactions. For example, transcription regulation is triggered by cellular signals and achieved with spatiotemporal coordinations of various protein-nucleic acid interactions (PNI) and protein-protein interactions (PPI) regulated in a number of ways, e.g., protein oligomerization, Liquid liquid phase separation (LLPS) and PTM [5,6]. They have evolved into interwoven regulatory mechanisms in modulating the interactions, for instance, repressing PNI of Histone H1 by inhibiting LLPS formation with phosphorylation, regulating binding activity between scaffold proteins and peptide with oligomerization and LLPS formation, and regulating localization of self-assembled transmembrane protein by phosphorylation of amino acid in domain [7,8,9,10,11,12,13,14].

Here, we introduce ARabidopsis Transcription regulatory Factor Domain-domain interaction Analysis Tool- LLPS, Oligomerization, GO analysis (ART FounDATion-LOG), a usefoul toolkit with integrative resources on key properties modulating PPI and PNI: protein oligomerization, LLPS, protein domain characterization such as structural prerequisites of a particular interaction (e.g. transmembrane, domain linker, intrinsically disordered regions (IDR), DMI, coiled-coil region) and post-translational modification (PTM) (Figure 1) [15,16,17,18,19]. The program consists of seven main modules where each module was built based on existing databases: protein/ protein assembly (prot-assembly) oligomerization module (ProtCAD), features of binding interface of DDI and domain characterization (3did, Plant-PrAS, qPTMplants), LLPS formations (DrLLPS), GO/PO analysis, TF binding profile module (Cis-BP), PPI module (String), TF target module (TF2DNA, Yu et al.), in addition to sequence and structure annotation databases such as UniProt, Pfam, InterPro, PDB, and TAIR [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. Analysis on LOG and domain characterization provides important information on diverse properties of prot assemblies, which help elucidate regulatory mechanisms. Integration of LOG in studying regulation of DDI and domain-nucleotide interaction (DNI) in GRN will significantly enhance prediction power of AI in assessing impact of genetic variants on phenotypic differences. The program is a versatile tool to study a wide spectrum of biological research subjects, for example, gene expression regulation such as TF-target binding activities and prediction on functions of transmembrane proteins within the context of localization.

To demonstrate how to use the program, we performed three simple tests. In the first test, we selected the entries common to oligomerization module and DDI module, divided the data into LLPS and non-LLPS datasets, and developed AI models. The accuracy of the model was 82%, 87%, and 91% in predicting the oligomerization types in LLPS dataset, the LLPS subcellular types, and the oligomerization types in non-LLPS dataset, respectively. In the second test, we selected LLPS factors interacting with TFs and involving in LLPS, merged LOG and TF-target data, and predicted oligomerization types with the accuracy of 97%. Prediction class distribution table was generated with Weka program [37]. In the last test, we selected LLPS factors involved in signaling pathway, made a list of DNA binding motif types of TFs that the factor interacted with, and predicted the motif types with the accuracy of 84%. Association study revealed that DNA binding motif types had strong associations with physicochemical properties, protein secondary structural properties, and important domain features of TFs, while LLPS functional type (scaffold, regulator, client) had strong associations with a category of binding interface of DDI and of DNA binding motif types of TFs. There were no structure data of in-vivo proteins available in PDB. The oligomerization types were predicted based on the prot assemblies grown in-vitro, mainly without LLPS formation. Relatively high associations between categories of DDI and of LLPS functional types suggested that the differences of features of binding interfaces from different domains might reflect the intrinsic properties of domain organization in LLPS formation, which had consistency with correlation analysis in the second test; Domain dimer feature from ProtCAD and DDI interface feature from 3did, LLPS types, LLPS functional types, and number of domains of LLPS factors had relatively high correlation with oligomerization types (corr > 0.6). A considerable portion of regulatory SNPs is related to LOG and domain features [38]. DDI/DNI analysis will provide important clues on what causes phenotypic differences by comparing genetic variants in unknown elements to known regulatory network in the program. The program codes and datasets of ART FounDATion-LOG are available for download at www.artfoundation.kr and sourceforge.

2. Methods

2.1. Database Modules

ART FounDATion-LOG consists of 7 DB modules containing features extracted from existing databases. It also included a simple model for analyzing TF-target interaction and a rough sketch for detecting nucleotide-containing ligand binding motifs in protein based on NBDB [39].

Oligomerization module: ProtCAD provided prot-assembly information derived from PDB entries. PDB contained multiple plausible in-vitro structures of prot assemblies that formed homo- or hetero-oligomers by oligomerizing by themselves or with other proteins. It provided information such as stoichiometries and symmetries of clusters belonging to a ProtCAD entry (GroupID) with the same Pfam architecture [20]. In this paper, we used the term “homo” to refer to a cluster with the same sequence(s) and had only one letter “A” in stoichiometry while “hetero” to refer to a cluster with the different sequences, e.g. “AB”, in stoichiometry. To distinguish those without symmetry, we used C1_obligate_hetero_single_oligomer_obligate to refer to a cluster with only one type: C1 molecules from multiple sequences (e.g. C1-A2BC) and CMA to a cluster having C1-A1. The term “oligomer” was to refer to a cluster with 2 or higher number, e.g. C2 or D3 in symmetry. We extracted 19 variables which were mainly frequencies of domains (single domains and domain dimers). Pfam dimer was to refer to two domains not necessarily maintaining the sequential order in actual sequence; For example a prot-assembly with Dom1 Dom2 Dom3 in sequence had three possible combinations of domain dimers: Dom1Dom2; Dom1Dom3; Dom2Dom3. Most of variables were vectors containing frequencies in 14 or 11 different oligomerization types, for example, those from an entry with one Pfam assignment or with multiple Pfams, from molecules of the symmetry with an even number (e.g. C2, D2) or an odd number (C3, D3), from a single gene or multiple genes, or from 16 different subcellular locations. Two variables contained a single digit which was the maximum number of symmetries (Table S1).

2.2. DDI Module

3did provided sequences of interfaces from DDI and from DMI [21]. We grouped entries of the same domain members into two clusters, one with the non-redundant (nr) sets (e.g. Dom1Dom2) and the other with redundant sets (e.g. Dom1Dom2Dom1). We extracted 13 features to measure the differences between two clusters. They mainly represented the chemical properties such as the number of interacting domains, the number of interacting motifs, the number of amino acids in each fragment where a fragment indicating sub-region of consecutive amino acids without gaps larger than 3, and ProtCAD value- maximum of symmetries. They also contained information such as mean and standard deviation of values in each cluster as well as a size of memberships of cluster. The ANOVA test was performed regardless of the normality of data because the magnitude of the differences between two clusters was particularly important, but no non-parametric statistics with the capacity was available.

Plant-PrAS provided genome-wide analysis of proteins: the grand average of hydrophobicity (GRAVY), isoelectric point (pI), binary presence/absence values of solubility, low complexity, protein secondary structural properties (b-sheets, IDRs, signal peptide(s), transmembrane helices, disulfide (S–S) bonds and domain linkers), N/O-glycosylation sites, ubiquitination sites, functional regions (PASS), peptide type (chloroplast transit peptide, mitochondrial targeting peptide, secretory pathway signal peptide), and subcellular location (E.R., chlo, mito, cysk, cyto, nucl, plas, extr, golg, pero, and vacu). The domain features derived from Plant-PrAS database will be referred as “Plant-PrAS features”. qPTMplants provided PTM information: Glycation, Lysine, Methylation, N-glycosylation, N-terminal, O-GlcNAcylation, Oxidation, Persulfidation, Phosphorylation, S-cyanylation, S-nitrosylation, S-sulfenylation. The domain features derived from qPTMplants database will be referred as “PTM features”. We mapped them to domains and inter-domain regions based on protein domain positions provided by TAIR [25].

2.3. LLPS Module

DrLLPS included approx. 40 distinct biomolecular condensates (Balbiani body, Cajal body, Centrosome/Spindle pole body, Chromatin, Chromatoid body, Cleavage body, DDX1 body, DNA damage foci, Droplet, Gemini of cajal body, Germ plasm/Polar granule, Histone locus body, Insulator body, Microtubule, Mitochondrial RNA granule, Neuronal granule, Nuage, Nuclear pore complex, Nuclear speckle, Nuclear stress body, Nucleolus, OPT domain, Others, Paraspeckle, P-body, PcG body, Pericentriolar matrix, Perinucleolar compartment, P granule, PML nuclear body, Postsynaptic density, Pyrenoid matrix, Receptor cluster, Sam68 nuclear body, siRNA body, Spindle apparatus, Sponge body, Stress granule, TAM body, U body) [24]. LLPS associated proteins usually were involved in formation of multiple condensates. In the database, a protein was classified according to the associations with the condensates, which resulted in 265 possible LLPS types (e.g. a protein only specialized in PML body formations, a protein involved in a number of LLPS- nucleoli, nuclear speckles, cajal body, centrosome, etc.). In addition, DrLLPS provided functional types of LLPS proteins: Client, Regulator, and Scaffold; In this paper, they were referred to LLPS factors. We created two different levels of variables; one with gene as unit and the other with domain dimer as unit. In gene level, we extracted entire domains belonging to a protein, and counted frequencies of the domains in LLPS types. In the same veins, we extracted the functional type of protein and repeated the process. Considering that we had prot assemblies in different modules and that some of large LLPS factors might have evolved from multiple genes, we created a variable to include partial match to larger molecules in LLPS types and counted the number of Pfam assignments of the larger molecules. In domain dimer level, we made lists of possible domain dimers, and calculated the frequencies of domain dimers. Mainly, variables were vectors containing the frequencies in different LLPS types. RNA binding domain, DNA binding domain (DBD), DMI from 3did, and domains with low complexity regions, with disordered regions, with repeats, with coiled-coil structures, with phosphorylation sites, and with active sites such as residues responsible for catalysis, which will be referred as “special flags”. They were created based on Pfam, D²P², and DrLLPS [24,33,40].

2.4. GO Analysis Module

TAIR provides GO and PO data [26,27]. We created five categories- GO analysis, signaling pathway, gene association, PO- anatomy gene, and PO- temporal gene. We retrieved 4 types of subcategories; the first one included following attributes involving signaling pathway, for example hormones, response to light, and osmosensing, the second one included 34 major GO terms for GO analysis such as cell communication and response to abiotic stimulus, and the third one included words related to regulatory roles in the annotation, e.g. enhancer, suppressor, chaperone, and activator, and the last sub-category was original attribute of the database, e.g., acts_upstream_of_negative_effect, or part_of. We created vector variables containing frequency information of categories in the same way as those in oligomerization module or LLPS module.

2.5. TF-Target Module

TF information was provided by Cis-BP for human, A. thaliana, and peach [28]. Cis-BP predicted sequence preferences of TFs and measured correlations between DBD sequence similarities and DNA sequence preferences. We counted the number of how many types of DNA motifs TF bound and how many TFs DNA motif bound. The information were incorporated in order to effectively search for TF and TF targets in respect to their relationship to LOG.

2.6. PPI Module

String and TcoF-DB databases provided PPI information for human, A. thaliana, and peach [29,41]. In addition, AtRegNet and Interactome 2.0 provided information for A. thaliana [25,42]. LPInsider and NPInter databases provided interaction between proteins and RNAs [43,44]. The RNAs were grouped according to a type such as lncRNA and miRNA. Pfam assignment into transcript was provided by GenBank, Gencode, TAIR, InterPro, and Pfam [25,33,34,45,46].

2.7. TF to Target Module

TF binding sites in targets were provided by TF2DNA database for human and by Yu et al (2016) for A. thaliana [25,30]. The binding sites were mapped to gene features by bedmap program [47]. Gencode GFF file and Ensembl GFF file were used for human and for A. thaliana, respectively [46,48]. Gene features included CDS, exon, UTR, intron, upstream-, downstream regions and binding frequencies in each feature were counted.

2.8. Proof of Concept of Search Algorithms

We provided a rough sketch to study interactions between DBD and target motif and between nucleotide-containing ligand and protein's ligand binding motif. To reduce dimensions of variables, amino acids were grouped according to polarity and charge of their side chain (Table 1). Cysteine, glycine, histidine, and proline were considered to have special properties. Each of them made a single membership group.

Using new amino acid group letters, the frequencies of trimers in DBDs were generated: PPP, PPN, PNP, ..., RHR, RRH.

For each DBD, Cis-BP provided ambiguous DNA sequences of the binding motifs in target genes. DNAs and ambiguous DNAs were reassigned to DNA group letters (Table 2). Trimers of DNA group letters and their frequencies in binding motifs were generated.

NBDB provided protein motifs (conserved sequence profiles) that interacted with 24 nucleotide-containing ligands (AMP, ADP, ATP, GMP, GDP, GTP, CTP, CoA, Acelyl-CoA, FMN, F-420, FAD(H), NAD(H), NADP, cyclic nucleotides and dinucleotides, cAMP, cGMP, c-di-AMP, c-di-GMP, and other biologically-relevant cofactors (SAM, PPS, PAP, PLP, ThPP, THD) [39]. We converted NBDB member sequences such as ENAGDTEAPT into new amino acid group letters, and created vector variables containing frequencies of the trimers of the new amino acid group letters per member sequence. Combinations of atoms and moieties belonging to 24 ligands were assigned into 11 groups: RBP, RBPF, RBPN, RBPS, RBPSO, RBSO, TOP, OP, TP, RPF, RPFO (R:ribose, B:base, P:phosphate, F:Flavin, N:Nicotinamide, S:Sulfur, and O:other moiety).

2.9. Demonstration of the Program Usage

2.9.1. Predictions of Important Features- LLPS Subcellular Type, LLPS Factor Type, GO Type, Gene Association Type, Signaling Pathway Type, Subcellular Location, Oligomerization Feature, and Oligomerization Type

We selected entries with multiple Pfam IDs, extracted features from LOG modules except those containing information on oligomerization types of domain dimers or of ProtCAD entries, converted each numeric vector variable to a categorical variable by mapping or clustering algorithms- EM, MakeDensityBasedClusterer, and SimpleKMeans in WEKA program. Most of the vector variables were sparse and mapped to categorical variables without applying clustering methods. The clustering algorithms were applied to three vector variables in LLPS module. Cluster memberships were the values of the categorical variables. We created datasets with entries common to oligomerization and to DDI modules, divided the sets into LLPS and non-LLPS based on LLPS module (Dataset_llps, Dataset_llps2, Dataset_none_llps).

2.9.2. Predictions of Oligomerization Type and Correlations Analysis

We selected TFs with multiple Pfam IDs from TF-target modules, retrieved TF interacting proteins from PPI module, selected only those proteins having LLPS properties, and retrieved information of the LLPS factors from LOG and TF-target modules (Dataset_tf_llps_factor). We excluded variables containing information on oligomerization types of the prot-assembly entries.

2.9.3. Prediction of DNA Binding Motif Types and Association Study

We selected target motif types of TFs that interacted with LLPS factors by physical contacts from the second test, retrieved information from LOG modules as the first and the second test, and additionally domain characterization information from DDI module (Dataset_co_tf_pras_ptm).

To datasets in 1, 2, and 3, we applied classification algorithms in WEKA program (Supplementary Data). We applied FAMD with the package ‘FactoMineR’ in R [49]. FAMD outputs contained information on coordinates of data projected in principal dimensions, cos2, and contrib of variables where cos2 was the quality of representation on principal dimension space and contrib was the contribution to principal dimensions. We applied associations function in dython module from Python to calculate Pearson's correlations [50]. In the third test, we used hotspot algorithm in Weka library to conduct association study to identify the variables having strong associations with DNA binding motif types and with LLPS functional type (scaffold, regulator, client). Detailed information on the algorithms used in Weka library is in the Prediction folder of supplementary data.

3. Results

3.1. Prediction of Important Features

AI models had the average accuracy of 87% in predicting the LLPS types on LLPS dataset (Prediction_llps.llps.type.txt). AI model based on the LLPS dataset without three clustered variables had the accuracy of 82% in predicting oligomerization type, which was higher than that based on all the variables by 25% (Prediction_llps.txt). Three oligomerization types- homo_hetero_moderate_oligomer_obligate, homo_obligate_monomer_oligomer_moderate, and homo_obligate_monomer_obligate- belonged to only prot assemblies in LLPS dataset and not to those in non-LLPS dataset. In the entire dataset, there was only one protein with homo_obligate_monomer_obligate type. Formation and dissolution of LLPS are dynamic and correlate with concentrations of proteins and nucleotides. Therefore, it seemed reasonable to have higher occurrences of the prot assemblies which formed both oligomers and monomers- those ending with “monomer_oligomer_moderate”- in LLPS dataset. Physicochemical and structural properties of LLPS remain elusive. If domain rearrangement occurs within LLPS, it might provide an explanation for higher occurrences of prot assemblies starting with "homo_hetero_moderate", which may have domain sets flexible enough to accommodate domain binding interfaces of both intra-molecule and inter-molecules. We had nested variables belonging to intra_inter_freq that was derived from ProtCAD; they indicated the frequencies of domain sets observed within the same gene and across multiple genes. Those starting with homo_hetero_moderate had non-zero numbers in the nested variables. It should be also noted that 1417 prot assemblies were included in LLPS dataset while only 114 in non-LLPS dataset. AI model developed based on non-LLPS dataset had higher accuracy- 95% (Prediction_none_llps.txt, Figure S1).

3.2. Prediction and Extraction of Important Features from TF-LLPS Factor Data

The accuracy in predicting dominant oligomerization type was 97% (Prediction_tf_llps_factor, Figures S2 and S3). In addition, oligomerization type prediction distributions, in other words probabilities of oligomerization types a prot-assembly formed, were calculated by AI algorithms in Weka library; An example is given in Table S2. Caution should be taken in selecting what prot assemblies to be included in a dataset, if you plan to analyze results from oligomerization type prediction distribution. There were several entries that had the same domain membership but different oligomerization types in the dataset. We selected one oligomerization type at random and included it in the dataset for each of these entries. The distribution had the actual oligomerization types in sequential order from the highest probability. The dataset was made on LLPS factors that interacted with TFs by physical contacts based on String database. However, only a few of them had several entries with multiple oligomerization types. Therefore, an estimation on the credibility of the accuracy needs to be made. Correlation analysis showed that domain dimer feature from ProtCAD, features of binding interface of DDI from 3did, flags, LLPS types, LLPS functional types, and number of domains of LLPS factors had relatively high correlation with oligomerization types (corr > 0.6). We applied FAMD and plotted coordinates of the variables in the first and the second principal dimensions; variables of TF's target motifs were located close to the variables from oligomerization and LLPS modules (Figure 2).

Note: Attributes in the vertical column in left is Y axis and those in the horizontal row on top is X axis in each scatter plot. Labels of X-axis- A:protcad_llps_one_pfam_groupid_sing_intra_inter_freq_sum; B:protcad_llps_one_pfam_groupid_sing_oli_stat_sum; C:protcad_llps_one_pfam_groupid_sing_even_odd_sum; D:protcad_llps_one_pfam_groupid_sing_all_max_sum; E:protcad_llps_dimer_sing_oli_sum; F:protcad_llps_dimer_sing_intra_inter_freq_sum; G:protcad_llps_dimer_multi_oli_sum; H:protcad_llps_dimer_multi_intra_inter_freq_sum; I:protcad_llps_one_pfam_sing_intra_inter_freq_sum; J:protcad_llps_one_pfam_multi_intra_inter_freq_sum; K:protcad_llps_one_pfam_sing_oli_stat_sum; L:protcad_llps_one_pfam_multi_oli_stat_sum; M:protcad_llps_cell_loc; N:threedid_llps; O:flag_llps; P:ANA_dimer; Q:ANA_onepfam; R:TEMP_dimer; S:TEMP_onepfam; T:CONDEN_FUNC_LLPS_onepfam; U:llpsOrNot; V:llpsChkUniprot; W:llpsType; X:llpsUniprot; Y:co_tf_fam_val; Z:co_tf_motif_fam_vall; A1:plant_prAS_data_0; B1:plant_prAS_data_1; C1:plant_prAS_list_data_0; D1:plant_prAS_list_data_1; E1:plant_prAS_dom_linker_data_0; F1:plant_prAS_dom_linker_data_1; G1:plant_prAS_dom_linker_list_data_0; H1:plant_prAS_dom_linker_list_data_1; I1:ptmfeature_data_0; J1:ptmfeature_data_1; K1:ptmfeature_list_data_0; L1:ptmfeature_list_data_1; M1:ptmfeature_dom_linker_data_0; N1:ptmfeature_dom_linker_list_data_0; O1:query_plant_prAS_data; P1:query_plant_prAS_list_data; Q1:query_plant_prAS_dom_linker_data; R1:query_plant_prAS_dom_linker_list_data; S1:query_ptmfeature_data; T1:query_ptmfeature_list_data; U1:co_tf_per_tf_val_0; V1:co_tf_per_tf_val_1; W1:co_tf_per_tf_val_2; X1:co_tf_per_tf_val_3; Y1:co_tf_per_tf_val_4; Z1:co_tf_per_tf_val_5; A2:co_tf_per_tf_val_6; B2:co_tf_per_tf_val_7; C2:co_tf_per_tf_val_8; D2:co_tf_per_tf_val_9; E2:co_tf_per_tf_val_10; F2:co_tf_per_tf_val_11; G2:co_tf_per_tf_val_12; H2:co_tf_per_tf_val_13; I2:co_tf_per_tf_val_14; J2:co_tf_per_tf_val_15; K2:co_tf_per_tf_val_16; L2:co_tf_per_tf_val_17; M2:co_tf_per_tf_val_18; N2:co_tf_per_tf_val_19; O2:co_tf_per_tf_val_20; P2:co_tf_per_tf_val_21; Q2:co_tf_per_tf_val_22; R2:co_tf_per_tf_val_23; S2:co_tf_per_tf_val_24; T2:co_tf_per_tf_val_25; U2:co_tf_per_tf_val_26; V2:co_tf_per_tf_val_27; W2:co_tf_per_tf_val_28; X2:co_tf_per_tf_val_29; Y2:co_tf_per_tf_val_30; Z2:co_tf_per_tf_val_31; A3:co_tf_per_tf_val_32; B3:co_tf_per_tf_val_33; C3:co_tf_per_tf_val_34; D3:co_tf_per_tf_val_35; E3:co_tf_per_tf_val_36; F3:co_tf_motif_val; Y-AXIS; co_tf_motif_fam_val.

Explanation of variable names:

one_pfam/onepfam: Per one domain; dimer: Per domain dimer (e.g., Dom1 and Dom2); multi: Per Pfam architechture of GroupID (e.g., Dom1, Dom2, Dom3); intra_inter_freq_sum: sum of frequencies of domain(s) in three different groups of prot-assemblies; the first, the second, and the third is where the number of the Uniprot genes is less than, the same as, and the greater than the stoichiometry provided in each GroupID, respectively; oli_stat_sum: sum of frequencies of occurrences in the oligomerization type (e.g., C1_obligate_monomer_obligate, C1_obligate_hetero_single_oligomer_obligate, homo_obligate_monomer_oligomer_moderate, homo_obligate_oligomer_obligate, hetero_obligate_monomer_oligomer_moderate, hetero_obligate_oligomer_obligate, homo_hetero_moderate_monomer_obligate, homo_hetero_moderate_monomer_oligomer_moderate, homo_hetero_moderate_oligomer_obligate); even_odd_sum: sum of frequencies of occurrences in two different groups of symmetries; the first and the second is odd (e.g., C3) and even (e.g., C4), respectively.

cell_loc: presence/absence in 16 different subcellular locations; threedid_llps: differences between non-redundant and redundant sets created based on 13 variables; flag_llps: frequences of special flags; ANA: frequences of 295 different PO IDs; TEMP: frequencies of 64 different PO IDs; CONDEN_FUNC_LLPS_onepfam: frequencies of LLPS functional types: client, regulator, scaffold; llpsType: frequencies in 265 different LLPS types; llpsChkUniprot: if uniprot genes are associated with LLPS or not in each llpsType; llpsUniprot: the maximum number of Pfam assignment of associated uniprot genes in each llpsType; co_tf_fam_val: frequencies of PPI with TF families (BHLH,BZIP,C2H2_ZF,CSD,E2F,HOMEODOMAIN,HSF,MADF,MYB/SANT,NAC/NAM,SOX,TBP,TCP,TCR/CXC,UNKNOWN); co_tf_motif_fam_vall: frequencies of TF motif families (MS02_2.00,MS10_2.00,MS11_2.00,MS13_2.00,MS18_2.00,MS21_2.00,MS27_2.00,MS28_2.00,MS31_2.00,MS33_2.00,MS42_2.00,MS46_2.00,MS56_2.00,MS57_2.00,MS59_2.00,MS62_2.00,MS63_2.00,MS64_2.00); For 37 TFs, (AT1G09770,AT1G55520,AT1G75080,AT2G17870,AT2G23380,AT2G23740,AT2G41130,AT3G02150,AT3G12810,AT3G13445,AT3G17609,AT3G19510,AT3G24140,AT3G24520,AT3G28730,AT3G44460,AT3G47620,AT3G48160,AT3G48430,AT3G52300,AT3G56770,AT3G56850,AT4G02020,AT4G02640,AT4G16780,AT4G29000,AT4G34530,AT4G35580,AT4G37790,AT5G04240,AT5G11260,AT5G22220,AT5G22290,AT5G28770,AT5G46690,AT5G51910,AT5G63420), plant_prAS and ptmfeature data were created; prAS_data: values created based on gravy,pi,local,local2,solblity,ubiqui,glycosyl,lowcomp,beta_sheet,disorder,signal,trans_mem,S–S bond,dom_link,pass,o_glycosyl, where loc1 includes chloroplast transit peptide, mitochondrial targeting peptide, secretory pathway signal peptide and loc2 includes E.R.,chlo,mito, cysk, cyto, nucl, plas, extr, golg, pero, vacu; ptmfeature_data: values created based on Glycation,Lysine,Methylation,N-glycosylation,N-terminal,O-GlcNAcylation,Oxidation,Persulfidation,Phosphorylation,S-cyanylation,S-nitrosylation,S-sulfenylation; query: values of LLPS factor; co_tf_per_tf_val: Per each of 37 TFs, domain associated values were calculated. For details, please refer to the codes provided (Table S1).

3.3. Prediction and Association Study on LLPS Factor Interacting with TFs

The accuracy in predicting DNA binding motif types of TFs that the LLPS factor interacted with was 97% (Prediction_co_tf_pras_ptm, Figure 2). Association study with Hotspot algorithm showed that DNA binding motif types had strong associations with Plant-PrAS features of TFs, which is a higher ranked variable created based on all of the Plant-PrAS features. Those serving as both scaffold and regulator had coiled-coil regions. All of the presence/absence value of multiple oligomerization type had a perfect match with binary variable- DMI flag. It should be noted that some LLPS factors had the same domain architectures, which might have biased the data.

4. Discussions

The Plant-PrAS features had associations with properties of different functional types of LLPS. For example, LLPS factors that contained WD40 domains served as scaffold. Some of their partner TFs had domain linker, S–S bond, IDR, beta-sheet, low complexity regions, glycosylation, ubiquitination, all of which also belonged to the LLPS factors (scaffold) themselves, except domain linker. Interestingly, TFs interacting with the LLPS factors all had phosphorylation. TFs containing WD40 involved in transcription activation of anthocyanin synthesis related structural genes in barley [51]. A number of different LLPS scaffolds and regulators seemed to manage coordinated interactions for anthocyanin synthesis, transport, and storage where Natural Deep Eutectic Solvent (NADES) was speculated to be used as inert solvent, which suggested highly complex regulatory processes [52]. As physicochemical properties of liquid condensates remain largely unknown, in-vitro experiments from liquid condensates may encounter problems of partial information. As interaction of liquid condensate is relatively new research topic, standard methods in molecular biology and downstream analysis may require implementation of new protocols and algorithms. Although current technology may have limitations in providing complete information, it may offer practical information for biomarker development. Structural properties contributing to LLPS formation or satisfying constraints imposed by the LLPS, which pose an impact on DNA binding sites, may be roughly estimated with comparative study using AI models based on a large number of factors indicating cellular processes retrieved from numerous databases. GO analysis showed that proteins containing WD40 are part of histone deacetylase complex, nuclear pore, vesicle coat, ubiquitin ligase complex, preribosome, and spliceosome, and enable the following: DNA-binding transcription factor, histone binding, kinase binding, protein heterodimerization and homodimerization activity, kinase activity, phosphatase regulator, ribosome binding, signaling receptor complex, structural molecule activity and transcription cis-regulatory region in A. thaliana (Table S3). WD40 may form important structural platforms for proteins involving in epigenetic activities, which are strongly correlated with genetic variants leading to diseases and phenotype diversity [53,54,55,56].

The Plant-PrAS features and PTM features had supplemented the limited representational power of flags. For instance, homeobox, bZIP, and TCR TF families included gene members that interacted with LLPS factors and had domains with coiled-coil regions. Coiled-coil regions showed strong associations with DBD domains and domains containg Intero annotation called "activity". The canonical coiled-coil has structural motif- heptad repeat [57]. The domains involving in oligomerization such as leucine zipper, N-terminal of homeobox, and helix-loop-helix (HLH) also contain repeats. Although a considerable portion of TF families has various repeats in oligomerization domains, repeats flag is only equipped with the capacity to detect domains defined as repeats. Addition of flags with detection capacity for domains containing such repeats and half sites will improve the program performance. Considering cellular in-vivo environment is dynamic with a large number of factors constantly changing, it will increase accuracy if we make predictions based on multiple variables from modules rather than a single variable such as oligomerization type (e.g. homo_obligate_monomer_oligomer_moderate) in predicting interaction modes in PPI. The proof of concept of search algorithms in the method section can be implemented in PNI analysis with additions of structural elements of nucleotides such as repeats and G-quadruplex, which remains for further research [58,59,60,61,62].