DOI: 10.5061/dryad.83bk3jb0q
title: Data from: Nascent transcription reveals regulatory changes in extremophile fishes inhabiting hydrogen sulfide-rich environments

description: Regulating transcription allows organisms to respond to their environment,
 both within a single generation (plasticity) and across generations
 (adaptation). We examined transcriptional differences in gill tissues of
 fishes in the Poecilia mexicana species complex (family Poeciliidae),
 which have colonized toxic springs rich in hydrogen sulfide (H2S) in
 southern Mexico. There are gene expression differences between sulfidic
 and non-sulfidic populations, yet regulatory mechanisms mediating this
 gene expression variation remain poorly studied. We combined capped-small
 RNA sequencing (csRNA-seq), which captures actively transcribed (i.e.,
 nascent) transcripts, and traditional messenger RNA sequencing (mRNA-seq)
 to examine how variation in transcription, enhancer activity, and
 associated transcription factor binding sites may facilitate adaptation to
 extreme environments. csRNA-seq revealed thousands of differentially
 initiated transcripts between sulfidic and non-sulfidic populations, many
 of which are involved in H2S detoxification and response. Analyses of
 transcription factor binding sites in promoter and putative enhancer
 csRNA-seq peaks identify a suite of transcription factors likely involved
 in regulating H2S-specific shifts in gene expression, including several
 key transcription factors known to respond to hypoxia. Our findings
 uncover a complex interplay of regulatory processes that reflect the
 divergence of extremophile populations of P. mexicana from their
 non-sulfidic ancestors and suggest shared responses among evolutionarily
 independent lineages.
