2.1 Identification of DNase I hypersensitive sites (DHSs) from DNase-seq data

   Reads were mapped to unmasked genomes using Bowtie2 (v2.2.8) with ≤ 1 mismatch in seed alignment allowed. Only uniquely mapped reads were kept for peak calling. The hotspot regions and peaks of DNase I hypersensitive sites were called using Hotspot (v4.1.0) under default settings (threshold: FDR 0.01).

2.2 Identification of genomic TF footprints from DNase-seq data

   Reads were mapped as mentioned above. The genomic TF footprints were identified based on the pipeline and tool (fp2012) described in the paper.

2.3 Identification of peaks of TF binding and histone modification from ChIP-seq data

   Reads were mapped as mentioned above. Macs2 (v2.1.1) was employed to call peaks from ChIP-seq data with options of "--keep-dup all". For histone modification, "--broad" was added to merge closed peak regions. The top 600 peaks of TF binding in each sample were extracted to discover TF binding motifs using MEME-ChIP with motif length spanning from 7bp to 20bp.

2.4 Mapping nucleosome positions from MNase-seq data

   Reads were mapped as mentioned above. DANPOS (v2.2.2) was employed to search for nucleosome peaks with pair-end mode (-m 1).

2.5 Calculation of sequence conservation

   Genome sequences and gene annotation data for 135 species collected in PlantRegMap were extracted and genomes with low quality assembly (less than 75% of bases in sequences longer than 10 kb) were removed and species were grouped by taxonomy. The groups with less than 4 species were removed, which remained 63 species in 7 groups for conservation analysis. For each species, RepeatMasker (v4.0.6) with RepBase (20150807) was used to mask transposons and simple repeats. Then pairwise alignments were performed in each group using lastz (v1.03.73) with the parameters of "--inner=2000 --xdrop=9400 -gappedthresh=3000 --hspthresh=2200". Overlapped alignment blocks were joined using axtChain and gaps were filled using chainNet in UCSC Kent Utilities, which produced 322 pairs of alignments with single coverage in target sequences. Multiple alignments were performed based on phylogenetic trees using multiz (v11.2). Non-conserved models were estimated based on the four-fold degenerate sites and conservation (conserved elements, PhastCons and PhyloP scores) was calculated using RPHAST (v1.6). For each species the parameters were tuned until the ratio of conserved elements in aligned CDS regions was near 65% according to the manual. Conserved elements and PhastCons scores were calculated with "viterbi" mode and PhyloP scores were calculated with "CONACC" mode.

2.6 Identification of TFBS with transcriptional regulatory functions (FunTFBS)

   Given the genomic position of specific TFBS, the motif frequencies of each base pair are extracted according to genomic sequence and the binding motif of corresponding TF. The PhyloP scores in the same region are extracted with the same order with motif frequencies (the PhyloP scores will be reversed for binding in negative strands). Pearson correlation is calculated between motif base frequencies and the absolute value of PhyloP scores. If the correlation test is significant (p ≤ 0.05) with correlation score higher than 0.5, the corresponding binding site will be treated as functional TFBS. (see details)

2.7 Collection and inference of transcriptional regulatory interactions

   We collected and inferred regulatory interactions in the following ways:

• Literature mining and manual curation: In our previous study, we mined and manually curated 1 431 functional confirmed transcription regulatory interactions from literature (ATRM), this high-confidence resource were integrated now.
• Identifying regulatory interactions from ChIP-seq: A potential interaction was assigned between a TF and a gene if there was at least one peak located in the promoter region (TSS +500bp ~ -100bp) of the gene. Only the potential interactions verified in ≥ 2 replications were kept.
• Motif: FIMO was used to scan the promoter region (TSS +500bp ~ -100bp) of a genes for binding sites using the high-quality binding motifs of a TF with threshold 1e-5, and a potential interaction was assigned if there was at least one binding site in the promoter of the gene.
• Motif + footprint (DGF): An interaction was assigned only if more than 60% length of a binding site overlapping with the footprints identified above.
• Motif + conserved elements (CE): An interaction was assigned only if more than 50% length of a binding site overlapping with the conserved elements identified above.
• FunTFBS: An interaction was assigned only if there was at least one functional binding site (identified by the FunTFBS method) in the promoter of the gene.

3.1 Binding site prediction

   Using the sets of high-quality, non-redundant binding motifs of TFs for 156 species with whole genome sequences (see details), FIMO is employed to scan the input sequences for TF binding sites.

3.2 Regulation prediction

   This tool is used to infer potential regulatory interactions between TF and input genes, and finds the TFs which possess over-represented targets in the input gene set. As described above, the sets of high-quality binding motifs of TFs and FIMO are used to scan TF binding sites in the promoters, and an interaction is assigned if there are one or more binding sites of a TF in the promoter of a gene. The genome assembly and annotation version described in the data source are used to extract promoter sequences when the coordinates of promoters/genes or gene sequences are input. The regulations predicted from the promoter regions (defined by users) of all annotated genes will be used as the background for enrichment analysis. In addition, refined regulation prediction using previous identified footprints is supported now.

3.3 GO enrichment

   Based on the sets of non-redundancy GO annotation for 165 species generated above, topGO (v2.22) is employed and Fisher’s exact tests are performed to find the significant enriched GO terms for the input gene set. Genome annotation IDs, UniProt AC/ID, Entrez Gene IDs and symbols are supported, and a built-in ID mapping tools is used to mapping the input IDs to genome annotation IDs (see id mapping file here). In the ontology structure plot, rectangles indicate the significant terms. Rectangle color represents the relative significance, ranging from dark red (most significant) to bright yellow (least significant).

3.4 TF enrichment

   Based on pre-generated regulatory interactions integrated from literature mining, ChIP-seq, motif combined footprint and motif mentioned above, Fisher’s exact tests are performed to find the TFs possess significantly over-presented target genes in the input gene set.

3.5 FunTFBS

   This package is designed for users to filter for functional TFBS using the FunTFBS method in any interesting species with custom TF binding motifs and conservation data.