Estimating cis and trans contributions to differences in gene regulation

In 1961, Jacob and Monod developed a theory of gene reg- ulation in which they distinguished local effects (cis) from distal regulation (trans) (1). Their work immediately raised the question of the relative contributions of these two regu- lation modalities (2, 3). One approach to assessing whether cis or trans regulation is responsible for differences in gene expression between strains or species is to compare differ- ences in expression of genes in parents to allele-specific dif- ferences in F1 hybrids. This approach was explored in (4), who used crosses of C57BL/6J and CAST/Ei mice to study regulatory mechanisms that could explain differences in gene expression between parental strains. Their approach was de- veloped in (5, 6), who used pyrosequencing to study the reg- ulation differences between D. melanogaster and D. simu- lans. With the advent of RNA-seq, genome-wide scans were possible, and (7) examined RNA-seq from F1 crosses of C57BL/6J and CAST/EiJ to tease apart cis and trans contri- butions to gene expression differences between the parental strains. Similarly, (8) performed such an RNA-seq analysis using Drosophila lines. Formally, the idea of using crosses to study cis and trans contributions to differences in gene expression between strains or species is as follows: consider a gene with ex- pressionXP 1 in a homozygous strain 1,XP 2 in a homozy- gous strain 2, and expression XH1 for the haplotype from strain 1 in the F1 cross of 1 and 2, and expression XH2 for the haplotype from strain 2 in the F1 cross of 1 and 2. Let RP =log ( XP 1 XP 2 ) and RH =log ( XH1 XH2 ) . That is, RP corresponds to the log fold-change difference in expression between the two parental strains, and RH to the log-fold- change difference between the expression of the hybrid hap- lotypes in the F1 offspring. The connection between RP , RH, and regulation is as follows: consider that a gene can be regulated via cis, trans, or both (Fig. 1A). Gene expres- sion measurements in the parents and hybrid (Fig. 1B), can be used to infer the nature of regulations underlying the differ- ence in expression in the parental strains (Fig. 1C). Specifi- cally, amending the classification of (6), we have: • conserved: No change in gene expression indicating there has been no change in regulation, i.e. RP = 0 andRH = 0, which impliesRP−RH = 0. • cis: The relative difference in gene expression between the parents is the same as between the haplotypes in the hybrid indicating that the difference in parents is due to local cis effects, i.e.RH̸= 0,RP̸= 0andRP−RH = 0, which impliesRP =RH, arises from changes only in cis-regulatory elements. • trans: Gene expression from the two haplotypes in the hybrid is the same, indicating that differences between the parents resulted from non-local trans regulation, i.e.RH = 0andRP−RH̸= 0, which impliesRP̸= 0, arises from a change only intrans-regulatory elements. • cis + trans: RH ̸= 0and RP−RH ̸= 0with sgn(RH) = sgn(RP−RH) arises as a result of change in both cis- and trans-regulatory elements with changes in cis and trans contributing to changes in gene expression between strains in the same direction. • cis ×trans: RH ̸= 0and RP−RH ̸= 0with sgn(RH)̸=sgn(RP−RH) arises as a result of com- pensatory change in both cis- and trans-regulatory el- ements with changes in cis and trans contributing to changes in gene expression between strains in the op- posite direction. This classification corrects (6), which fail to properly as- sign regulation differences to cis + trans when bothRH < 0 andRP < 0. This issue, and the relationships between RP andRH in general, can be visualized as lines and regions in a 2D plot (5), as illustrated in Fig. 1D(i). While this direct representation of RP and RH is useful, a quantitative as- sessment of the gene regulatory modalities reflected in RP and RH requires a biologically meaningful notion of dis- tance between points in Fig. 1D(i). Consider, for example, the situation where XP 1 = 2XP 2 andXH1 = 2XH2, i.e. a 2-fold change in gene expression between the parents due to solely to cis regulation which corresponds to the point (1,1) in Fig. 1D(i). The distance from this point to the origin (con- served), is √ 2, whereas the same 2-fold difference in gene expression in in the parents due solely to trans regulation withRP = 1 andRH = 0 is distance 1 from the origin. This imbalance can be corrected via a linear transformation. Hallgrimsdottir et al. | bioRχiv | July 14, 2024 | 1–6 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 16, 2024. ; https://doi.org/10.1101/2024.07.13.603403doi: bioRxiv preprint Fig. 1. Geometry of parental and hybrid expression ratios. A) The cis, trans, cis + trans, and cis×trans types of regulation. B) Homozygous parents with haplotypesP1 andP2, counts of RNA molecules for a genes (XP 1 andXP 2 respectively), and the haplotypes in anF1 hybrid along with counts for the gene (XH1 andXH2). C) Differences in regulation between the parents is reflected in distinct ratios between countsXP 1,X P 2 andXH1,X H2. D.i) Illustration of how regulation differences emerge in log-fold changesRP andRH, D.ii) Linear transformation ofRP =log2 (XP 1 XP 2 ) andRH =log2 (XH1 XH2 ) to yield orthogonal cis and trans coordinates, D.iii) Representation of proportion cis in real projective space P1.

