DOI: 10.5061/dryad.9w0vt4bmk
title: Data from: A comparative analysis of stably expressed genes across diverse angiosperms exposes flexibility in underlying promoter architecture

description: Promoters regulate both the amplitude and pattern of gene expression—key
 factors needed for optimization of many synthetic biology applications.
 Previous work in Arabidopsis found that promoters that contain a TATA-box
 element tend to be expressed only under specific conditions or in
 particular tissues, while promoters which lack any known promoter
 elements, thus designated as Coreless, tend to be expressed more
 ubiquitously. To test whether this trend represents a conserved promoter
 design rule, we identified stably expressed genes across multiple
 angiosperm species using publicly available RNA-seq data. Comparisons
 between core promoter architectures and gene expression stability revealed
 differences in core promoter usage in monocots and eudicots. Furthermore,
 when tracing the evolution of a given promoter across species, we found
 that core promoter type was not a strong predictor of expression
 stability. Our analysis suggests that core promoter types are correlative
 rather than causative in promoter expression patterns and highlights the
 challenges in finding or building constitutive promoters that will work
 across diverse plant species.
description: <em>RNA-seq dataset
 processing</em> <em>(Relevant files:
 0_Slurm_Pipeline)</em> RNA-seq atlases were located in the NCBI Sequence
 Read Archive (SRA) database. The references for the datasets can be found
 in Supplemental Table S1. The individual datasets were retrieved using
 sratoolkit-3.0.1 prefetch followed by fasterq-dump functions.
 Fastqc-0.11.9 were used to generate a QC report for each dataset.
 Trimmomatic-0.39 were used for adaptor and low quality ends trimming using
 the following settings: ‘SLIDINGWINDOW:4:20 MINLEN:36’. ILLUMINACLIP files
 TruSEq3-PE-2.fa was supplied for paired end data and TruSEq3-SE.fa were
 supplied for single end data. Reference transcriptome were downloaded from
 the Ensembl Plants (http://plants.ensembl.org/index.html) for <em>Arabidopsis thaliana</em>, <em>Camelina sativa</em>, <em>Cucumis melo</em>, <em>Glycine max</em>, <em>Phaseolus vulgaris</em>,<em> Pisum sativum</em>, <em>Vigna unguiculata</em>, <em>Sorghum bicolor</em>, <em>Zea mays</em>, <em>Solanum lycopersicum</em>, <em>Actinidia chinensis</em>, <em>Triticum aestivum</em> and Phytozome (https://phytozome-next.jgi.doe.gov) for <em>Arachis hypogaea</em>, <em>Cicer arietinum</em>, and <em>Solanum tuberosum</em> (Cunningham et al., 2021; Goodstein et al., 2012). An index file was generated and the reads aligned and counted using Kallisto-0.44.0 with ‘-o counts -b 500’. For single end data, Fragment Length and Standard Deviation were required, but the information is difficult to locate, and so a default value of ‘-l 200 -s 20’ were used across the board. Another Fastqc was performed on the trimmed files, and a final MultiQC-1.13 were run on the entire folder encompassing all the log files that Fastqc, Trimmomatic, and Kallisto generated. The MultiQC report was inspected to ensure the trimming step improved read quality and there were no major warnings. <em>Normalizing count, Calculating CV and Percent Ranking</em> <em>(Relevant files: 1_Metadata_from_RUNselector.Rmd, 2_MOR_Normalization.Rmd)</em> Using an R script, the raw counts for each species were normalized using the DESeq2 package using a metadata file curated from the original study for the RNA-seq datasets. The coefficient of variation across all samples for a given atlas was used as a metric for stability for each gene, and the percentile ranking for each gene was calculated. The geometric mean for each gene was also calculated across all samples.  <em>Extracting intergenic region and 5’UTR</em> <em>(Relevant files: 3_ExtractPromUTR(ALL_Transcripts).ipynb, 8_ExtractPromUTR(Orthologs).ipynb)</em> Gff3 annotation files and reference genomes were downloaded from Ensembl or Phytozome depending on where the reference transcriptomes were retrieved from. 40% of transcripts were selected from the total transcriptome and their intergenic region and 5’UTR were extracted from the Gff3 annotation. Intergenic region and 5’UTRs of identified orthologs were extracted in a similar manner. <em>Labeling core promoter types</em> <em>(Relevant files: 4_Label_Promoters.Rmd, 9_Motif_Scan.Rmd, 10_Octamer_Scan.ipynb)</em> Motif Scan: Intergenic regions and 5’UTR sequences are trimmed to only regions to be scanned for each core promoter types: TATA box (-100 to TSS), Y patch (-100 to +100), and Inr (-10 to +10).  Intergenic regions shorter than 100bps were excluded from analysis. Each regions were scanned for their respective motifs according using motif files as well as methods outlined in (Jores et al., 2021). A motif is considered to be present when the relative motif scores are above 0.85.  Octamer Scan: Intergenic regions and 5’UTR sequences were trimmed based on the positions relative to the TSS outlined in Yamamoto et al. 2009 (TATA, −45 to −18; Y Patch, −50 to +50; CA, −35 to −1; GA, −35 to +75). Each region was scanned for the presence of octamer motifs from the TATA, Y patch, GA, and CA lists outlined in Yamamoto et al. 2009. If the specified region contained at least one motif for a given promoter type, it was labeled as positive. <em>Ortholog Analysis</em> <em>(Relevant files: 5_At_gene_ranking.Rmd, 6_Identifying_orthologs.Rmd, 7_Processing_orthologs.Rmd)</em> The <em>Arabidopsis</em> transcriptome was filtered to only include primary transcripts, and mitochondria as well as chloroplast transcripts were removed. Top 5% stable genes by CV, bottom 5% stable genes by CV and a random set of 1343 genes (5%) were randomly selected. Using biomaRt in R, the Ensembl and Phytozome databases were queried for orthologs for the selected set of <em>Arabdiopsis</em> genes for each species (Durinck et al., 2009). Orthologs from <em>Arachis hypogaea</em>, <em>Cicer arietinum</em>, and <em>Solanum tuberosum</em> were retrieved from Phytozome, and the rest of the species from Ensembl. For analysis in Figure3B, significance test of done by ANOVA followed by Tukey’s HSD. For each target gene that matched to an <em>Arabidopsis</em> transcript, only the highest expressing transcript was kept. If an <em>Arabidopsis</em> transcript retrieved more than one orthologs from a target species, these pairs of orthologs were removed from analysis. We only kept orthologous gene groups that had a “change” in expression pattern, defined as crossing the 50<sup>th</sup> percentile CV, in two target species, and the remaining candidates were manually mapped onto the phylogenetic tree to identify gene groups that had changes in expression pattern that are consistent with the tree. This means having changes in expression pattern that are mostly found in the same clade. Gene trees were built for these candidates using blast-align-tree (https://github.com/steinbrennerlab/blast-align-tree) and the candidate lists were further trimmed based on the gene trees to ensure a 1:1 relationship between all members in the gene group. The dataset contains all the necessary scripts to transform the data as described in the manuscript and perform the analysis in the paper.
description: Scripts are written in R and Python.
description: # A comparative analysis of stably expressed genes across diverse
 angiosperms exposes flexibility in underlying promoter architecture ***
 The dataset contains raw counts from RNA-seq atlases (RNA-seq datasets on
 multiple tissues throughout development) from fifteen species: Arabidopsis
 thaliana, Camelina sativa, Cucumis melo, Glycine max, Phaseolus vulgaris,
 Pisum sativum, Vigna unguiculata, Sorghum bicolor, Zea mays, Solanum
 lycopersicum, Actinidia chinensis, Triticum aestivum, Arachis hypogaea,
 Cicer arietinum, and Solanum tuberosum. Expression stability (coefficient
 of variation) and mean expression strength (geometric mean) were
 calculated for each transcript to determine expression pattern. Core
 promoter elements were scanned for each transcripts to determine promoter
 core type. Relationships between core promoter elements and expression
 pattern were explored. We were able to map gene expression pattern onto
 core promoter type in multiple genomic contexts. We found that while
 TATA-box-containing promoters are consistently over-represented in
 selectively expressed genes, Coreless promoters are only associated with
 stable expression in eudicots but not monocots. Additionally, by
 identifying orthologous gene groups within these species, we were able to
 track changes in core promoter type and expression pattern for groups of
 evolutionarily related promoters. We found that stably expressed genes are
 also more likely to have orthologs in other species compared to unstably
 expressed genes, and the orthologs tend to retain similar expression
 patterns. Lastly, we show that changes in core promoter types do not
 explain changes in expression pattern. This evolution-guided approach
 reveals design rules surrounding core promoter architecture and expression
 patterns. ## Description of data and file structure There are three main
 folders: ``` - Data: Contains all necessay input data to run the scripts
 as well as the output data from the scripts. Files in the folders are
 reference by relative path in the scripts. - Summary.csv: unformatted
 version of Supplemental_Table_S3 - Raw_Counts: Output from Kallisto broken
 down by species - Promoters: - Y_40pcent_Z: Subset of promoters used to
 generate Figure2 - Y_orthologs_Z: Subset of promoters used to generate
 Figure4 - X_Y_annotation: details and genomic position of the promoter
 regions after parsing the gff3 - X_Y_sequence: actual sequence of promoter
 region - Orthologs: Intermediate files generated for ortholog related
 analysis - Normalized_Counts: Raw_Counts processed by DESeq2 - Metadata:
 Metadata regarding experimental setup for the RNAseq datasets. Used to
 perform DESeq2 - Genomes: Reference transcriptome and gene annotation
 (gff3) organized by species - !!! Please note: due to licencing reasons,
 files in this folder is uploaded as supplemental information
 (Data_Genome.zip) on Dryad. Please download the file and unzip the folders
 into Data/Genomes/ and the scripts should work as intended !!! -
 Core_Promoters: Analysis regarding core promoter types -
 40pcent_promoters: Analysis for Figure2 and Supplemental_Figure_S1 -
 X_MotifScores: numerical value of predicted motif scores for each promoter
 - X_MotifsLabeled: Binary True/False for each motif based on motif scores
 with CV and geometric mean information - Orthologs: Final Candidates -
 Motif_results.csv: Motif scoring result (numeric value) -
 Octamer_results.csv: Octamer presence or absence (True/False) - Figures:
 Output folder for figures generated from the analysis. - Scripts: All the
 scripts necessay to perform the analysis. The datasets are numbered by the
 sequence they need to be run to recreate the analysis, and each script
 explicitly refers to necessay input files in the Data folder. -
 0_Slurm_Pipeline are slurm/shell scripts performed on UW's Hyak
 computing cluster, taking in NCBI SRA numbers and generating transcript
 count tables ``` ### Details for 0\_Slurm\_Pipeline This folder contains
 slurm/shell scripts ``` - Each step (with the exception of
 5_kallisto_index.slurm) is performed in parallel - Output files are saved
 under folder [Species_name] - The scripts always take three arguments: 1.
 path to list_of_SRA.txt 2. SE (single end) or PE (paired ends) 3.
 Species_name (For the purpose of directing to the correct input/output
 folder) Not all the arguments are used for all scripts but they're
 always included for consistency sake - !!!Before running, the scripts will
 have to be modified to direct slurm to your folder of choice as they are
 coded currently with absolute path!!! ``` * 1\_prefetch\_parallel.slurm:
 Prefetch data given a list of NCBI SRA numbers * modules: sratoolkit 3.0.1
 * usage example: 1\_prefetch\_parallel\.slurm \[list\_of\_SRA\.txt\] \[SE
 or PE\] \[Species\_name\] * 2\_fasterq-dump\_parallel.slurm: Download
 fasta files given a list of NCBI SRA numbers * modules: sratoolkit 3.0.1 *
 usage example: 2\_fasterq\-dump\_parallel\.slurm \[list\_of\_SRA\.txt\]
 \[SE or PE\] \[Species\_name\] * 3\_fastqc\_parallel.slurm: fastqc on
 downloaded fasta files * modules: java 17.0.5, fastqc 0.11.9 * usage
 example: 3\_fastqc\_parallel\.slurm \[list\_of\_SRA\.txt\] \[SE or PE\]
 \[Species\_name\] * 4\_trimmomatic\_parallel.slurm: performs trimmomatic
 QC step * modules: java 17.0.5, trimmomatic 0.39 * usage example:
 4\_trimmomatic\_parallel\.slurm \[list\_of\_SRA\.txt\] \[SE or PE\]
 \[Species\_name\] * 5\_kallisto\_index.slurm: Kallisto requires building
 an index from the reference transcriptome * modules: kallisto 0.44.0 *
 usage example: The script has to be modified to direct slurm to the folder
 output\_species\_name/reference\_transcriptome\.fa\.gz\, and the script
 can be run by calling 5\_kallisto\_index\.slurm \[list\_of\_SRA\.txt\]
 \[SE or PE\] \[Species\_name\] * 6\_kallisto\_parallel.slurm: running
 Kallisto * modules: kallisto 0.44.0 * usage example:
 6\_kallisto\_parallel\.slurm \[list\_of\_SRA\.txt\] \[SE or PE\]
 \[Species\_name\] * 7\_consolidate\_parallel.slurm: Copies individual
 count files (abundance.tsv) into the same folder (Counts) * modules: None
 * usage example: 7\_consolidate\_parallel\.slurm \[list\_of\_SRA\.txt\]
 \[SE or PE\] \[Species\_name\] * 8\_fastqc\_parallel\_postTrim.slurm:
 Another fastqc after the trimming step to ensure the trim was successful *
 modules: java 17.0.5, fastqc 0.11.9 * usage example:
 8\_fastqc\_parallel\_postTrim\.slurm \[list\_of\_SRA\.txt\] \[SE or PE\]
 \[Species\_name\] ### Details for 1\_Metadata\_from\_RUNselector.Rmd R
 packages: ``` - dplyr 1.1.2 ``` Purpose: * Generate DESeq2 compatible
 metadata files Inputs: ``` - Run-selector information from NCBI database
 or additional information gathered from the original study - Additional
 information can be located in ../Data/Metadata/ -
 ../Data/Metadata/[Species_name].txt ``` Outputs: ``` - DEseq2 compatible
 metadata files - ../Data/Metadata/[Species_name].csv ``` ### Details for
 2\_MOR\_Normalization.Rmd R packages: ``` - tidyverse 2.0.0 - DESeq2
 1.38.3 ``` Purpose: * Takes output from Kallisto (in ../Data/Raw\_Counts)
 and normalize according to the experimental setup (Metadata) using
 DESeq2's Normalization function. Inputs: ``` - Run-selector
 information from NCBI database or additional information gathered from the
 original study - Additional information can be located in
 ../Data/Metadata/ - ../Data/Metadata/[Species_name].txt ``` Outputs: ``` -
 DEseq2 compatible metadata files - ../Data/Metadata/[Species_name].csv ```
 ### Details for 3\_ExtractPromUTR(ALL\_Transcripts).ipynb Python packages:
 ``` - gffpandas.gffpandas 0.1.0 - pandas 2.0.0 ``` Purpose: * This samples
 40% of the total transcriptome of a species and first extract the
 intergenic region and 5'UTR annotation according to the associated
 gff3 file, and then use the annotation to extract actual sequences from
 the reference genome. Parsing GFF3 annotation: Inputs: ``` - List of
 transcripts (../Data/Normalized_Counts/[Genus_species_Summary].csv) - Gff3
 annotation file (../Data/Genomes/[Genus_species]/*.gff3) ``` Outputs: ```
 - Position information for intergenic region and 5'UTR
 (../Data/Promoters/[Species]_40pcent_annotation.csv) ``` Extracting
 Sequences: Inputs: ``` - Position information for intergenic region and
 5'UTR (../Data/Promoters/[Species]_40pcent_annotation.csv) -
 Reference genomes (../Data/Genomes/[Genus_species]/*.fa) ``` Outputs: ```
 - Extracted Sequences
 ("../Data/Promoters/[Species]_40pcent_sequence.csv") ``` Note:
 the codes uses pd.str.contains() to search for seach gene in a loop, which
 isn't terribly efficient once the file gets big and might take hours
 to run. Parsing the Attributes might be a better way to go. ### Details
 for 4\_Label\_Promoters.Rmd R packages: ``` - tidyverse 2.0.0 -
 universalmotif 1.16.0 - ggplot2 3.4.2 - showtext 0.9.6 ``` Purpose: *
 Label the "promoter" regions of the 40% extracted genes as TATA,
 Inr, Ypatch, or Coreless Input files: ``` - Extracted promoter sequences
 stored in ../Data/Promoters/[Genus_species]_40pcent_sequence.csv -
 CPEs.meme files for the motifs stored in ../Data/Core_Promoters/ (obtained
 from Jores et al.) - maize_TATA.meme is generated by modifying the meme
 file above and changing the TATA motif to a maize specific motif published
 in (Mejía-Guerra et al. 2015) ``` Output files: ``` - Scores from motif
 scanning are stored in
 ../Data/Core_Promoters/40pcent_promoters/[Genus_species]_MotifScores.csv -
 Scores are summarized to presence or absence of motif and merged with CV,
 geom_mean data are stored in
 ../Data/Core_Promoters/40pcent_promoters/[Genus_species]_MotifsLabeled.csv
 - Figure2 is saved in ../Figures/Figure2.svg - Supplemental_Figure_S1_A is
 saved in ../Figures/Supplemental_Figure_S1_A.svg - Supplemental_Figure_S2
 is saved in ../Figures/Supplemental_Figure_S2.svg -
 Supplemental_Figure_S1_B is saved in
 ../Figures/Supplemental_Figure_S1_B.svg ``` ### Details for
 5\_At\_gene\_ranking.Rmd R packages: ``` - tidyverse 2.0.0 ``` Purpose: *
 Extract the top 5% stable, bottom 5% stable, and 5% random Arabidopsis
 genes for homology analysis Inputs: ``` -
 ../Data/Normalized_Counts/Arabidopsis_thaliana_Summary.csv (Generated from
 2_MOR_Normalization.Rmd) ``` Output: ``` -
 ../Data/Orthologs/arabidopsis_gene_set.csv ``` ### Details for
 6\_Identifying\_orthologs.Rmd R packages: ``` - biomaRt 2.54.1 - tidyverse
 2.0.0 ``` Purpose: * Takes the top 5% stable, 5% unstable, and 5% random
 genes from Arabidopsis and fetch their orthologs via biomaRt inputs: ``` -
 Arabidopsis gene sets generated from step 5_At_gene_ranking.Rmd
 (../Data/Orthologs/arabidopsis_gene_set.csv) ``` outputs: ``` - Raw output
 from Phytozome is saved in ../Data/Orthologs/Phytozome_orthologs.csv -
 Entries fetched from different species versions are collapsed as the same
 species in ../Data/Orthologs/Phytozome_orthologs_condensed.csv - Output
 from Ensembl is saved in ../Data/Orthologs/Ensembl_orthologs.csv - Merged
 orthologs with arabidopsis CV and geom_mean information is stored in
 ../Data/Orthologs/All_orthologs.csv - ../Figures/Figure3.svg ``` ###
 Details for 7\_Processing\_orthologs.Rmd R packages: ``` - tidyverse 2.0.0
 ``` Purpose: * Taking the ortholog list, filtering down the interesting
 candidates with flipped expression patterns in at least two species.
 inputs: ``` - ../Data/Orthologs/All_orthologs.csv generated from
 6_Identifying_orthologs.Rmd ``` outputs: ``` - Candidates that 1. Changed
 expression pattern in at least 2 species 2. Had changes that mapped onto
 phylogenetic tree in a reasonable manner were selected for and written to
 ../Data/Orthologs/High_Confidence_Candidates.csv - Candidates that also
 had been verified by blast-align tree pipeline (see manuscript) were
 selected for with CV and geom_mean information attached in
 ../Data/Orthologs/Verified_Candidates.csv ``` ### Details for
 8\_ExtractPromUTR(Orthologs).ipynb Python packages: ``` -
 gffpandas.gffpandas 0.1.0 - pandas 2.0.0 ``` Purpose: * Similar to
 3\_ExtractPromUTR(ALL\_Transcripts).ipynb, parse gff3 file for just the
 Verified\_Candidates list, and extract promoters. inputs: ``` -
 ../Data/Orthologs/Verified_Candidates.csv generated from
 7_Processing_orthologs.Rmd - Gff3 files from
 ../Data/Genomes/[Genus_species]/*.gff3 - reference genome from
 ../Data/Genomes/[Genus_species]/*.fa ``` outputs: ``` - extracted gff3
 info into ../Data/Promoters/[Genus_species]_orthologs_annotation.csv (from
 extract_prom_utr function) - extracted sequences into
 ../Data/Promoters/[Genus_species]_orthologs_sequence.csv (from
 extract_sequence function) ``` ### Details for 9\_Motif\_Scan.Rmd R
 packages: ``` - readr 2.1.4 - tidyverse 2.0.0 - tidyr 1.3.0 - tibble 3.2.1
 - stringr 1.5.0 - universalmotif 1.16.0 ``` Purpose: * First merge
 Verified\_Candidates list with promoter info, then label the candidate
 promoters as TATA, Coreless, Ypatch, Inr similar to
 4\_Label\_Promoters.Rmd * The code is modified from Jores et al., Nat.
 Plants 2021 Inputs: ``` - Verified_candidates list from
 ../Data/Orthologs/Verified_candidates.csv generated by
 7_Processing_orthologs.Rmd - "promoter" and UTR sequence from
 ../Data/Promoters/[Genus_species]_orthologs_sequence.csv generated by
 8_ExtractPromUTR(Orthologs).ipynb Note the "promoter" here is up
 to 2kb of the intergenic region as defined in
 8_ExtractPromUTR(Orthologs).ipynb ``` Outputs: ``` - Verified_candidates
 list merged with sequence info:
 ../Data/Orthologs/Verified_candidates_seq.csv - Motif scan merged with
 Verified_candidates list in
 ../Data/Core_Promoters/Orthologs/Motif_results.csv ``` ### Details for
 10\_Octamer\_Scan.ipynb Python packages: ``` - pandas 2.0.0 - Biopython
 1.81 ``` Purpose: * Scans input promoter sequences for TATA, Y patch, CA,
 GA octamer motifs * The octamers are defined by Yamamoto et al., 2009
 Inputs: ``` - Verified_candidates_seq list from
 ../Data/Orthologs/Verified_candidates_seq.csv generated by
 8_ExtractPromUTR(Orthologs).ipynb Note the "promoter" region in
 the dataframe is up to 2kb of intergenic region as defined in
 8_ExtractPromUTR(Orthologs).ipynb ``` Outputs: ``` - Output result from
 motif scanning: ../Data/Core_Promoters/Orthologs/Octamer_results.csv ```
 ### Details for 11\_Append\_Geom\_mean\_PR.Rmd R packages: * tidyverse
 2.0.0 Purpose: * Finish building the final dataframe: 1. Merge Motif
 screen results with Octamer scan results 2. Append geom\_mean PR info
 Inputs: ``` - Motif scan output from
 ../Data/Core_Promoters/Orthologs/Motif_results.csv - Octomer scan output
 from ../Data/Core_Promoters/Orthologs/Octamer_results.csv ``` Outputs: ```
 - ../Data/Summary.csv is written. It is turned into xlsx and conditionally
 formatted to become Supplemental_Table_S3 ``` ## Sharing/Access
 information Data was derived from the following sources: RNA-seq Datasets
 were retrived from NCBI via their BioProject number: PRJNA268115,
 PRJNA314076, PRJNA194429, PRJNA231618, PRJDB6414, PRJNA291488,
 PRJNA413872, PRJNA79597, PRJNA221782, PRJNA277076, PRJNA389300, SRA558272,
 SRA558514, SRA558539, PRJNA171684, SRP010680, PRJNA507622, PRJNA691387,
 PRJEB25639, PRJEB2430 Reference transcriptomes and reference genomes were
 downloaded from the Ensembl Plants (http://plants.ensembl.org/index.html)
 for Arabidopsis thaliana, Camelina sativa, Cucumis melo, Glycine max,
 Phaseolus vulgaris, Pisum sativum, Vigna unguiculata, Sorghum bicolor, Zea
 mays, Solanum lycopersicum, Actinidia chinensis, Triticum aestivum. and
 Phytozome (https://phytozome-next.jgi.doe.gov) for Arachis hypogaea, Cicer
 arietinum, and Solanum tuberosum.