description: <em>Sample collection</em> Adult
 female fish were collected via seine net from sulfidic and non-sulfidic
 sites in the Pichucalco, Puyacatengo, and Tacotalpa drainages of the Río
 Grijalva basin in southern Mexico (<em>N</em> = 5–6 per site;
 electronic supplementary material, table S1) as part of a previous study
 by Kelley <em>et al.</em> (Kelley et al., 2016). In the
 Pichucalco drainage, sulfidic springs are inhabited by the highly endemic
 extremophile <em>Poecilia sulphuraria </em>(sulphur molly;
 figure 1b, top) (Alvarez del Villar, 1948), but non-sulfidic springs are
 populated by <em>P. mexicana </em>(figure 1b, bottom). In the
 Puyacatengo and Tacotalpa drainages, <em>P. mexicana</em>
 inhabits both the sulfidic and non-sulfidic sites (figure 1c–d). Fish were
 sacrificed, and their gill arches were extracted, preserved in RNAlater
 (Invitrogen, Inc), and stored at -80 ºC. Gill tissues were selected
 because they are directly exposed to H<sub>2</sub>S in the
 water and were previously shown to have significant gene expression
 differences between the sulfidic and non-sulfidic ecotypes (Kelley et al.,
 2016; Passow et al., 2017). While inter-drainage
 variation in adaptation to H<sub>2</sub>S has been
 demonstrated in <em>Poecilia</em> (Greenway et al., 2020;
 Pfenninger et al., 2014), we were most interested in the broad patterns of
 gene regulation across ecotypes, as there is pre-established convergence
 in the differential expression of H<sub>2</sub>S-related genes
 in all sulfidic populations (Kelley et al., 2016; Passow et al., 2017).
 For this reason, we compared gill tissues sampled from all sulfidic fish
 to those from all non-sulfidic fish combined across drainages.
 <em>csRNA library preparation and
 sequencing</em> csRNA-seq was generated from a
 subset of the fish samples described above (<em>N </em>= 6
 fish per population; electronic supplementary material, table S1). Total
 RNA was extracted from gill tissues using Qiagen’s miRNeasy Mini Kit
 following all standard protocols. Total RNA concentrations were estimated
 with the Qubit RNA HS Assay Kit and the Agilent 2100 BioAnalyzer using the
 RNA 6000 Nano Kit. Short RNAs (~20–60 nucleotides) were isolated using
 size selection with gel electrophoresis. For each sample, a 10% aliquot of
 the size-selected RNA suspension was reserved for use as an “input
 library” (Duttke, Beyhan, et al., 2022). These aliquot are sequenced and
 used in the HOMER pipeline to control for exonic contamination when
 identifying transcription initiation sites (Duttke et al., 2019). For the
 remaining size-selected RNA for each sample, 5’ 7-methylguanosine capped
 RNA was isolated following protocols from Duttke <em>et
 al</em>. (Duttke et al., 2019), and this constituted our csRNA
 libraries. csRNA and the reserved input aliquots were converted into cDNA
 libraries for each sample, amplified with PCR, size selected to remove
 primer dimers and purified with gel electrophoresis, pooled, and
 single-end 75-bp reads were sequenced using an Illumina NextSeq
 500. <em>Analysis of  differential transcription
 initiation using csRNA-seq</em> We used FastQC
 (v. 0.11.9) to examine the quality of the raw csRNA and input RNA reads
 (Andrews, 2010). Subsequent analyses were carried out using the csRNA-seq
 analysis pipeline in HOMER (v. 4.11) (Duttke et al., 2019). We used
 HOMER’s <em>trim</em> to remove Illumina TruSeq adapters from
 the 3’ end of reads, discard reads less than 20 bp, and trim the ends of
 reads below a Phred score of 20. We indexed the <em>P.
 mexicana</em> reference genome (GenBank assembly accession:
 GCF_001443325.1) (Warren et al., 2018) using STAR (v. 2.7.6a) and then
 aligned the reads with STAR, outputting one primary alignment per read
 (Dobin et al., 2013). Tag directories were created for the uniquely mapped
 csRNA libraries and the input libraries for each sample using HOMER’s
 <em>makeTagDirectory</em>, including the
 <em>-single</em> flag to account for the large number of
 scaffolds in the reference genome and <em>-fragLength </em>set
 to 30 bp to reflect the average length of the nascent RNAs captured by the
 csRNA-seq protocol (Duttke et al., 2019). A third tag directory was
 created for each sample using the mapped mRNA libraries (described below)
 and used to control for false positives and exonic contaminants in the
 next step. The csRNA, input, and mRNA tag directories were input into
 HOMER’s <em>findcsRNATSS.pl</em> to identify sites of
 transcription initiation, hereafter called peaks. This step of the HOMER
 pipeline also determines the stability of initiated transcripts by
 assessing whether stable mRNA transcripts map within -100 bp to +500 bp of
 each csRNA-seq peak. Unstable, nascent transcripts are not processed into
 stable mRNAs and are thus determined by the lack of stable mRNA
 transcripts mapping to the adjacent region (see (Duttke et al., 2019) for
 additional detail). We retained peaks with at least seven csRNA-sequencing
 reads per 10 million reads for further analysis following the default
 parameters. Peaks were annotated using the GTF file for <em>P.
 mexicana</em> [GenBank assembly accession: GCF_001443325.1,
 converted from GFF using <em>gffread </em>(v. 0.9.9)]. HOMER’s
 <em>mergePeaks</em> was used to create a non-redundant list of
 peaks across all samples, condensing peaks that directly overlapped. The
 merged peaks file was used with homerTools
 <em>annotatePeaks.pl</em> and the <em>-strand +
 -fragLength 1</em> options required for csRNA-seq to quantify the
 raw read counts for each sample, associate each peak with the nearest
 annotated gene from the GFF, and identify the region of the genome they
 mapped to relative to this gene [i.e., in an exon, intergenic, intron,
 promoter or transcription start site (promoter-TSS), or transcription
 termination site (TTS)]. The promoter-TSS region is defined as the region
 1000 bp upstream through 100 bp downstream of a gene’s annotated TSS in
 the GTF file. We the identified differentially initiated peaks in R using
 edgeR (v. 4.0.12) (Robinson et al., 2010) by comparing all sulfidic
 samples to all non-sulfidic samples while including drainage and species
 as covariates. Significantly differentially initiated peaks were defined
 as those with a false discovery rate (FDR) &lt; 0.05. Variation in
 transcription initiation across samples was plotted using a principal
 components analysis (PCA), using the 500 peaks with the highest average
 signal across all samples using the R package <em>limma</em>
 (v. 3.46.0). Hierarchically clustered heatmaps of csRNA-seq peaks were
 generated with the R package <em>pheatmap</em> (v.b1.0.12)
 using the top 5,000 peaks based on overall number of counts across samples
 and subsequently with only significantly differentially initiated
 peaks. We used a curated a list of genes with Gene
 Ontology (GO) terms related to H<sub>2</sub>S detoxification
 or response (GO:0070221, GO:0006790, GO:0044273, GO:0070813, GO:0000096,
 GO:0006090, GO:0006749, GO:1901687, GO:0000098, GO:0008272, GO:0006534,
 GO:0044272) to carry out a Fisher’s exact test using
 <em>fisher.test()</em> in base R (v. 4.0.3), which tested
 whether peaks in or near genes related to H<sub>2</sub>S
 detoxification and response were overrepresented in the list of
 significantly differentially initiated peaks.
 <em>mRNA library preparation and
 sequencing</em> Samples were prepared and
 sequenced as part of a previous study by Kelley <em>et
 al.</em> (Kelley et al., 2016). In brief, total RNA was extracted
 from gill tissues (<em>N</em> = 17 sulfidic fish, <em>N
 </em>= 18 non-sulfidic fish; electronic supplementary material,
 table S1), poly-A messenger RNA (mRNA) were isolated, and libraries were
 sequenced using an Illumina HiSeq 2000 using paired-end 101-bp
 reads. <em>mRNA differential gene
 expression</em> The quality of all raw reads was
 assessed using FastQC. Trim Galore! (v. 0.4.2) was used to remove the
 first 11 bp from the 5’ ends of reads (<em>--clip_R1 11 --clip_R2
 11</em>) and perform adapter trimming (<em>--stringency
 6</em>) but not quality trimming (<em>--quality 0</em>)
 (Krueger, 2014). A second round of Trim Galore! was used to quality trim
 the ends of reads with a Phred score below 24 (<em>--quality
 24</em>) and remove reads less than 50 bp long (<em>--length
 50</em>). HISAT2 <em>hisat2-build</em> (v. 2.1.0) was
 used to index the <em>P. mexicana</em> reference genome with
 the mitochondrial genome appended (GenBank accession: KC992995.1)
 (Pfenninger et al., 2014). Trimmed reads were mapped using the
 <em>--downstream-transcriptome-assembly</em> flag for
 compatibility with StringTie (v. 2.0.3), <em>--fr</em>
 specifying the mate pair orientation, and read group information included
 for sample, read group identifier, and platform (Kim et al., 2019).
 SAMtools <em>view</em> (v. 1.9) was used to convert the
 resulting SAM files to BAM files, <em>sort</em> to sort the
 BAM files by coordinate, and <em>merge</em> to combine samples
 with more than one technical replicate (Li et al., 2009). StringTie (v.
 2.0.3) with the <em>-eB</em> flag was used to generate read
 coverage tables in GTF format for each sample (Pertea et al., 2016). An
 edited version of the prepDE.py script from StringTie (see code
 repository) was used to generate gene counts matrices, resulting in counts
 for nuclear genes only. Genes were filtered by a minimum
 counts-per-million (cpm) of at least five in the smallest library and to
 have non-zero counts in at least five samples. Using edgeR (v. 4.0.12),
 libraries were normalized using the trimmed mean of M-values (TMM, the
 default) and quasi-likelihood F-tests were used to compare differential
 gene expression in sulfidic samples compared to non-sulfidic samples with
 species and drainage included as covariates (McCarthy et al., 2012;
 Robinson et al., 2010). <em>Nascent transcription
 near regions with stable gene expression </em> We
 investigated if patterns in differential transcript initiation between
 ecotypes matched patterns in nearby gene expression. Peaks with
 differential transcript initiation from csRNA-seq (FDR &lt; 0.05) were
 intersected with differentially expressed genes from mRNA-seq (FDR
 &lt; 0.05) based on their gene annotation to identify genomic regions
 that showed both changes in nascent transcription and stable gene
 expression between ecotypes. We then examined directionality to see if
 csRNA and mRNA datasets were both upregulated in sulfidic samples compared
 to non-sulfidic samples, both downregulated, or divergent. Those that were
 upregulated in both datasets or downregulated in both datasets are
 hereafter referred to as the “parallel upregulated subset” and “parallel
 downregulated subset”, respectively. <em>TF
 binding site enrichment in differentially initiated promoter
 regions</em> We identified TF binding sites in
 differentially initiated csRNA-seq peaks in the promoter-TSS regions of DE
 genes, focusing specifically on the upregulated and downregulated subsets
 defined above. We used homerTools <em>findMotifsGenome.pl</em>
 to extract FASTA sequences spanning a 200-bp region from 150 bp upstream
 through 50 bp downstream of the peak center, which should approximate the
 site of transcription initiation. We chose this narrow window in order to
 conservatively identify only the primary TF binding sites involved in
 initiating gene expression. We used CiiiDER (Gearing et al., 2019) to
 identify significantly enriched TF binding motifs in these regions using
 the JASPAR CORE 2022 non-redundant vertebrate transcription factors
 database (Sandelin et al., 2004). We used a maximum deficit of 0.15 to
 allow for non-exact motif matching, and used a stringent p-value cutoff of
 0.01 for classifying significantly enriched binding sites. A set of
 promoter-TSS regions that were not differentially initiated (FDR &gt;
 0.05 and log2 fold-change &lt; 0.5) and were not associated with a
 differentially expressed gene between ecotypes were used as the background
 for these enrichment analyses (n = 7,443).
 <em>Identification and analysis of putative enhancer
 regions</em> csRNA-sequencing can capture
 enhancer RNAs, small unstable transcripts that result from active enhancer
 regions (Duttke, Montilla-Perez, et al., 2022; Lim et al., 2021). Enhancer
 RNA abundance often correlates with the expression of an enhancer’s target
 gene(s) (Arnold et al., 2020; Azofeifa et al., 2018); we therefore
 correlated differential initiation of putative enhancer csRNA-seq peaks
 with the expression of nearby genes to identify putative enhancers
 exhibiting differential activity between ecotypes and to infer genes that
 may be targeted by these enhancers. csRNA-seq peaks putatively
 corresponding to enhancer regions were defined as peaks annotated as
 distal from nearby genes (&gt;500 bp from the promoter region of the
 nearest gene) and being comprised of unstable transcripts (i.e.,
 transcripts present in csRNA but not mRNA libraries). This list of peaks
 was then intersected with the results of differential initiation analysis
 to identify putative enhancers with significantly increased or decreased
 activity in the sulfidic ecotype. We then identified
 all significantly differentially expressed genes with a TSS located less
 than 500 kilobase pairs (kbp) from the center of an ecotype-responsive
 putative enhancer peak by first expanding all ecotype-responsive putative
 enhancer peaks by 500 kbp in both directions using the
 <em>slop</em> tool within bedtools (v. 2.30.0) and then
 intersecting these expanded regions with a bedfile of gene TSS positions
 using bedtools <em>intersect</em>. The resulting pairs of
 putative enhancers and potential target genes were filtered to retain
 those with the same direction of differential initiation/expression in the
 sulfidic ecotype (i.e., upregulated csRNA-seq peak and upregulated gene
 expression, and vice versa), and the log2 fold-change values of csRNA-seq
 and mRNA-seq were correlated using linear regressions in R. To focus
 downstream analyses only on putative target genes with high correlations
 to putative enhancer regions, pairs of putative enhancers and target genes
 were filtered to retain only those falling within one standardized
 residual of the regression line.  For the resulting
 putative enhancer peaks that were highly correlated with the expression of
 their inferred target gene(s), fasta sequences were extracted using the
 getfasta tool from bedtools. A background set of 1,000 peaks were randomly
 selected from all putative enhancer peaks that were not differentially
 initiated between ecotypes (i.e., distal csRNA-seq peaks comprised of
 unstable transcripts and with an FDR &gt; 0.05 in differential
 initiation analyses). CiiiDER was again used to identify significantly
 enriched TF binding motifs in these regions using the JASPAR CORE 2022
 non-redundant vertebrate transcription factors database. We used a maximum
 deficit of 0.15 to allow for non-exact motif matching, and used a
 stringent p-value cutoff of 0.01 for defining significantly enriched
 binding sites. References Alvarez
 del Villar, J. (1948). <em>Descripción de una nueva especie de
 Mollienisiacapturada en Baños del Azufre, Tabasco (Pisces,
 Poeciliidae)</em>. Andrews, S. (2010).
 <em>FastQC:  A Quality Control Tool for High Throughput Sequence
 Data [Online]</em>.
 http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
 Arnold, P. R., Wells, A. D., & Li, X. C. (2020).
 Diversity and emerging roles of enhancer RNA in regulation of gene
 expression and cell fate. <em>Frontiers in Cell and Developmental
 Biology</em>, <em>7</em>, 377.
 Azofeifa, J. G., Allen, M. A., Hendrix, J. R., Rubin, J. D.,
 & Dowell, R. D. (2018). Enhancer RNA profiling predicts
 transcription factor activity. <em>Genome Research</em>,
 <em>28</em>(3), 334–344. Dobin, A., Davis,
 C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, P.,
 Chaisson, M., & Gingeras, T. R. (2013). STAR: ultrafast universal
 RNA-seq aligner. <em>Bioinformatics</em>,
 <em>29</em>(1), 15–21. Duttke, S. H.,
 Beyhan, S., Singh, R., Neal, S., Viriyakosol, S., Fierer, J., Kirkland, T.
 N., Stajich, J. E., Benner, C., & Carlin, A. F. (2022). Decoding
 transcription regulatory mechanisms associated with coccidioides immitis
 phase transition using total RNA. <em>Msystems</em>,
 <em>7</em>(1), e01404-21. Duttke, S. H.,
 Chang, M. W., Heinz, S., & Benner, C. (2019). Identification and
 dynamic quantification of regulatory elements using total RNA.
 <em>Genome Research</em>, <em>29</em>(11),
 1836–1846. Duttke, S. H., Montilla-Perez, P., Chang, M.
 W., Li, H., Chen, H., Carrette, L. L. G., Guglielmo, G. de, George, O.,
 Palmer, A. A., & Benner, C. (2022). Glucocorticoid
 receptor-regulated enhancers play a central role in the gene regulatory
 networks underlying drug addiction. <em>Frontiers in
 Neuroscience</em>, <em>16</em>, 858427.
 Gearing, L. J., Cumming, H. E., Chapman, R., Finkel, A. M.,
 Woodhouse, I. B., Luu, K., Gould, J. A., Forster, S. C., &
 Hertzog, P. J. (2019). CiiiDER: A tool for predicting and analysing
 transcription factor binding sites. <em>PloS One</em>,
 <em>14</em>(9), e0215495. Greenway, R.,
 Barts, N., Henpita, C., Brown, A. P., Arias Rodriguez, L., Rodriguez Pena,
 C. M., Arndt, S., Lau, G. Y., Murphy, M. P., Wu, L., Lin, D., Tobler, M.,
 Kelley, J. L., & Shaw, J. H. (2020). Convergent evolution of
 conserved mitochondrial pathways underlies repeated adaptation to extreme
 environments. <em>Proc Natl Acad Sci U S A</em>.
 https://doi.org/10.1073/pnas.2004223117 Kelley, J. L.,
 Arias-Rodriguez, L., Patacsil Martin, D., Yee, M. C., Bustamante, C. D.,
 & Tobler, M. (2016). Mechanisms Underlying Adaptation to Life in
 Hydrogen Sulfide-Rich Environments. <em>Mol Biol Evol</em>,
 <em>33</em>(6), 1419–1434.
 https://doi.org/10.1093/molbev/msw020 Kim, D., Paggi,
 J. M., Park, C., Bennett, C., & Salzberg, S. L. (2019).
 Graph-based genome alignment and genotyping with HISAT2 and
 HISAT-genotype. <em>Nat Biotechnol</em>,
 <em>37</em>(8), 907–915.
 https://doi.org/10.1038/s41587-019-0201-4 Krueger, F.
 (2014). <em>Trim Galore!: A wrapper tool around Cutadapt and FastQC
 to consistently apply quality and adapter trimming to FastQ
 files</em> (0.3.7). Li, H., Handsaker, B.,
 Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G.,
 & Durbin, R. (2009). The sequence alignment/map format and
 SAMtools. <em>Bioinformatics</em>,
 <em>25</em>(16), 2078–2079. Lim, J.-Y.,
 Duttke, S. H., Baker, T. S., Lee, J., Gambino, K. J., Venturini, N. J.,
 Ho, J. S. Y., Zheng, S., Fstkchyan, Y. S., & Pillai, V. (2021).
 DNMT3A haploinsufficiency causes dichotomous DNA methylation defects at
 enhancers in mature human immune cells. <em>Journal of Experimental
 Medicine</em>, <em>218</em>(7), e20202733.
 McCarthy, D. J., Chen, Y., & Smyth, G. K. (2012).
 Differential expression analysis of multifactor RNA-Seq experiments with
 respect to biological variation. <em>Nucleic Acids Res</em>,
 <em>40</em>(10), 4288–4297.
 https://doi.org/10.1093/nar/gks042 Passow, C. N.,
 Henpita, C., Shaw, J. H., Quackenbush, C. R., Warren, W. C., Schartl, M.,
 Arias‐Rodriguez, L., Kelley, J. L., & Tobler, M. (2017). The roles
 of plasticity and evolutionary change in shaping gene expression variation
 in natural populations of extremophile fish. <em>Molecular
 Ecology</em>, <em>26</em>(22), 6384–6399.
 Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., &
 Salzberg, S. L. (2016). Transcript-level expression analysis of RNA-seq
 experiments with HISAT, StringTie and Ballgown. <em>Nature
 Protocols</em>, <em>11</em>(9), 1650–1667.
 Pfenninger, M., Lerp, H., Tobler, M., Passow, C., Kelley, J. L.,
 Funke, E., Greshake, B., Erkoc, U. K., Berberich, T., & Plath, M.
 (2014). Parallel evolution of cox genes in H2S-tolerant fish as key
 adaptation to a toxic environment. <em>Nat Commun</em>,
 <em>5</em>, 3873. https://doi.org/10.1038/ncomms4873
 Robinson, M. D., McCarthy, D. J., & Smyth, G. K. (2010).
 edgeR: a Bioconductor package for differential expression analysis of
 digital gene expression data. <em>Bioinformatics</em>,
 <em>26</em>(1), 139–140. Sandelin, A.,
 Alkema, W., Engström, P., Wasserman, W. W., & Lenhard, B. (2004).
 JASPAR: an open‐access database for eukaryotic transcription factor
 binding profiles. <em>Nucleic Acids Research</em>,
 <em>32</em>(suppl_1), D91–D94. Warren, W.
 C., Garcia-Perez, R., Xu, S., Lampert, K. P., Chalopin, D., Stock, M.,
 Loewe, L., Lu, Y., Kuderna, L., Minx, P., Montague, M. J., Tomlinson, C.,
 Hillier, L. W., Murphy, D. N., Wang, J., Wang, Z., Garcia, C. M., Thomas,
 G. C. W., Volff, J. N., … Schartl, M. (2018). Clonal polymorphism and high
 heterozygosity in the celibate genome of the Amazon molly. <em>Nat
 Ecol Evol</em>, <em>2</em>(4), 669–679.
 https://doi.org/10.1038/s41559-018-0473-y