Geometry. To decouple the effects of cis and trans regula- tion onRP andRH we begin by noting that if the difference in parental expression is solely due to cis regulation, then RP =RH, or equivalently RP−RH = 0 (orange vertical line in Fig. 1D(ii)). If the difference in parental expression is solely due to trans regulation, thenRH = 0 (blue horizon- tal line in Fig. 1D(ii)). Therefore, the transformation from Fig. 1D(ii) to Fig. 1D(i) is obtained by (RP RH ) = (1 1 0 1 ) (RP−RH RH ) . (1) Thus, the inverse transformation from the coordinate system in Fig. 1D(i) to the coordinate system in Fig. 1D(ii) is given by [1 1 0 1 ]−1 = [1 −1 0 1 ] , (2) i.e., the transformation inverts the mapping of the vertical axis to the diagonal cis line. The determination of whether a difference in parental gene expression is due to cis or trans can now be understood to be an assessment of whether the line passing through the origin and a point (RP,RP−RH) is a perturbation (due to noise) of the line ∆ trans, the line ∆ cis or sufficiently far away from the axes in Fig. 1D(ii) to merit a designation of∆ 2 | bioRχiv Hallgrimsdottir et al. | .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 16, 2024. ; https://doi.org/10.1101/2024.07.13.603403doi: bioRxiv preprint Fig. 2. Contributions of cis and trans regulation to differences in gene expression between yeast strains. A) Results of (9) B) Comparison of the results of (9) to the assignment of genes based on hypothesis testing, C) Reanalysis of the data from (9). cis + trans or ∆ cis×trans (the designations are enumerated in Supp. Table 1). In other words, the sufficient statistic is a point in real projective space P1 (Fig. 1D(iii)), and the pro- portion of the difference in gene expression between parents that can be attributed to cis can be understood to be a scaling of the angle of the line through the origin correspond to the point in P1, i.e. proportion cis = ⏐⏐⏐⏐ 2 πtan−1 ( RH RP−RH )⏐⏐ ⏐⏐. (3) Statistics. The determination of whether a measurement (RP,RH) reflects a difference in gene expression between parents due to cis or trans regulation, or both, requires a sta- tistical assessment (10). Specifically, hypothesis tests can be used to reject a null hypothesis of a difference in gene ex- pression being due solely to trans or solely to cis. Geomet- rically, as evident from the linear transformation on log-fold changes, these two tests correspond to testing whether one can reject the null hypotheses that(RP,RP−RH) is located on thex- andy- axes respectively. We performed such hypothesis tests for data from (9), which consists of bulk RNA-seq performed on hybrid and parental strains of genetically divergentSaccharomyces cere- visiae (see Methods, Supp. Figs. 1,2). Briefly, (9) generated 323 hybrid crosses from 26 parental yeast isolates derived from diverse environmental conditions. In (9), genes were classified according to the representation shown in Fig. 1D(i). First, both RP and RH were tested for statistically signifi- cant differences from 0 and assigned "null" (our conserved) if both tests failed to reject the null hypothesis; then, changes between parental and hybrid allelic ratios were tested for significance. For genes that passed significance thresholds, regulatory assignments were made by dividing the untrans- formed 2D plane into cones (Fig. 2A). We reassigned genes based on our hypothesis tests as derived from the transformed coordinate system (Fig. 2C), with one test with the null hy- pothesis of cis regulation and one with the null hypothesis of trans regulation, thereby putting the two regulation strategies on equal footing. Compared to the original study, we found major differences in assignment (Fig 2B, Supp. Fig. 3A,3B, Supp. Table 2). Interestingly, whereas (9) conclude that "the transcrip- tome is globally buffered at the genetic level mainly due to trans-regulatory variation in the population", we find that a considerable amount of difference in gene expression can be attributed to cis (supp. Fig. 3C), with the difference Hallgrimsdottir et al. | bioRχiv | 3 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 16, 2024. ; https://doi.org/10.1101/2024.07.13.603403doi: bioRxiv preprint Fig. 3. Proportion cis. A) Comparison of slope vs. angle in determining proportion cis. B) Ratio between slope and angle. C) Examples of four genes displaying high variance in proportion cis across cell types. due to (9) deriving assignments in the untransformed coor- dinate system using a statistical testing procedure that treats cis and trans regulation asymmetrically. Specifically (9) re- port 57,253 cases where gene expression difference is due to trans regulation (Supp. Fig. 3A). We find 51,627 cases that can be assigned to trans (Supp Fig. 3B). These numbers are similar; however, (9) report 2,804 cases assigned to cis (Supp. Fig. 3A), whereas we find 17,112 (Supp.‘Fig. 3B). Furthermore, we find 4,807 cases assigned to cis + trans (Supp. Fig. 3B) versus 1,727 in (9) (Supp. Fig. 3A). Overall, there is a marked difference between our results and those of (9) (Supp. Fig. 3C). Proportion cis. In addition to naturally revealing the ap- propriate hypothesis tests to conduct for attribution of gene expression difference in parents to cis or trans, as previ- ously discussed the linear transformation we propose leads directly to a meaningful measure of the proportion of dif- ference in gene expression that can be attributed to cis reg- ulation (or trans) (Equation 3). To illustrate this, we ex- amined and re-analyzed a dataset of gene expression from human-chimpanzee hybrids (11), calculating the proportion cis according to Equation 3. In (11), the proportion cis was calculated using slope, as is natural to do when working in the untransformed coordinate framework, whereas the cor- rect calculation (Equation 4) uses angle in the transformed coordinate system. While the absolute differences are small (Fig. 3A), with a maximum difference of 0.045 (see Meth- ods), the relative difference is large when the proportion cis is small (Fig. 3B), and can be as high as 57% (see Methods). Moreover, our computation provides a biologi- cally interpretable measure of proportion cis. We found sev- eral genes in (11) exhibiting high variance in proportion cis (Fig. 3C), and identified interesting differences between cell types (Supp. Fig. 4 – 7).

The use of crosses between strains to identify the nature of differential regulation is a powerful tool for genetics stud- ies that is particularly relevant now that single-cell RNA-seq can be used for cell type resolution. Moreover, whereas orig- inal studies were limited by gene expression measurement technologies to a handful of genes, genome-wide single-cell RNA-seq assays can complement genome-wide eQTL stud- ies. We have shown that geometric considerations reveal the need for applying a linear transformation prior to visu- alisation and data analysis. The linear transformation also highlights independent axes that lead naturally to hypothesis tests for classifying genes according to the type of regulation underlying differences in gene expression between parents. While we have focused on explaining differences between 4 | bioRχiv Hallgrimsdottir et al. | .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 16, 2024. ; https://doi.org/10.1101/2024.07.13.603403doi: bioRxiv preprint parents that fall into five categories (conserved,cis, trans, cis + trans, cis×trans), our approach can be extended to finer classifications such as in (9). We note that in the hypothesis testing framework, referring to a gene as having differences in expression explained by cis or trans is technically incor- rect. This is because the rejection of the cis null hypothesis only shows that the difference in gene expression is not due solely to cis regulation. It does not mean that the difference in gene expression can, or should be, attributed solely to trans. The same is the case for rejecting the trans hypothesis. In Fig. 2, our coloring of genes as cis or trans is therefore not precise, but we have done so to facilitate comparisons to prior work. Our statistical tests also depend on the assumption that read counts are binomially distributed. This assumption is standard in the absence of biological replicates, however in practice biological replicates should be obtained, and they can and should be used to better assess the extent of technical variation. An extension of this work to that case is possible and can be based on modeling counts with negative binomial distributions, as is done in methods such as W ASP (12). Finally, we note that our framework is general and can be applied to phenotypes other than gene expression. For example, with single-cell RNA-seq data it could be used to assess the regulation mechanisms underlying differences in variance in gene expression. Such extensions will be partic- ularly interesting to explore in conjunction with complemen- tary modalities (13). Data and Code Availability All code to download data and generate the main and supplementary figures is available at https://github.com/pachterlab/HCP_2024, executable in Google Colaboratory notebooks.