description: # Data from: Nascent transcription reveals regulatory changes in
 extremophile fishes inhabiting hydrogen sulfide-rich environments This
 repo contains code for all data processing and analysis used in our study
 of transcription and gene expression in H2S-adapted populations of *P.
 mexicana*. ## Description of Data and File Structure ###
 Table\_S5\_01.30.24.xlsx * **csRNA Peak ID:** The identifier assigned to
 the csRNA-seq peak. * **Chromosome:** The chromosomal scaffold in the *P.
 mexicana* reference genome the peak mapped to. * **Start:** The start
 coordinate of the peak. * **End:** The end coordinate of the peak. *
 **Strand:** Strand (+/-) the read mapped to. * **Annotation:** The type of
 genomic element the read mapped to [e.g., intergenic, intron, exon,
 promoter-transcription start site (TSS), transcription termination site
 (TTS)]. * **Nearest Gene:** The name of the closest annotated gene to the
 peak. * **Distance to TSS (bp):** The distance in base pairs (bp) to the
 nearest annotated TSS. * **Log2 Fold Change:** The log2-fold change in
 transcript initiation levels comparing sulfidic fish to non-sulfidic fish.
 * ***p-value***: Significance. * **FDR:** The false discovery rate. *
 **Human Homologue:** The name of the closest human homologues to the *P.
 mexicana* genes from a BLAST search. * **Human Protein Annotation:** A
 detailed description of the human homologue. ### Table\_S9\_01.30.24.xlsx
 * **Gene Name:** The name of the closest annotated gene to the peak. *
 **Human Homologue:** The name of the closest human homologues to the *P.
 mexicana* genes from a BLAST search. * **Human Protein Annotation:** A
 detailed description of the human homologue. * **Log2 Fold Change:** The
 log2-fold change in transcript initiation levels comparing sulfidic fish
 to non-sulfidic fish. * ***p-value***: Significance. * **FDR**: The false
 discovery rate. ### Table\_S10\_02.02.24.xlsx * **Peak Group:**
 “Upregulated Promoter-TSS Peaks” indicate TFBS enriched in H2S-upregulated
 csRNA-seq peaks in the promoter-TSS of H2S-upregulated genes.
 “Downregulated Promoter-TSS Peaks” indicate TFBS enriched in
 H2S-downregulated csRNA-seq peaks in the promoter-TSS of H2S-downregulated
 genes. * **TF Motif ID:** Binding site motif ID from the JASPAR database
 (Mathelier, Fornes et al. 2016). * **TF Name:** TF that binds to the
 motif. The :: separates co-binding TFs. * **TF Family:** Family of
 transcription factors to which a TF belongs. * **TF Class:** Class of
 transcription factors to which a TF belongs. * **Deficit:** The difference
 between the MATCH score (Kel et al. 2003) of a TF binding site motif and a
 maximum score of 1. TFs are displayed at the deficit with the lowest
 *p*-value with a maximum deficit of 0.15. A deficit of zero equals an
 exact match. * **Number of Foreground Peaks with Motif:** Number of
 differentially initiated promoter-TSS peaks in a given Peak Group (see
 column 1) in which a given TF binding site motif was identified. *
 **Number of Background Peaks with Motif:** Number of peaks in the
 background set of non-differentially initiated promoter-TSS peaks in which
 a given TF binding site motif was identified.   * **Proportion of Target
 Sites:** The proportion of differentially initiated promoter-TSS peaks in
 a given Peak Group (see column 1) in which a given TF binding site motif
 was identified. * ***p-value***: *p*-value resulting from a Fisher’s exact
 test of TF binding site enrichment in the foreground peaks compared to the
 background peaks. * **Log2 Enrichment:** The larger the value (&gt; 0),
 the more overrepresented the motif in the foreground compared to the
 background. The smaller the value (&lt; 0) the more underrepresented. *
 **Significance Score**: -log10(*p*-value) such that higher Significance
 Score values indicate lower *p*-values. Used for plotting. ###
 Table\_S11\_02.02.24.xlsx * **TF Motif ID:** Binding site motif ID from
 the JASPAR database (Mathelier, Fornes et al. 2016). * **TF Name:** TF
 that binds to the motif. The :: separates co-binding TFs. * **Log2
 Enrichment:** The larger the value (&gt; 0), the more overrepresented the
 motif in the foreground compared to the background. The smaller the value
 (&lt; 0) the more underrepresented. * **Enrichment *p*-value:** *p*-value
 resulting from a Fisher’s exact test of TF binding site enrichment in the
 foreground peaks compared to the background peaks. * **Peak Group:**
 “Upregulated Promoter-TSS Peaks” indicate TFBS enriched in H2S-upregulated
 csRNA-seq peaks in the promoter-TSS of H2S-upregulated genes.
 “Downregulated Promoter-TSS Peaks” indicate TFBS enriched in
 H2S-downregulated csRNA-seq peaks in the promoter-TSS of H2S-downregulated
 genes. * **TF Encoding Gene:** Gene in the *P. mexicana* genome encoding a
 given TF. * **TF Encoding Gene Log2 Fold Change:** Log2 transformed
 fold-change value of a given TF-encoding gene from differential gene
 expression analysis. Positive values indicate upregulation in the sulfidic
 ecotype. * **TF Encoding Gene *p*-value:** Uncorrected p-value from
 differential gene expression analysis. * **TF Encoding Gene FDR:** False
 discovery rate from differential gene expression analysis. * **Human
 Homologue**: The name of the closest human homologues to the *P. mexicana*
 genes from a BLAST search. ### Table\_S12\_02.05.24.xlsx * **csRNA Peak
 ID:** The identifier assigned to the csRNA-seq peak. * **Chromosome:** The
 chromosomal scaffold in the *P. mexicana* reference genome the peak mapped
 to. * **csRNAseq Peak Start:** The start coordinate of the peak. *
 **csRNAseq Peak End:** The end coordinate of the peak. * **csRNAseq Peak
 Strand:** Strand (+/-) the read mapped to. * **csRNAseq Log2 Fold
 Change:** The log2-fold change in transcript initiation levels comparing
 sulfidic fish to non-sulfidic fish, with positive values indicating higher
 expression in the sulfidic ecotype. * **csRNAseq Peak FDR:** The false
 discovery rate from analysis of differential initiation in the csRNAseq
 peak * **Inferred Target Gene**: Putative target genes defined by being
 within ± 500 kbp from the putative enhancer csRNA-seq peak, significantly
 differentially expressed, and with a high positive correlation between the
 target gene and putative enhancer log2 fold change values (see methods). *
 **Human Homologue:** The name of the closest human homologues to the *P.
 mexicana* genes from a BLAST search. * **Human Protein Annotation:** A
 detailed description of the human homologue. * **Target Gene Log2 Fold
 Change:** The log2-fold change in gene expression comparing sulfidic fish
 to non-sulfidic fish, with positive values indicating higher expression in
 the sulfidic ecotype. * **Target Gene FDR:** The false discovery rate from
 analysis of differential gene expression. ### Table\_S13\_02.05.24.xlsx *
 **Peak Group:** “Upregulated Enhancer-Gene Pairs” indicate TFBS enriched
 in H2S-upregulated putative enhancer csRNA-seq peaks inferred to target
 H2S-upregulated genes. “Downregulated Enhancer-Gene Pairs” indicate TFBS
 enriched in H2S-downregulated putative enhancer csRNA-seq inferred to
 target H2S-downregulated genes. * **TF Motif ID:** Binding site motif ID
 from the JASPAR database (Mathelier, Fornes et al. 2016). * **TF Name:**
 TF that binds to the motif. The :: separates co-binding TFs. * **TF
 Family:** Family of transcription factors to which a TF belongs. * **TF
 Class:** Class of transcription factors to which a TF belongs. *
 **Deficit:** The difference between the MATCH score (Kel et al. 2003) of a
 TF binding site motif and a maximum score of 1. TFs are displayed at the
 deficit with the lowest *p*-value with a maximum deficit of 0.15. A
 deficit of zero equals an exact match. * **Number of Foreground Peaks with
 Motif:** Number of differentially initiated putative enhancers in a given
 Peak Group (see column 1) in which a given TF binding site motif was
 identified. * **Number of Background Peaks with Motif:** Number of peaks
 in the background set of non-differentially initiated putative enhancer
 peaks in which a given TF binding site motif was identified.  *
 **Proportion of Target Sites:** The proportion of differentially initiated
 putative enhancer peaks in a given Peak Group (see column 1) in which a
 given TF binding site motif was identified. * ***p-value***: *p*-value
 resulting from a Fisher’s exact test of TF binding site enrichment in the
 foreground peaks compared to the background peaks. * **Log2 Enrichment**:
 The larger the value (&gt; 0), the more overrepresented the motif in the
 foreground compared to the background. The smaller the value (&lt; 0) the
 more underrepresented. * **Significance Score**: -log10(*p*-value) such
 that higher Significance Score values indicate lower *p*-values. Used for
 plotting. ### Table\_S14\_02.05.24.xlsx * **TF Motif ID:** Binding site
 motif ID from the JASPAR database (Mathelier, Fornes et al. 2016). * **TF
 Name:** TF that binds to the motif. The :: separates co-binding TFs. *
 **Log2 Enrichment:** The larger the value (&gt; 0), the more
 overrepresented the motif in the foreground compared to the background.
 The smaller the value (&lt; 0) the more underrepresented. * **Enrichment
 *p*-value:** *p*-value resulting from a Fisher’s exact test of TF binding
 site enrichment in the foreground peaks compared to the background peaks.
 * **Peak Group:** “Upregulated Enhancer-Gene Pairs” indicate TFBS enriched
 in H2S-upregulated putative enhancer csRNA-seq peaks inferred to target
 H2S-upregulated genes. “Downregulated Enhancer-Gene Pairs” indicate TFBS
 enriched in H2S-downregulated putative enhancer csRNA-seq inferred to
 target H2S-downregulated genes. * **TF Encoding Gene:** Gene in the *P.
 mexicana* genome encoding a given TF. * **TF Encoding Gene Log2 Fold
 Change:** Log2 transformed fold-change value of a given TF-encoding gene
 from differential gene expression analysis. Positive values indicate
 upregulation in the sulfidic ecotype. * **TF Encoding Gene *p*-value:**
 Uncorrected p-value from differential gene expression analysis. * **TF
 Encoding Gene FDR:** False discovery rate from differential gene
 expression analysis. * **Human Homologue:** The name of the closest human
 homologues to the *P. mexicana* genes from a BLAST search. ### mRNA.zip *
 This zip file contains the raw mRNA-seq read counts
 (gene_count_matrix-2019-10-01_no_STRG.csv) located within the raw_counts
 directory. This file is generated by the code in section **2. mRNA-seq
 trimming and mapping** in the detailed readme below. Each row of this
 table corresponds to a gene (with the first column containing gene
 identifiers), and each column corresponds to a sample. Values in the table
 represent raw, untransformed mRNA-seq read counts for a given gene in a
 given sample.  csRNA.zip * This zip file contains relevant output files
 generated from initial processing of csRNA-seq data via HOMER as described
 in section **1. cRNA-seq QC, mapping, peak calling and annotation** of the
 detailed readme below. These files are organized in the following three
 directories: * **06_homertools_findcsrnatss**: this directory contains
 output files of the readme subsection **1. cRNA-seq QC, mapping, peak
 calling and annotation: Identify transcription initiation peaks** below. *
 This includes six files per sample as follows (see the HOMER documentation
 [here](http://homer.ucsd.edu/homer/ngs/csRNAseq/) for additional detail):
 * *[sampleID].tss.txt* : primary output file, contains positions of TSS
 clusters and annotation information * *[sampleID].stats.txt* : contains
 general information about the run, the command used, parameters selected
 for the analysis. * *[sampleID].tss.txt.freq.tsv* : Nucleotide frequencies
 relative to the primary TSS of each valid TSS cluster. *
 *[sampleID]*.alltss.txt: All putative TSS clusters, not just the valid
 ones, in case you want to go back and look through the others. *
 *[sampleID].inputDistribution.txt* : if you used an input/control sample,
 this will show the CDFs for the threshold optimization. *
 *[sampleID].rnaDistribution.txt* : same as the input, but for an RNA-seq
 control. * Additionally, the **bed** directory contains a bedfile of all
 csRNA-seq peaks identified via HOMER for each sample (i.e., the contents
 of the [sampleID].tss.txt file for each sample transformed into standard
 bedfile format). * **07_homer_mergepeaks**: this directory contains output
 files of the readme subsection **1. cRNA-seq QC, mapping, peak calling and
 annotation:** **Merge peaks** below: * merged.tss.txt: a set of
 non-redundant, merged peaks generated by comparing all csRNA-seq peaks
 identified across all samples.  * merged.tss.bed: same as merged.tss.txt,
 but modified to adhere to a standard bedfile format. *
 allintersects.pmex.txt: file indicating which individual peaks from each
 sample were merged together in the merged.tss.txt/bed files above.  *
 08_homer_annotatepeaks: this directory contains the raw csRNA-seq counts
 for each merged peak region generated in subsection **1. cRNA-seq QC,
 mapping, peak calling and annotation:** **Annotate peaks** below: *
 rawcounts.txt: table containing peak annotation information (see HOMER
 documentation [here](http://homer.ucsd.edu/homer/ngs/annotation.html) for
 additional information) in Columns 1-19, and raw csRNA-seq count values
 for each individual in columns 20-31.  ### Code files * All code files
 ending in .R should be opened with R or Rstudio and are described in
 detail in the Code and Analysis README below. --- ## Sharing/Access
 Information * Newly-generated csRNA-seq data * Raw data available at NCBI
 BioProject: PRJNA743555 * Previously published mRNA-seq data ([Kelley et
 al.
 2016](http://mbe.oxfordjournals.org/content/early/2016/03/01/molbev.msw020)) * Raw data available at NCBI BioProject: [PRJNA290391](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA290391/) --- ## Code and Analysis README - Table of Contents 1. [cRNA-seq trimming, mapping, peak calling, and peak annotation](#1-crna-seq-qc-mapping-peak-calling-and-annotation) 2. [mRNA-seq trimming and mapping ](#2-mrna-seq-trimming-and-mapping) 3. [Analyses of differential transcript initiation (csRNA)](#3-analyses-of-differential-transcript-initiation-csrna) 4. [Analyses of differential gene expression (mRNA)](#4-analyses-of-differential-gene-expression-mrna) 5. [Joint analyses of csRNA and mRNA data](#5-joint-analyses-of-csrna-and-mrna-data) * a. [Characterization of transcript initiation annotations](#a-characterization-of-transcript-initiation-annotations) * b. [Identification of putative enhancer regions](#b-identification-and-analysis-of-putative-enhancer-regions) 6. [TF binding site enrichment analyses](#6-tf-binding-site-enrichment-analyses) * [a. Promoter-TSS csRNA peaks](#a-promoter-tss-csrna-peaks) * [b. Putative enhancers](#b-putative-enhancers) --- ## 1. cRNA-seq QC, mapping, peak calling and annotation These analyses were run on WSU's HPC ([Kamiak](https://hpc.wsu.edu/)). For generalizability, simplified commands are presented here rather than the specific SLURM scripts used to run these commands on Kamiak. Note: The variable `$sample` should be replaced with the full sample name. ### Trim with homerTools *trim* ```bash homerTools trim \ -3 AGATCGGAAGAGCACACGTCT \ -mis 2 \ -minMatchLength 4 \ -min 20 \ -q 20 \ $sample.fastq.gz ``` ### Map csRNA-seq with STAR Reference genome and annotations in GFF format were downloaded from NCBI (Accession: [GCF_001443325.1](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_001443325.1/)). ```bash # Index genome for use with STAR STAR \ --runThreadN 7 \ --runMode genomeGenerate \ --genomeDir genomeDir \ --genomeFastaFiles GCF_001443325.1_P_mexicana-1.0_genomic.fna \ --sjdbGTFfile GCF_001443325.1_P_mexicana-1.0_genomic.gff \ --sjdbOverhang 83 \ --genomeSAindexNbases 13 # Map Reads STAR \ --genomeDir genomeDir \ --runThreadN 24 \ --readFilesIn $sample.fastq.gz.trimmed \ --outFileNamePrefix $sample. \ --outSAMstrandField intronMotif \ --outMultimapperOrder Random \ --outSAMmultNmax 1 \ --outFilterMultimapNmax 10000 \ --limitOutSAMoneReadBytes 10000000 ``` ### Make tag directories Create tag directories for csRNA-seq, input, and mRNA-seq data (see below for details on mRNA-seq processing) ```bash makeTagDirectory \ $sample\_tagDir \ $sample.Aligned.out.sam \ -genome $reference_genome \ -checkGC \ -unique \ -fragLength 30 \ -single ``` ### Identify transcription initiation peaks Peaks are clusters of transcription initiation, called using the "full approach" of the HOMER pipeline. This incorporates mRNA tag directories into the csRNA peak calling to help eliminate false positives like exomic contaminants. ```bash # Convert GFF to GTF format gffread \ -T GCF_001443325.1_P_mexicana-1.0_genomic.gff \ -o GCF_001443325.1_P_mexicana-1.0_genomic.gtf findcsRNATSS.pl \ $sample\_csRNA_tagDir \ -o $sample \ -i $sample\_input_tagDir \ -rna $sample\_rna_tagDir \ -gtf $annotations \ -genome GCF_001443325.1_P_mexicana-1.0_genomic.fna \ -ntagThreshold 7 ``` ### Merge peaks Produce a non-redundant list of peaks from all samples; 12 samples in this study. ```bash mergePeaks \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ $sample.tss.txt \ -strand \ &gt; merged.tss.txt ``` ### Annotate peaks Quantify csRNAseq read counts and annotations for merged peaks. ```bash annotatePeaks.pl \ merged.tss.txt \ GCF_001443325.1_P_mexicana-1.0_genomic.fna \ -gtf GCF_001443325.1_P_mexicana-1.0_genomic.gtf \ -strand + \ -fragLength 1 \ -raw \ -d $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ $sample\_csRNA_tagDir \ &gt; rawcounts.txt ``` --- ## 2. mRNA-seq trimming and mapping The variable $sample should be replaced with the full sample name. ### Trim with TrimGalore ```bash # Trim adapters trim_galore \ --output_dir /path/to/outdir \ --quality 0 \ --fastqc \ --fastqc_args "--nogroup --noextract --outdir /path/to/outdir" \ -path_to_cutadapt /path/to/cutadapt \ --illumina \ --stringency 6 \ --clip_R1 11 \ --clip_R2 11 \ --paired \ --gzip ${sample}_1.fq.gz ${sample}_2.fq.gz # Quality trim trim_galore \ --output_dir /path/to/outdir \ --quality 24 \ --fastqc \ --fastqc_args "--nogroup --noextract --outdir /path/to/outdir" \ -path_to_cutadapt /path/to/cutadapt \ --illumina \ --stringency 6 \ --length 50 \ --paired \ --gzip ${sample}_1_val_1.fq.gz ${sample}_2_val_2.fq.gz ``` ### Map RNA-seq with HISAT2 ```bash # Index the genome hisat2-build \ -f GCF_001443325.1_P_mexicana-1.0_genomic_with_mito_sequence.fa \ GCF_001443325.1_P_mexicana-1.0_genomic_with_mito_sequence # Map hisat2 \ -q \ --threads 4 \ -k 5 \ --downstream-transcriptome-assembly \ --fr \ --rg SM:${sample} \ --rg ID:${rgid} \ --rg PL:ILLUMINA \ -x GCF_001443325.1_P_mexicana-1.0_genomic_with_mito_sequence \ -1 ${sample}_1_val_1_val_1.fq.gz \ -2 ${sample}_2_val_2_val_2.fq.gz \ -S ${sample}_HISat_out.sam \ --summary-file ${sample}_statistics.txt ``` ### Convert, sort, and merge SAM/BAM files ```bash # Convert to BAM samtools view \ -bSh ${sample}_HISat_out.sam &gt; \ ${sample}_HISat_out.bam # Sort by coordinate samtools sort \ -@ 8 \ ${sample}_HISat_out.bam \ -o ${sample}_sorted.bam # Merge BAM files from same individual samtools merge \ -r \ ${sample}_merged_sorted.bam \ ${sample}_sorted.bam \ ${sample}_sorted.bam ``` ### Build gene count matrices Note: modified prepDE.py script ([prepDE-2019-10-01_edited.py](https://github.com/blairperry/pmex_csrnaseq/blob/main/prepDE-2019-10-01_edited.py)) is provided in this repository. ```bash # Generate GTF files for each sample using StringTie and modified version of the associated prepDE.py script mkdir -p ballgown stringtie \ ${sample}_sorted.bam \ -o ballgown/${sample}/${sample}.gtf \ -p 8 \ -G GCF_001443325.1_P_mexicana-1.0_genomic_with_mito_sequence-edited.gff \ -B \ -e # Generate gene count matrix grep -v "STRG" gene_count_matrix.csv &gt; gene_count_matrix_no_STRG.csv ``` --- ## 3. Analyses of differential transcript initiation (csRNA) The following R script contains code used to: * Filter csRNA-seq count data * Perform analysis of differential transcript initiation between sulfidic and non-sulfidic populations using edgeR * Plot heatmaps of csRNA-seq peaks for Figure 2. * Assess differential initiation in candidate H2S detox genes Link to Rscript: 3_csRNA_DiffInitiation_01.15.24.R --- ## 4. Analyses of differential gene expression (mRNA) The following R script contains code used to: * Filter mRNA-seq count data * Perform analysis of differential gene expression between sulfidic and non-sulfidic populations using edgeR Link to Rscript: 4_mRNA_DiffExpression_01.15.24.R --- ## 5. Joint analyses of csRNA and mRNA data ### a. Characterization of transcript initiation annotations The following R script contains code used to: * Summarize annotations of differentially initiated peaks (i.e., whether peaks are located in promoters, exons, etc.) * Identify differentially initiated peaks (csRNA) associated with differentially expressed genes (mRNA) * Test whether specific subsets of differentially initiated peaks (e.g., DI peaks associated with DE genes, DI peaks associated with non-DE genes) exhibit unique distributions of peak annotations * Generate barplots and alluvial plots for Figure 3. Link to Rscript: 5a_peakAnnotationExploration_01.15.24.R ### b. Identification and analysis of putative enhancer regions The following R script contains code used to: * Parse and summarize csRNA peak annotations * Identify putative enhancer regions defined as csRNAseq peaks that are distal from protein coding genes and comprised of unstable transcripts Link to Rscript: 5b_1_putEnhIdentification_01.16.24.R Using the bedfile of all putative enhancer regions identified above, the following commands were used to expand each region by 500 kbp in each direction and intersect with a bedfile of gene TSS positions. Genes found to be located within 500 kbp of a putative enhancer are considered potential target genes that are evaluated further below. ```bash # Sort bedfile of putative enhancer regions bedtools sort -i putativeEnhancerPeaks_pmex_01.12.24.bed &gt; putativeEnhancerPeaks_pmex_01.12.24.sort.bed # Expand each putative enhancer region by 500 kbp in both directions bedtools slop -b 500000 -i putativeEnhancerPeaks_pmex_01.12.24.sort.bed -g ../../data/reference/chrNameLength.txt &gt; putativeEnhancerPeaks_pmex_01.12.24.500kbSlop.bed # Intersect with a bed file of gene TSS positions bedtools intersect -a putativeEnhancerPeaks_pmex_01.12.24.500kbSlop.bed -b ../../data/reference/tss.bed -wa -wb &gt; putativeEnhancerPeaks_pmex_01.12.24.500kbSlop.tssIntersect.txt ``` The following R script contains code used to: * Identify putative enhancer regions with evidence of differential initiation (i.e., differential activity) in csRNA-seq edgeR analyses described above. * Identify putative enhancers and their putative target genes with positively correlated log2 fold-changes (i.e., shared upregulation or downregulation in the sulfidic ecotype compared to nonsulfidic) * Run linear regression and characterize enhancer-gene pairs that fall within 1 standardized residual of the resulting regression line (i.e., highly correlated enhancer-gene pairs) * Generate enhancer-gene correlation dot plot for Figure 5. * Identify candidate H2S detox genes putatively targeted by a putative enhancer. Link to Rscript: 5b_2_putEnhAnalysis_02.02.24.R --- ## 6. TF binding site enrichment analyses ### a. Promoter-TSS csRNA peaks HOMER's findMotifsGenome tool was used to to get regions 50 bp upstream through 150 bp downstream of all peaks, removing any with high repeat content. ```bash findMotifsGenome.pl merged.tss.txt GCF_001443325.1_P_mexicana-1.0_genomic.fna ./ -size -150,50 -dumpFasta # Result fasta is target.fa ``` The following R script contains code used to: * Generate fastas of promoter-TSS csRNAseq peaks that are either: * upregulated in H2S and located in the promoter of an upregulated gene * downregulated in H2S and located in the promoter of a downregulated gene * Generate a fasta of promoter-TSS csRNAseq peaks that are NOT differentially initiated and are located in the promoter of a gene that is not differentially expressed - this is used as the background in TFBS enrichment analyses * Read in and parse result tables from CiiiDER TFBS enrichment analysis * Generate dotplots of top enriched TF binding sites for Figure 4. * Identify TF genes found to be differentially expressed between ecotypes and enriched for binding sites in promoter-TSS peaks. Link to Rscript: 6a_promTSS_TFBS_01.17.24.R ### b. Putative enhancers The following R script contains code used to: * Generate bedfiles for putative enhancers that are either: * upregulated in H2S and inferred to target an upregulated gene * downregulated in H2S and inferred to target a downregulated gene * Generate a bedfile of 1000 putative enhancers that are not differentially initiated and do not target a differentially expressed gene. * Export the bedfiles so that bedtools can be used to pull out fasta sequences (*see commands below*) * Read in and parse result tables from CiiiDER TFBS enrichment analysis * Generate dotplots of top enriched TF binding sites for Figure 5. * Identify TF genes found to be differentially expressed between ecotypes and enriched for binding sites in putative enhancers. Link to Rscript: 6b_enhancer_TFBS_02.04.24.R The following command was used to extract fasta sequences for upregulated and downregulated putative enhancers for use in CiiiDER TFBS enrichment analyses: ```bash mkdir fasta bedtools getfasta -fi ../../data/reference/GCF_001443325.1_P_mexicana-1.0_genomic.fna -bed corr_enhancers_upreg_01.19.24.bed -name &gt; fasta/corr_enhancers_upreg_01.19.24.fa bedtools getfasta -fi ../../data/reference/GCF_001443325.1_P_mexicana-1.0_genomic.fna -bed corr_enhancers_dwreg_01.19.24.bed -name &gt; fasta/corr_enhancers_dwreg_01.19.24.fa bedtools getfasta -fi ../../data/reference/GCF_001443325.1_P_mexicana-1.0_genomic.fna -bed background_enhancers_01.15.24.bed -name &gt; fasta/background_enhancers_01.19.24.fa ```