IH and LP were funded, in part, by NIH 5UM1HG012077- 02. MC was funded by a NSF graduate research fellowship under Grant No. 2139433. This work was additionally sup- ported by the Caltech Bioinformatics Resource Center.

1. François Jacob and Jacques Monod. Genetic regulatory mechanisms in the synthesis of proteins. Journal of molecular biology, 3(3):318–356, 1961. 2. Rita JS Phillips. A cis-trans position effect at the a locus of the house mouse. Genetics, 54 (2):485, 1966. 3. John S Kovach, Antonio O Ballesteros, Marilyn Meyers, Marco Soria, and Robert F Gold- berger. A cis/trans test of the effect of the first enzyme for histidine biosynthesis on regula- tion of the histidine operon. Journal of Bacteriology, 114(1):351–356, 1973. 4. Christopher R Cowles, Joel N Hirschhorn, David Altshuler, and Eric S Lander. Detection of regulatory variation in mouse genes. Nature genetics, 32(3):432–437, 2002. 5. Patricia J Wittkopp, Belinda K Haerum, and Andrew G Clark. Evolutionary changes in cis and trans gene regulation. Nature, 430(6995):85–88, 2004. 6. Christian R Landry, Patricia J Wittkopp, Clifford H Taubes, Jose M Ranz, Andrew G Clark, and Daniel L Hartl. Compensatory cis-trans evolution and the dysregulation of gene expres- sion in interspecific hybrids of drosophila. Genetics, 171(4):1813–1822, 2005. 7. Angela Goncalves, Sarah Leigh-Brown, David Thybert, Klara Stefflova, Ernest Turro, Paul Flicek, Alvis Brazma, Duncan T Odom, and John C Marioni. Extensive compensatory cis- trans regulation in the evolution of mouse gene expression. Genome research, 22(12): 2376–2384, 2012. 8. Andreas Massouras, Sebastian M Waszak, Monica Albarca-Aguilera, Korneel Hens, Wiebke Holcombe, Julien F Ayroles, Emmanouil T Dermitzakis, Eric A Stone, Jeffrey D Jensen, Trudy FC Mackay, et al. Genomic variation and its impact on gene expression in drosophila melanogaster. PLoS genetics, 8(11):e1003055, 2012. 9. Andreas Tsouris, Gauthier Brach, Joseph Schacherer, and Jing Hou. Non-additive genetic components contribute significantly to population-wide gene expression variation. Cell Ge- nomics, 4(1), 2024. 10. C Joel McManus, Joseph D Coolon, Michael O Duff, Jodi Eipper-Mains, Brenton R Grav- eley, and Patricia J Wittkopp. Regulatory divergence in drosophila revealed by mrna-seq. Genome research, 20(6):816–825, 2010. 11. Kenneth A Barr, Katherine L Rhodes, and Y oav Gilad. The relationship between regulatory changes in cis and trans and the evolution of gene expression in humans and chimpanzees. Genome Biology, 24(1):207, 2023. 12. Bryce Van De Geijn, Graham McVicker, Y oav Gilad, and Jonathan K Pritchard. Wasp: allele- specific software for robust molecular quantitative trait locus discovery.Nature methods, 12 (11):1061–1063, 2015. 13. Ibai Irastorza-Azcarate, Alexander Kukalev, Rieke Kempfer, Christoph J Thieme, Guido Mastrobuoni, Julia Markowski, Gesa Loof, Thomas M Sparks, Emily Brookes, Kedar Nath Natarajan, et al. Extensive folding variability between homologous chromosomes in mam- malian cells. bioRxiv, pages 2024–05, 2024. 14. Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Sté- fan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey,˙Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Hen- riksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Al- gorithms for Scientific Computing in Python. Nature Methods, 17:261–272, 2020. doi: 10.1038/s41592-019-0686-2. 15. Y oav Benjamini and Y osef Hochberg. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal statistical society: series B (Methodological), 57(1):289–300, 1995. Hallgrimsdottir et al. | bioRχiv | 5 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 16, 2024. ; https://doi.org/10.1101/2024.07.13.603403doi: bioRxiv preprint

Acquisition and preprocessing of data from yeast crosses. Processed RNA-seq data from (9) were obtained from http://1002genomes.u-strasbg.fr/files/Diallel_RNAseq/ASE in the file ‘Datafile2_ase_sum_20230609.tab.’ This file included reported values for allelic expression in parents and expression per parental allele in hybrids for 323 unique parent-hybrid trios. Details on the estimation of allelic expression are described in the original study (9). As integer value counts are required for the binomial and the two sample binomial ratio statistical tests, reported values were rounded to the nearest integer before statistical tests were performed. In (9), 7 categories were considered (reverse, attenuating, reinforcing, compensatory, cis only, trans onlyand null). In order to maintain consistency with (6), we collapsed these into four categories: cis, trans, cis + trans, and cis×trans as follows: reverse→cis×trans, cis only→cis, trans only→trans, null→conserved, and attenuating, reinforcing and compensatory →cis×trans if sgn(RP−RH) = sgn(RH) and cis + trans otherwise. Acquisition of data from human-chimpanzee crosses. The reportedlog2 fold changes between human and chimpanzee parental cell lines (RP ) and hybrid alleliclog2 fold changes for human-chimpanzee hybrid cell lines (RH) were obtained from (11). The reportedlog2 fold changes were used without modification to recalculate proportion cis as proportion cis = ⏐⏐ ⏐⏐ 2 πtan−1 ( RH RP−RH )⏐⏐ ⏐⏐. (4) The provided data also included reported proportion cis per gene per cell type, defined as proportion cis (11) = |RH| |RH|+|RP−RH|, (5) which were directly compared to Equation 4. The maximum difference between x := |RH| |RH|+|RP−RH|and 2 πtan−1 ( RH RP−RH ) is 0.045atx = 0.239andx = 0.761. The maximum of the ratio is 1.57and occurs atx≈0.00006. Statistical tests. Binomial test. To test the null hypothesis that there is no difference between expression of parental alleles in hybrids (or that the difference in regulation is purely trans) was performed on rounded integer counts from (9) using the function scipy.stats.binomtest (14). This tests for the probability of having observed a value at least as extreme as k suc- cesses given probabilityp of success and a total ofN trials by summing over binomial probabilities: P(k;N,p ) = (N k ) pk(1−p)N−k. (6) If XH1 is the allelic expression in the hybrid of one parental allele and XH2 is the allelic expression in the hybrid of the other parental allele, the binomial test was performed with k =XH1,N =XH1 +XH2,p = 0.5and a two-sided alternative hypothesis. Two sample binomial ratio test. To test the null hypothesis that there is no difference in the ratio of expression of alleles in parents and allelic expression in hybrids (or that the difference in regulation between parents is purelycis), we performed a two sample binomial ratio test in which we assumed both XP 1 andXH1 are sampled from a binomial distribution with the same probability of successps. Using ps = XP 1 +XH1 XP 1 +XP 2 +XH1 +XH2 , (7) we calculated the grid of probabilities over possibleP1 andH1 values: P(XP 1 =xP 1,XH1 =xH1;NP,NH,ps) = (NP xP 1 ) pxP 1 s (1−ps)NP−xP 1· (NH xH1 ) pxH1 s (1−ps)NH−xH1, (8) whereNP =XP 1 +XP 2 is the total number of counts from parents andNH =XH1 +XH2 is the total number of counts from the hybrid. We then calculated the probability of having observed values at least as extreme as the observedXP 1 andXH1. To account for multiple testing, we corrected both binomial and two sample binomial significance values using the Benjamini-Hochberg correction for false discovery rates (15). We rejected the null hypothesis for tests with false discovery rates less than 0.05. 6 | bioRχiv Hallgrimsdottir et al. | .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted July 16, 2024. ; https://doi.org/10.1101/2024.07.13.603403doi: bioRxiv preprint