text:
Article | Open | Published:

# Gene-by-environment interactions in urban populations modulate risk phenotypes

## Abstract

Uncovering the interaction between genomes and the environment is a principal challenge of modern genomics and preventive medicine. While theoretical models are well defined, little is known of the G × E interactions in humans. We used an integrative approach to comprehensively assess the interactions between 1.6 million data points, encompassing a range of environmental exposures, health, and gene expression levels, coupled with whole-genome genetic variation. From 1000 individuals of a founder population in Quebec, we reveal a substantial impact of the environment on the transcriptome and clinical endophenotypes, overpowering that of genetic ancestry. Air pollution impacts gene expression and pathways affecting cardio-metabolic and respiratory traits, when controlling for genetic ancestry. Finally, we capture four expression quantitative trait loci that interact with the environment (air pollution). Our findings demonstrate how the local environment directly affects disease risk phenotypes and that genetic variation, including less common variants, can modulate individual’s response to environmental challenges.

## Introduction

Environmental exposures, coupled with genetic variation, influence disease susceptibility, and deconstructing their respective contributions remains one of the principal challenges in understanding complex diseases1,2,3,4,5,6,7. Individuals with different genotypes may respond differently to environmental variation and generate an array of phenotypic landscape8,9,10,11,12,13,14. Such gene-by-environment interactions are thought to be pervasive and may be responsible for a large fraction of the unexplained variance in heritability and disease risk9,15,16. Yet, disease risk, owing to either environmental exposures and/or their interactions with genotype, remains poorly understood2,17,18.

Canada’s precision medicine initiative, the Canadian Partnership for Tomorrow Project (CPTP: http://www.partnershipfortomorrow.ca) is a cohort comprising over 315,000 Canadians, and captures over 700 variables, ranging from longitudinal health information to environmental exposures, to determine genetic and environmental factors contributing to chronic disease. The program includes the Quebec regional cohort, CARTaGENE, which has enrolled over 40,000, predominantly French-Canadian (FC) individuals between 40 and 70 years of age19,20,21, to date. Drawing from this founding population of individuals with largely French ancestry, we selected 1007 individuals to determine mechanisms by which genomes, the environment, and their interactions contribute to phenotypic variation. After attributing a regional and/or continental ancestry to each individual using genome-wide polymorphism data, we are able to capture the effect of different environmental exposures on gene expression and health-related traits, while simultaneously controlling for genetic relatedness and migration. Further, in order to capture gene-by-environment interactions through eQTL analyses, we combine whole-transcriptome RNA-Sequencing profiles with whole-genome genotyping and extensive fine-scale environmental exposure data.

## Results

### Population history reveals a fine-grained regional structure

Individuals selected for analyses include those living across different regions in Quebec: Montreal, the largest urban center in the Quebec province (MTL, 4500 individuals/km2); Quebec City, a smaller urban center (QUE, 1140 ind/km2); and Saguenay-Lac-Saint-Jean, a less urbanized region (SAG, 800 ind/km2). Differences in the regional environment within and across these cities, including ambient pollutant concentrations, are known to be associated with various health outcomes22,23. The majority of the Quebec population is of FC descent; a group of individuals descending from French settlers that colonized the Saint-Lawrence Valley from 1608 to the British conquest of 175924. Despite considerable expansion, the population remained linguistically and religiously isolated while remote regions were colonized by small numbers of settlers, such as SAG25,26 and contributed to the establishment of subpopulations. These sequential population bottlenecks impacted the genome of FCs through increasing the relative deleterious mutations load27, while reducing overall genetic diversity in the population relative to the European population28. Using high-density whole-genome genotyping assays (Illumina Omni 2.5), we confirm that FCs (n = 689) form a distinct genetic cluster relative to those of European descent (n = 136) (Fig. 1a, Supplementary Fig. 1a–c), as has been previously observed27. Within this FC group, we capture fine-scale regional genetic variation across Quebec (Fig. 1c, b and Supplementary Fig. 1d), consistent with Quebec settlement history and local ancestry.

### Ancestry contributes marginally to regulatory variation

We were particularly interested in the extent to which individual regional-ancestry and regional-environment account for transcriptional variation in the Quebec population. In an attempt to reduce batch effects in our RNA sequencing experiment, the sampling protocol was standardized across all clinics and all manipulations were performed in the same laboratory. Furthermore, participant’s fasting blood samples were collected by CARTaGENE between 9 a.m. and 11 a.m. Individuals were randomized across sequencing lanes to reduce false associations with traits owing to sequencing differences across lanes. Corrections to mitigate remaining batch effects, unwanted technical and biological variation in gene expression were applied (Supplementary Fig. 4)29 (Methods). Using whole-genome genotyping, we are able to distinguish between “FC-locals” and “FC-regional migrants” (Fig. 1b, Supplementary Fig. 1d, Supplementary Table 2). We define “FC-locals” as individuals of regional ancestry identical to the region they reside in and “FC-regional migrants” as FC immigrants from a different regional ancestry. Among FC-locals, an increasing number of genes are significantly differentially expressed between Mtl- vs Que-locals, n = 505, Que- vs Sag-locals n = 2167, up to n = 6649 and Mtl- vs Sag-locals (Fig. 2a) (value < 0.05/15,632, log-fold change > 0.5). Additionally, a greater number of genes are differentially expressed between individuals having the same regional ancestry but who reside in different regions (FC-locals vs FC-regional migrants with the same genetic ancestry, but residing in different regions), and we find this pattern in nearly all pairwise comparisons of this nature (Fig. 2b). On the other hand, when we performed comparisons between FC-locals and FC-regional migrants, we find very few differentially expressed genes in nearly all comparisons (Fig. 2c).

We replicate these findings by performing comparisons of Europeans and FC-locals residing within the same region and find very few differentially expressed genes between them (Fig. 2d, Supplementary Fig. 5). The lack of differentially expressed genes is not attributable to differences in statistical power as we are able to identify up to 75% of our differentially expressed genes using only 30% of our FC individuals (n = 200) (Supplementary Fig. 6). Furthermore, results are consistent after performing differential expression analyses between regions using a resampling-based method (1000 replicate permutations for each pairwise comparison between regions), thus reducing the possibility that undetected sampling differences between regions, or outlier individuals, drive those patterns. Differentially expressed genes between regions are enriched for genes implicated in oxygen and gas exchange, G-protein-coupled receptors, and inflammatory response (Supplementary Fig. 7, Supplementary Table 3). Although we initially captured both genotypic and transcriptional variation correlated with geographic structure among the FC subpopulations, these results indicate that shared regional environmental exposures influence peripheral blood expression profiles to a greater extent than regional or local (and continental) ancestry, and point to potential critical exposures contributing to pathways, phenotypic variation, and possibly disease development.

### Environment shapes regulatory profiles and clinical traits

To test whether environmental exposures contribute to the geographic variation associated with transcriptional profiles and clinically relevant phenotypes across Quebec, a large collection of fine-scale environmental data (Supplementary Fig. 8 and 9, Supplementary Table 4): satellite-land-use regression models (particulate matter 2.5 (PM2.5) and nitrogen dioxide (NO2)), community land-based measures (ozone (O3) and sulfur dioxide (SO2) for air pollution) are collated. Community level estimates of socio-economic indices (social and material deprivation, population density), and built environment features (greeness, food availability, walkability, park density, street network) are also incorporated. A total of 12 environmental exposures, all of them measured or estimated at the level of three-digit postal code (these areas in Canada include several houses or a neighborhood, and their sizes are inversely proportional to population density) (Supplementary Table 4), are included. Indeed, we observe that these environmental exposures capture broad environmental correlates and variance across the Quebec province (Supplementary Fig. 9). In an attempt to pinpoint if particular environmental exposures contribute more to the gene expression differences across regions, we use this fine-scale information for the analytical treatment of individual exposures specific to the individual and ignore broader geographic sampling categories (i.e., regions).

We find that the expression profiles of differentially expressed genes between regions are largely associated with gradients of annual ambient air composition across Quebec (Fig. 3, Supplementary Fig. 9). A north–south urbanization gradient is indeed reflected by higher annual concentrations of PM2.5 and NO2 in downtown Montreal (data derived from satellite-land-used regression models), however, higher concentrations of SO2 and O3, (land-based measures) are observed in SAG (Supplementary Figs. 8 and 9). The higher annual concentrations of SO2 in SAG, a smaller urban center, are related to the presence of several large industrial complexes22,30. It is widely known that ambient air pollution covaries with season, and we account for blood collection date in our models (Supplementary Fig. 4c). However, we cannot fully exclude a possible residual contribution of season on gene expression patterns.

We apply coinertia analyses31 (CoIA) to our multidimensional data to capture associations between 57 clinical endophenotypes (Supplementary Table 5), environmental exposures (Supplementary Fig. 12), and expression levels of differentially expressed genes and their regulators (Supplementary Fig. 10). All phenotypes are standardized health tests captured by CARTaGENE, and all self-reported disease diagnostics were cross-validated with electronic health records of the participants19. Consistent with previously documented effects of air pollution on cardiac and respiratory traits32,33, we find that arterial stiffness measures, asthma and stroke prevalence, monocytes counts, low-density lipoprotein (LDL), respiratory function (FEV1), as well as liver enzyme levels (Alanine aminotransferase level (ALT), aspartate aminotransferase level (AST), and gamma-glutamyl transferase (GGT)) show the strongest associations with annual SO2 and O3 ambient levels (Supplementary Fig. 10). In our cohort, the gradient of SO2 exposure is associated with detectable detrimental effects on cardio-respiratory phenotypes, more so than ambient annual PM2.5 and NO2 levels (Supplementary Fig. 10), and is the environmental variable that has the highest replicability of association with gene expression (Fig. 3). As a result of these strong associations of annual SO2 ambient levels with detrimental cardio-respiratory phenotypes, and as O3 levels are more dependent on other various ambient factors (sunlight, other NOx emissions), we focus our high-resolution analyses on the participant’s weekly SO2 exposure.

We use a 2-week exposure to SO2 pollution, obtained from averaging over a 14-day window preceding the time-point of each individual blood sampling (Supplementary Fig. 14). The large temporal fluctuations in weekly SO2 ambient concentrations allow us to include individuals from SAG that were exposed to low levels of SO2 (despite SAG having high annual averages), and MTL individuals exposed to high levels of SO2 (despite MTL having lower annual averages), or vice-versa. In that way, we can single out the effect of the local environment itself, predominantly attributable to SO2 exposure, to the broader regional effect detected earlier. Using a robust resampling approach to balance the number of individuals in each category, we are able to identify with confidence 170 differentially expressed genes between high- and low-SO2-exposure individuals; these are also found to be differentially expressed between regions (Fig. 2a, Supplementary Table 6).

Furthermore, while multivariate models show that gene expression variation for those 170 genes is significantly associated with 2-week SO2, they do not show an association with smoking, socioeconomic status, or with most built environment characteristics (Supplementary Table 6). We perform a sensitivity analysis using MTL-only samples, thereby removing the potential influences of geographic region and regional ancestry. We replicate these associations with pollution, and the lack thereof, for smoking and socio-economic status (Supplementary Table 6). These results indicate that the regional effect on the gene expression is mostly associated with ambient air pollution, and less so, or not at all, with diseases, smoking, or the socio-economic factors that were measured. Those 170 differentially expressed genes are again enriched in oxygen-transport activities, and in several pathways involved in leukocyte migration during chronic inflammation, including CXCR chemokine activity and G-protein-coupled receptors (Supplementary Table 6). Circulating blood leukocytes can migrate to sites of tissue injury by responding to proinflammatory cues and are known to migrate through the blood flow to lung epithelial cells during inflammatory response34.

To disentangle the effects of SO2 exposure from the effects of region on gene expression, we conducted a sensitivity analysis and show that not only is this pattern observed across the whole Quebec province, but it also replicates within Montreal (Supplementary Table 7), suggesting that SO2 exposure, rather than the region itself, is modulating these associations. We find that the expression of the 170 DEGs (between high- and low-exposure to SO2) is also associated with four key clinical traits (Forced expiratory volume (FEV1), lung disease, liver enzymes, and arterial stiffness) (Supplementary Fig. 12). Additionally, when the effects of these four clinical traits are regressed out from gene expression, SO2 exposure remains significantly associated with gene expression (Supplementary Table 7). This suggests that SO2 exposure itself modulates some of the variation in gene expression, and this variation is not only associated with the underlying health status.

The four clinical traits that were found to be associated with differential gene expression (FEV1, lung disease, liver enzymes, and arterial stiffness), are consistently reported as influenced by air pollution by other studies35,36,37,38,39. Chronic diseases developing from these detrimental endophenotypes (asthma and cardiovascular diseases) are well documented to be associated with air pollution levels12,13,14,22,37. Gamma-glutamyltransferase (GGT) has been reported to occur in atherosclerotic plaques40, is elevated following pollution exposure41,42, and is predictive in a dose-dependent manner of cardiovascular risk43. Interestingly, we find GGT levels to be associated with the differentially expressed genes across SO2 exposure environments, in particular those genes enriched in blood coagulation and platelet regulation (Supplementary Fig. 12), Collectively, these results reveal associations between environmental pollutants, endophenotypic traits, as well as transcript levels, and that the type and direction of associations are consistent with detrimental effects of air pollution, or a correlated variable, on health status.

### Environment modulates the penetrance of genetic variants

Environmental factors not only directly affect phenotypic variation, but can also modulate associations between segregating genetic variants and phenotypes1,44,45. To discover gene-by-environment interactions in both FCs and Europeans, we identify eQTLs for which the effect size is modulated by exposure with one of four ambient air pollutants (env-eQTLs): PM2.5, NO2, O3, and SO2. First, we identify canonical eQTLs using 5,313,384 genotypes and show a high replication for proximal canonical eQTLs (cis-eQTLs) with previously discovered cis-eQTLs (Supplementary Table 8).

To identify gene-by-environment interactions with air pollution (env-eQTLs), we use a randomly generated discovery cohort (n = 416) to perform regressions of gene expression levels (eGenes) on cis-SNPs (eSNPs), pollution level, and the interaction between eSNP and pollution (see Supplementary Fig. 15, for a schematic representation of the procedure/design). We use a four-step process that accounts for multiple testing: (1) we compute Bonferroni-corrected p values, adjusting for the number of eSNPs tested for each gene, (2) we retain the lowest Bonferroni-corrected p value for each eGene and transform this set into q values46 to determine statistical significance (FDR < 0.05, to correct for the 15,632 total genes tested in the cohort). This results in the identification of ten unique significant eSNP–eGene pairs (with nine unique eGenes). (3) We then examine these candidate pairs in our replication cohort (n = 417), where four out of the ten pairs are significantly replicated (q value < 0.05) with the same direction of effect in both the discovery and replication cohorts. Last, (4) all four replicated eSNP–eGenes associations (eGenes, n = 3; eSNPs, n = 4) remain significant using empirical p value estimates through permutations on the combined cohort (n = 833 individuals) (Supplementary Table 9).

Following the application of this stringent filtering, we identify and replicate three eGenes (four eSNP–eGene pairs) for which air pollution (either PM2.5, NO2, SO2, or O3) modulates the association between the genotype of at least one eSNP and the eGene expression (Fig. 4, Supplementary Fig. 16, Supplementary Table 9). One eGene, atad2, is identified as interacting with both NO2 and SO2 ambient levels. zp3 is a glycoprotein interacting with proteins in the extracellular space. Interestingly, two eGenes are ATPases with epigenetic activities, regulating chromatin structure (smarca2) or assisting in chromatin and histone binding of transcription factors (atad2)47,48.

Among the significant env-eQTLs (FDR q value < 0.05 in discovery and replication cohorts) (Supplementary Fig. 16, Supplementary Table 9), we identify an interaction with NO2 and the SNP–gene pair rs10156534-smarca2 (Fig. 4a). Furthermore, we find a deletion (chr9: 3,177,272) in an enhancer downstream of smarca2 that is significant for an interaction with NO2 levels (Fig. 4c). SMARCA2 protein is part of the large chromatin remodeling complex SNF/SWI (Fig. 4b), and is required for the transcriptional activation of genes repressed by chromatin by mobilizing nucleosomes. The SNF/SWI complex is a tumor-suppressor gene complex and is also required to activate other tumor-suppressor genes. In addition, it has been found to be potentially contributing to a range of inflammatory diseases, including childhood asthma and systolic blood pressure. Interestingly, as discussed above, we find, in CARTaGENE, that spirometry phenotype (FEV1) and arterial stiffness, which are tightly linked to asthma and blood pressure respectively, are associated with differential expression of genes across regional environments. This suggests that environmental differences in air quality may act on the regulation of several genes and pathways and promote pro-inflammatory states which can lead to cardiorespiratory dysfunction.

The eSNP–eGene pair rs62518566-atad2 is an env-eQTL that interacts with SO2 and NO2 exposition (Supplementary Fig. 16d and e). ATAD2 protein belongs to a large family of ATPases that contains a bromodomain; that is, a protein domain that reads epigenetics marks on chromatin and affects gene regulation48. It is a regulator of chromatin dynamics and acts as a co-activator of estrogen and androgen receptors. atad2 is associated with several human diseases, and serves as a marker of poor prognosis in a variety of different cancers49,50.

### Variant frequency and the environmental impact on traits

Allelic frequency has an inverse relationship with phenotypic variation, and, in particular, on eQTLs susceptibility to environmental modifications. First, an inverse relationship between effect sizes on transcript abundances and lead eSNP minor allele frequencies (MAFs) is observed in our cohort (Supplementary Fig. 17). This pattern is consistent with natural selection acting to stabilize gene expression51,52,53. Second, using the estimated correlations from a CoIA analysis between all SNPs in cis of significant eGenes and endophenotypic traits, we test whether the size of the correlations are related to the MAF of the SNP. To do so, we classify the SNPs as common (MAF > 0.1), and less common (MAF between 0.05 and 0.1) and, for each endophenotype, calculate the odds ratios of observing less common variants (compared to common) for stronger endophenotypic associations. We find that less common variants are overrepresented for stronger associations between eSNPs and some endophenotypic traits (Supplementary Fig. 18). More specifically, respiratory (Asthma and FEV1/FVC ratio) and cardiovascular (Stroke, peripheral AIX) phenotypes show larger changes in values in individuals with less common variants at env-eQTL loci (Supplementary Fig. 18). These results suggest that SNP allele frequency is negatively correlated with endophenotypic trait changes when influenced by environmental perturbations, which is coherent with previous studies and theoretical predictions53,54.

Our findings illustrate that the impact of the geographic region of residence on the blood transcriptome overrides that of ancestry. Moreover, ambient air pollution exposures are likely contributing to this regional effect in Quebec and may explain the differences in some clinical traits among regions such as asthma prevalence. Fortunately, in Quebec, and in many parts of the developed world, air quality has improved since the 1980s30,55. However, there has been a sharp increase in anthropogenic pollution levels in many parts of Asia caused by the rapid industrialization and increased use of fossil fuel energies. In the context of global climate change, air pollution and hazardous air quality events are predicted to become more frequent and cause additional morbidity and mortality23. More broadly, our work shows how environmental exposures modulate gene expression directly, can act upon the penetrance of genetic variants, and can affect clinically relevant phenotypes in humans.

## Methods

### Contact for reagent and resource sharing

Further information and requests for reagents may be directed to the Biobank CARTaGENE which regulates the access to the data and biological materials (http://www.cartagene.qc.ca/en/contact-us).

### Study population

The study protocol was approved by the Ethical Review Board Committee of Sainte-Justine Research Center and all participants provided informed consent. CARTaGENE biobank comprises more than 40,000 participants aged between 40 and 60 years, recruited at random among three urban centers in the province of Quebec. CARTaGENE is a regional cohort within the Canadian Partnership for Tomorrow Project, including over 315,000 participants, with various measures obtained from blood parameters, biological function, disease history, lifestyle, and environmental factors19.

### Sample selection

For freeze 1, we selected 708 individuals from the CARTaGENE’s biobank samples with available Tempus Blood RNA Tubes (ThermoFisher Scientific) and Framingham risk scores, ensuring an equal representation of ages and gender. Two-hundred-and-ninety-two additional samples were subsequently selected from CARTaGENE (freeze 2) based on their RNA and complete arterial stiffness (AIx) measures availability. These samples were selected for having high AIx values as well as average AIx values to complement the first freeze of samples with the intention of achieving a broad range of arterial stiffness values across the complete study cohort. All samples were collected in the same year, with a standardized protocol in all sampling clinics19. All blood samples were collected in the morning, on fasting participants.

### Genotyping and QC

In total, 928 samples with RNA-Seq profiles that passed quality control (QC) thresholds were genotyped on the Illumina Omni2.5 array to obtain high-density SNP genotyping data. A total of 1,213,103 SNP were retained after filtering and QC (Hardy–Weinberg p value > 0.001, MAF > 5% and percent of missing data < 1%).

### RNA sequencing

Whole blood samples were collected from participants in 2010. Total RNA was isolated using the Tempus Spin RNA isolation kit (ThermoFisher Scientific) and a globin mRNA-depletion was performed using the GLOBINclear-Human kit (ThermoFisher Scientific). The quality and integrity of the RNA samples were verified using an Agilent Bioanalyzer 2100 and all samples had an RNA integrity number (RIN) > 7.5. A RIN above 7.5 is indicative of high quality RNA in the sample and for which RNA degradation is minimal, indicating optimal transport and preservation conditions. Our RIN threshold is more stringent than other large-scale consortium studying gene expression in tissues51,56. TruSeq RNA Sample Prep kit v2 (Illumina) was used to construct paired-end RNA-Seq libraries with 500 ng of globin-depleted total RNA. Recommended Illumina protocols were followed for quantification and quality control of RNASeq libraries prior to sequencing. Paired-end RNA sequencing was performed on a HiSeq 2000 platform at the Genome Quebec Innovation Center (Montreal, Canada). Sequencing was performed for freeze 1 (708 samples) using three samples per lane, and for freeze 2 (292 samples) using six samples per lane yielding about 60 million reads per samples. All RNA-seq experimental steps following blood draw were conducted in the same central laboratory, and samples were distributed randomly over sequencing lanes (Supplementary Fig. 3a, b), thereby reducing the introduction of experimental bias at these steps.

Reads were trimmed for adapters and bad quality bases first using Trim Galore and were then assembled to a reference genome (hg19, European Hapmap (CEU) Major Allele release) using STAR (v2.3.1z15)57 using the two-pass protocol, as recommended by the Broad Institute. The two-pass protocol consists in two consecutive mappings steps having the same set of parameters with only the reference that is optimized in the second mapping procedure. The first mapping is done using the reference gene definition coming from ENSEMBL (release 75). Then, using the splicing junction database files formed by the first pass mapping step for all the samples combined together and the same gene definition file, a second reference is indexed and optimized and is used for the second mapping step. The number of mismatches allowed across pair is five and a soft-clipping step that optimizes alignment scores is also done automatically by STAR. The PCR duplicates were conserved as it was shown that quantification of highly expressed genes were disproportionately affected by PCR duplicates removal58. Only properly paired reads were kept (using samtools59) for the analysis, according to STAR parameters. After these steps, HTseq (v.0.6.1p1)60 was launched separately on each alignment file using the same gene reference file that was used for the alignments.

All analyses downstream were conducted using R 3.1.2 and R 3.2.2 and Bioconductor R packages.

### Fine-scale population genetic structure within Quebec

To unveil finer scale patterns of population structure, i.e., differences between individuals with European ancestry versus individuals having a French Canadian ancestry, we also used ChromoPainter (v0.04)61, a haplotype-based method powerful enough to detect fine-scale genetic structure. Original genotyping data was used apart from singletons, yielding to 1,908,336 SNPs. Singletons were removed as they are non-informative for phasing and contribute to computation burden for the step of haplotypes sharing inference performed with ChromoPainter. Genotypic data was phased with SHAPEIT (v2.r644)62 using the HapMap genetic maps. Coancestry matrices were obtained from ChromoPainter with parameters estimation step done with ten iterations on four chromosomes only. ChromoPainter method performs a reconstruction of every individual genome using chunks of DNA donated by the other individuals and report matrices of the number and length of those chunks. We used the chunk count matrix to (1) run FineSTRUCTURE algorithm to build a tree (as recommended for large data set, we performed 10,000,000 burn-in and runtime MCMC iterations) (Supplementary Fig. 1D) and to (2) perform a PCA (Fig. 1a, Supplementary Fig. 1c). Regional ancestry for each FC was determined based on the three clusters obtained from the fineSTRUCTURE tree, (Supplementary Fig. 1d, Fig. 1b).

In agreement with Quebec settlement history, previous studies of the Quebec population28,63, and the fineSTRUCTURE tree, a PCA of FC individuals reveals groupings of sub-populations of individuals that follow a North–South structure (Fig. 1b, c). The founding event from French settlers followed by the subsequent colonization of remote regions has led to population differentiation among regions in Quebec28,63. By further restricting the group of individuals to be analyzed to only FC (n = 726) and considering their region of residence (either Quebec City, Montreal, and Saguenay) a PCA on the chunk count matrix reveals three groups corresponding to region of residence, with the Montreal and Quebec groups overlapping to a greater extent, in line with their greater geographic proximity (Fig. 1b, c). Those three groups were also recovered by the fineSTRUCTURE tree (Supplementary Fig. 1D). Considering all SNPs and the whole haplotypic structure is the key in seeing differences for those two metropolitan regions that have low differentiation. We further identified several participants with a regional ancestry discordant with their region of residence: an indication of recent regional migration of these participants across Quebec regions (Supplementary Table 2).

### Imputation

To increase the power for the association study with gene expression levels, variant imputation was conducted on 968 individuals for which the genotyping was available from the Illumina Omni2.5 array. We pre-phased the genotypes with SHAPEIT (v2.r64410)62 using the default parameters, on both the autosomes and the chromosome X. We filtered variants for MAF > 1% and Hardy–Weinberg value > 0.0001 and passed the haplotypes to IMPUTE2 (v2.2.2)64 to perform the imputation using the 1000 Genomes Phase I integrated haplotypes (Dec 2013). We used the parameters Ne = 11418 and call thresh = 0.9. We removed variants with a call rate <90%, MAF > 1%, and Hardy–Weinberg p value > 0.0001. A total of 9,157,622 variants passed the filters. Of these, 8,877,297 variants were found on the autosomes and included 779,579 insertion-deletion polymorphisms (indels) (8.78%) and 8,097,718 SNPs (91.22%). 280,325 variants were found on the chromosome X, which included 28,504 indels (10,16%) and 251,821 SNPs (89.84%).

To determine the ancestry of each individual from genotyping data, we carried out a principal component analysis (PCA) with SNPs pruned for LD (pairwise r2 > 0.2 and 50 SNPs window shifting every five SNPs) (Supplementary Fig. 1A), yielding 146,689 SNPs. The continental ancestry (African/European/Asian/Canadian/American/Middle-Eastern) of each individual was determined based on the PCA plot (Supplementary Fig. 1A) and verified as to whether it corresponds to self-reported ancestry based on the country of origin of four grandparents. If the country of origin of three out of four grandparents and the PCA continental grouping were concordant, the individual was assigned to a continental origin.

### RNA-sequencing filtering

Genes with counts-per-million below 0.5 in more than half of the cohort (505 individuals) were removed from the analysis for a total of 15,632 genes retained for all downstream analyses. Individuals that showed obvious outlier after visual inspection of principal component plots were removed (three individuals).

### Variables contributing to transcriptomic variation

The deep phenotyping of the CARTaGENE cohort allow for a thorough exploration of the biological and environmental factors that may influence genome-wide gene expression patterns. As most statistical procedures assume a normal distribution to the underlying data, we transformed the normalized counts from freeze 1 to a Gaussian distribution using a log2cpm transformation using edgeR. We summarize the gene expression levels by performing a PCA on the normalized expression matrix (ePCA). To identify variables that contribute to genome-wide gene expression variation, we performed a stepwise regression (stepwise search from both directions) on ePC1 and ePC2. Results of the stepwise regression are given in Supplementary Table S1, as well as the results from the replication analyses using freeze 2. We included the following low level endophenotypes in the stepwise procedure: set, region of residence, cell counts (lymphocytes, neutrophils, monocytes), arterial stiffness, age, and sex.

Using the freeze 1 data set of 708 individuals, we quantified the proportion of the variance in expression attributable to cell counts, age, sex, region, and arterial stiffness (Supplemental material) by using principal variance component analysis (PVCA), and found that the region of residence explains 16% of the variance in gene expression, while the effects of age, sex, and cell counts were much lower (Fig. 1e). These analyses were repeated on an additional 289 participants (freeze 2) and both of these effects were found to be replicated on expression profiles (Supplementary Table 1). Similarly, when combining transcriptional profiles for all individuals, we found that the region of residence explains 15% of the variance in gene expression both in FCs and in Europeans (Supplementary Fig. 2).

### Sampling site effect within region

The RNA extractions and library preparation were performed for all individuals in the same laboratory to reduce technical bias. However, participants were sampled across four different sampling sites inevitably situated within geographical regions where participants lived. Our experimental design was built in such a way that sequencing run was not correlated with region of residence (Supplementary Fig. 3a). To evaluate whether the sampling site has any effect on the RNA-Seq quantification data, we performed extensive analyses of the two sampling sites situated within Quebec City: St-Sacrement (STS, n = 136) and Enfant-Jesus (EF, n = 129). QUE individuals expression profiles from the combined data set show that individuals from STS and EF form a single cluster on a ePCA plot (Supplementary Fig. 3b). Furthermore, a variance component analysis (PVCA) was performed on the QUE individuals only and including sampling site as an explanatory variable shows that the sampling site explains <5% of the variance within QUE region, while freeze explains 15%, age 5%, and gender 2.5% (Supplementary Fig. 3c). In comparison, in FCs or Europeans, region of residence accounts for 15% of variance in gene expression. In addition, we performed a differential expression analysis between sites within a region (see details below) using permutations, and found that there are no genes differentially expressed between clinics within a region, supporting the absence of sampling differences between clinics affecting gene expression to a detectable and significant level.

### Correction for technical and biological unwanted variation

High quality RNA-sequencing of all 997 individuals reveals a similar geographic structure in transcriptional profiling than population structure from genotyping (Fig. 1c). Investigation of the variance associated to gene expression reveals that region of residency (variable of interest) explains about 16% (Supplementary Fig. 2a) of the variance regarding the population of origin (Supplementary Fig. 2b, c), but unwanted variables explain a certain proportion of the variance (Fig. 1e, Supplementary Fig. 2b, c).

RNA-Seq data generation, and expression data in general, are prone to technical biases which in some cases can mimic, or be confounded, with biological variation. The appropriate normalization pipeline in an RNA-Seq experiment will depend on the experimental design and the hypothesis being tested. Local sequence context can bias the uniformity of read counts along the genome, and sophisticated normalization pipeline may be necessary when comparing expression levels across genes65. Most experimental designs of RNA-Seq studies, like the one presented here, compares different groups of individuals to each other, therefore the normalization pipeline should rather focus on removing unwanted variation across individuals.

We removed the effects of hidden covariates potentially affecting expression levels using surrogate variable analysis (SVA)29. We used the SVA correction, retaining five surrogate variable, for the differential expression analyses, correcting for technical (i.e., runs, sets, number of reads) and biological (i.e., date of appointment, time of the year, sex, smoking status, cell counts) effects on gene expression (Supplementary Fig. 4). We performed the same stepwise regression approach as previously, but on the SVA corrected expression level matrices and show that we retained the variation associated with region, but removed any effects of cell counts and arterial stiffness that was present in the uncorrected expression levels (Supplementary Fig. 4, Supplementary Table 1). The corrections do not fully compensate for the effect of the freeze (technical), we therefore include this covariate in all subsequent analyses. Estimating the variance associated with hidden batch has been shown to remove variation associated with biological and technical factors and also increase the power to identify eQTLs58,66.

### Differential expression analysis

Because of the large proportion of the variance in gene expression explained by region of residence, we identified genes that are differentially expressed between pairwise comparisons between the FC-locals from the three regions (Montreal, Quebec, and Saguenay). Using edgeR67, we performed a differential gene expression analysis using the 15,632 genes that passed the QC filters established above. We performed the differential expression modeling using the following statistical model:

$${\mathrm{\mu }}_{{{{i}}\mathrm{g}}} = {\mathrm{\beta }}_{\mathrm{g}}{\mathrm{Rr}}_i + {\mathrm{\beta }}_{\mathrm{g}}{\mathrm{Ro}}_{{i}} + {{B}}_{{i}} + {{S}}_{\mathrm{g}} + \epsilon _{{{{i}\mathrm{g}}}}$$

where Rr is the region of residence, Ro the region of origin, B is the surrogate variable, representing the batch effect estimated by SVA, and S represent the freeze effect that is included in the final (see below for further details).

The significance level of the test was estimated as a gene value below the Bonferroni-corrected threshold of 3.20 × 10−6 (0.05/15,632). The SVA corrected expression levels retained the variation associated with region, but removed any effects of cell counts that was present in the uncorrected expression levels (Supplementary Table 1).

We performed a power analysis of our ability to detect differentially expressed genes with smaller samples sizes. Several of our comparisons of regional- or continental-migrants with FC-locals involve smaller number of individuals (Supplementary Table 2). We therefore assessed our ability to detect differentially expressed genes by performing differential expression analyses between groups for which we found large number of differentially expressed genes, but using a smaller subset of random individuals (without replacement) of each of these groups. We randomly selected 15 Mtl-locals and 15 Sag-locals, and performed the DGE analysis using the same model as above. We also performed the analysis using 50 Mtl-locals and 50 Sag-locals. In each case, we could identify differentially expressed genes which largely overlap with the differentially expressed genes detected in comparisons using all individuals (Supplementary Fig. 6). We observe that with an increasing number of individuals, our power to detect differentially expressed genes increases and that the identity of the differentially expressed genes detected in each of these comparisons largely overlap (Supplementary Fig. 6).

We further support the effect of region of residency on gene expression by performing differential gene expression analysis across regions using permutations that are even more robust to batch effects. The permutation-DGE analyses confirm that differences are the greatest between MTL and SAG. Similar permutation analyses also show that individuals living in the same region but sampled in different clinics have similar gene expression profiles (Supplementary Fig. 3B), supporting the absence, if not minor, of effects of sampling procedures on the gene expression across sampling clinics.

### Regional environmental effects on gene expression

We take advantage of the presence of individuals from different regional and continental origins in our cohort to disentangle further the effects of the genetic background and environmental influences on genome-wide gene expression. We first selected individuals of either FC and European continental ethnicity (Fig. 1a, Supplementary Fig. 1). A total of 798 individuals including 136 Europeans and 662 FC were selected for downstream analyses. We stratified the individuals according to their continental origin (FC vs Europeans), and further stratified the FCs into their assigned genetic ancestry (MTL, QUE, SAG) obtained from the fineSTRUCTURE analysis (Fig. 1b, Supplementary Fig. 1D). We then determined their region of residence (MTL, QUE, SAG) for a total of 12 ancestry-residence groups: we identified individuals for which their origin (Continental or regional) is discordant with the region they reside, which we refer to as continental- and regional-migrants respectively (Supplementary Table 2). We also identified FC individuals for which their regional origin is concordant to the region they reside, which we refer to as FC-locals (Mtl-FC-locals, Que-FC-locals and Sag-FC-locals). We performed the differential gene expression analysis pipeline as described above for different pairs of continental-migrants, regional-migrants, and FC-locals to disentangle the effects of the genetic background and the regional environment on genome-wide expression (Fig. 2). We selected 6649 genes that show differential expression (p value < 3.20 × 10−6) in the comparison between Mtl-FC-locals and Sag-FC-locals. Using the 12 origin-living groups and the 6649 genes, we performed an unsupervised clustering and visualized the groupings using a heatmap (Supplementary Fig. 5).

### Gene enrichment and reactome analyses

Gene enrichment analyses were performed using the topGO package in R, with a Fisher exact test. Differentially expressed genes between MTL-locals and SAG-locals were compared against the 15,632 genes expressed in the CARTaGENE cohort that were retained after QC filters (background). Reactome enrichment analyses were conducted with R the package reactomePA, and here again, the background set of genes was defined as the 15,632 genes expressed in blood that pass our filters (Supplementary Fig. 7 and Supplementary Table 3).

### Fine-scale environmental data

We obtained air quality measures in the year of sampling (2010) from either land-based stations (SO2, ozone) or national LUR models estimates (PM 2.5 and NO2) incorporating information from land use data and satellite remote sensing55,68,69,70. Built environment variables (street network, population density, food deserts, greenness, walkability) and social and material deprivation indicators were accessed through the Quebec government data portal (https://www.inspq.qc.ca/environnement-bati). All environmental data sources are described in Supplementary Table 4.

Environmental data was available at the three-digit postal code district level (i.e., Forward Sortation Area, FSA), or was reformatted to this geographic scale. Postal code districts in Canada are small geographic areas which assist in delivering mail. Postal codes are a series of six digits that identify a small geographic area in a municipality, usually grouping just a few houses together or a small neighborhood. Three-level digits are larger areas that include several houses, a small neighborhood, or a small village. The population of FSAs in Canada range from a few hundreds to tens of thousands of individuals. Three-digit postal code districts can be of different areas, and are smaller in densely populated areas, and larger in areas of low population density. Maps in Fig. 1c, d and Supplementary Fig. 9 depict three-digit postal code districts as thin gray lines areas, and each district is colored with the mean value of interest in each map. Each individual in the CARTaGENE cohort has a three-digit postal code district associated to it, referring to the location of its primary residence. We assigned fine-scale environmental measures to each individual based on its three-digit postal code.

### Coinertia analyses

Coinertia analysis (CoIA)31,71 is a multivariate statistical part of the large family of ordination methods, such as PCA, redundancy analysis (RDA), or canonical correlation analysis (CCA). CoIA is a general approach and existing methods such as the ones mentioned above appear as special cases of it31. These methods have been widely used in ecological research, including CoIA which has been more recently developed. Collectively, these methods allow for detecting an underlying data structure between two data tables. CoIA uses a combination of PCA and multivariate linear regressions to detect linear combinations of variables from one data table that explain the variance in the second data table. CoIA is more flexible than RDA or CCA, and overcomes their limitations by allowing for more variables than the number of samples to be tested31,71, which is generally the case in genome-wide scale analyses (i.e., more genes than individuals). This makes CoIA a method of choice to integrate data of diverse types, and of high-throughput like most omics data.

We first used CoIA analysis to reveal the common structure between differentially expressed genes (Fig. 3, Supplementary Fig. 11) and the fine-scale environmental data. We produced two separate principal component analyses (PCAs) based on continuous encoded matrices of both environmental and gene expression levels (normalized for library size and sequencing freeze). The data were centered and reduced to one unit of variance prior performing the PCA analysis. We conserved components for each PCA to explain 80% of the variance in the data. We imputed missing data only for the fine-scale environmental data set (there were no missing data in the gene expression matrix) using the function imputePCA from the R package missMDA. The coinertia analysis performs a double inertia analysis of each data set and then project the variables of the original environmental and gene expression data sets on the new co-inertia axes. Relationships between the two matrices were assessed by comparing the CoIA estimated from the real data set with the CoIA distribution estimated after bootstrapping. Two sets of 500 of CoIAs were computed independently between gene expression and fine-scale environmental data. Supplementary Fig. 11 depicts the resampling scheme. For each Group 1 or Group 2 (n = 497 for each) a total of 10,000× resampling of 200 individuals (without replacement) were performed. We performed a CoIA for each resampling step. We report the median value of the distribution of each environment–gene expression pair cross-tabulated values for each group. Gene enrichment were performed using gProfiler72, and using the 15,632 expressed genes that passed our filters in whole blood as the background gene set (Supplementary Table 3). We evaluated the significance of the correlations between the two matrices with a multivariate generalization of the Pearson correlation coefficient (RV coefficient) using a permutation test (RV-test) with 10,000 steps  from the R package ade4.

To identify clinically relevant endophenotypes that are associated with fine-scale environmental data, we performed a CoIA between 57 clinically relevant endophenotypes (Supplementary Fig. 10) and fine-scale environmental data. The 57 clinically relevant endophenotypes were selected to encompass physical measures (BMI, height, age, sex), most systems relevant to the human health (cardiovascular system, pulmonary functions, hepatic system, renal system, disease history, vision, immune system) and lifestyle measures (smoking status, alcohol consumption, nutrition, physical activity). All biochemical endophenotypes were measured in a single central laboratory. We resampled 10,000 times 493 individuals from the cohort, and performed CoIA at each step between endophenotypes and fine-scale environmental variables. We report the median value of the distribution of each environment–endophenotype pair cross-tabulated values (Supplementary Fig. 10).

To reveal possible associations between expression levels and endophenotypes, we then performed CoIAs with a similar resampling scheme between 12 selected endophenotypes that were the most strongly associated with air pollutants from Supplementary Fig. 10 (Stroke, Arterial stiffness measures, spirometry measures, Asthma, monocyte counts, LDL, AST, ALT, GGT) and differentially expressed genes (results shown in Supplementary Fig. 12).

### Exposure windows of weekly SO2

To increase our resolution in air pollution exposures, we used daily SO2 ambient levels measured in each three-digit postal code. We calculated the average exposure during the 14 days preceding the blood draw for each participant. This way, we reduce the effect of random fluctuations due to technical artifacts or short-term meteorological anomalies that may affect measurements. Also, changes in gene expression and biomarkers in blood following a pollution exposure has been documented as a relatively fast phenomenon, occurring after just a few days of exposure36. We then categorized the participants using a k-means algorithm73 into high exposure or low exposure categories (see details on the number of participants and cluster centers in the eQTL section below).

### DGE between high- and low-SO2 exposure

To find differentially expressed genes between high and low exposure individuals, we used the same approach as described above for identifying differentially expressed genes between regions, with the following modifications: given the unbalanced number of individuals in each category (108 high exposure vs 800 low exposure) of exposure, we resampled 100 times 108 individuals with replacement from the low- and high-exposure category and performed the DGE pipeline. We performed the SVA while retaining variation associated with SO2 exposure. We combined the results of DGE analysis in a list of 468 differentially expressed genes, and from these candidates, 170 genes were also identified as differentially expressed between regions (Fig. 2a). Those strong 170 candidates were used for enrichment, CoIA, and multivariate models. We also identified genes (transcription factors) that regulate our 170 differentially expressed genes (RDEGs) using cytoscape, and we used them in addition to the differentially expressed genes in the CoIA analyses.

### Multivariate models for SO2 exposure

In an effort to characterize the effects of confounding variables on pollution exposure, we applied multivariate models on gene expression levels. First, similar as in the differential gene expression analysis, we performed a SVA to remove unwanted variation of technical or unknown biological variables while retaining the variation around SO2 exposure. We then built multivariate models using the SO2, O3, and PM2.5 14-day exposures, as well as the remaining 9 non-pollution environmental exposures (Supplementary Fig. 9), as well as smoking status. Smoking status may indeed cause similar changes in endophenotypes as pollution exposure. We then selected the endophenotypes revealed by the CoIA as being the most associated with region and pollution exposure (Lung disease, Asthma, Stroke, monocyte counts, liver enzymes (AST, ALT, GGT), Arterial stiffness, spirometry tests, and lymphocyte counts), and tested whether any of these would explain variation in the 170 candidate genes. Furthermore, after having identified the health endophenotypes that are associated with gene expression in MTL and in the whole data set (FEV1, liver enzymes, lung diseases, and arterial stiffness, see Supplementary Figs. 10 and 12), we regressed out their effect from the expression of the 170 candidate genes, and run the multivariate models to test for the effects of environmental variables itself (results collated in Supplementary Table 7).

### env-eQTL analysis

Environmental factors not only directly affect phenotypic variation, but can also modulate associations between segregating genetic variants and phenotypes1,44,45. To discover gene-by-environment interactions, we identified eQTLs for which the effect size is modulated by environmental exposure to one of four ambient air pollutants (env-eQTLs): PM2.5, NO2, O3, and SO2. We categorized the participants using a k-means algorithm73 into two categories, high or low exposure, irrespective of the pollutant type (Table 1). A k-means algorithm attempts to partition the individuals into k groups (here, k = 2), such that the sum of squared Euclidean distances from points to the assigned centroid (cluster mean) is minimized.

We adopted a strategy (Supplementary Fig. 15) to randomly divide the CARTaGENE cohort into discovery and replication cohorts. During this process, for the discovery of eSNP–eGene pairs, we scan the genome at ± 500 kb of the TSS of gene to find all putative cis-eSNPs. We used the following model where gene expression (Y) is regressed on a given SNP (S), a given environmental air pollutant (E) and the interaction between S and E:

$${\mathrm{Model}}:{\mathrm{Y}}_{{{ijk}}}\sim {\mathrm{S}}_{{{ijk}}} + {\mathrm{E}}_{{{ijk}}}{\mathrm{ + S}}_{{{ijk}}}{\mathrm{E}}_{{{ijk}}}$$

The gene expression level, was normalized using an inverse normal transformation, and corrected for relatedness and other batch effects using the SVA R package (see above for further information). Here, we focussed solely on the p value associated with the Student’s t-statistic for the interaction term S ijk E ijk . We applied a Bonferroni correction to the interaction p values for SNP-wise multiple testing within gene and retained the most significant putative eSNP–eGene pair from each gene. We then assessed this set of “best” eSNP–eGene p values for significance across all 15,632 genes at the false discovery rate (FDR) threshold of 0.05 by transforming the set into q values46). This represented the set of significant discovery eSNP–eGene pairs to be tested in the replication set. We then reported the environmental eSNP–eGene pairs that were significant (replicated) in the replication cohort (q value < 0.05, adjusted for the ten pairs being tested) and had the same direction of effect in both cohorts (n = 4 out of the 10).

To provide support for the replicated environmental eSNP–eGene pairs that we reported as significant, we estimated “honest” empirical p values for the whole sample (discovery + replication) using permutation: for each eSNP–eGene pair we performed the same eQTL modeling (Y ijk  ~ S ijk  + E ijk  + S ijk E ijk ) and permuted the expression values (Y) before obtaining the test statistic (Student’s t) for the interaction term. By repeating this procedure 1000 times for each eSNP–eGene pair, we built null distributions to assess the original observed (not permuted) t-statistics. The empirical permutation p value for each eSNP–eGene pair was taken as the proportion of permutation t-statistics larger than the observed t-statistic (Supplementary Fig. 16, Supplementary Table 9).

### Impact of lower frequency variants

We performed a CoIA analysis between all eSNPs of significant eGenes and endophenotypic traits. To do so, we resampled 1000 times, without replacement, 420 individuals from the cohort, and performed a CoIA at each step between endophenotypes showing variation across environments and the eSNPs. The median value for each endophenotype-eSNP correlation from the 1000 CoIA was calculated. The CoIA results are the correlations between eSNPs and the endophenotypic traits values. We then tested whether the strength of these correlations between eSNPs and endophenotypic traits were related to the MAF of the eSNP by examining the odds ratio of observing less common variants (MAF between 0.05 and 0.1, compared to common variants of MAF > 0.1) for stronger endophenotypic associations (Supplementary Fig. 18). The MAF was estimated from the complete cohort data.

### Data availability

Genotyping, expression, health phenotypes, and exposure data used in this study are available from CARTaGENE (www.cartagene.qc.ca) or the CPTP portal (http://portal.partnershipfortomorrow.ca) upon request. The built environment data set is publicly available from the Quebec government data portal. The air pollution data set is available upon request to Air Health Effects division, Government of Canada. All environmental data sources are detailed in Supplementary Table 4.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## References

1. 1.

Gibson, G. The environmental contribution to gene expression profiles. Nat. Rev. Genet. 9, 575–581 (2008).

2. 2.

Rappaport, S. M. & Smith, M. T. Epidemiology. Environment and disease risks. Science 330, 460–461 (2010).

3. 3.

Ye, C. J. et al. Intersection of population variation and autoimmunity genetics in human T cell activation. Science 345, 1254665 (2014).

4. 4.

Wu, S., Powers, S., Zhu, W. & Hannun, Y. A. Substantial contribution of extrinsic risk factors to cancer development. Nature 529, 43–47 (2016).

5. 5.

Nelson, M. R. et al. An abundance of rare functional variants in 202 drug target genes sequenced in 14,002 people. Science 337, 100–104 (2012).

6. 6.

Grubert, F. et al. Genetic control of chromatin states in humans involves local and distal chromosomal interactions. Cell 162, 1051–1065 (2015).

7. 7.

Carr, E. J. et al. The cellular composition of the human immune system is shaped by age and cohabitation. Nat. Immunol. 17, 461–468 (2016).

8. 8.

Kilpeläinen, T. O. et al. Physical activity attenuates the influence of FTO variants on obesity risk: a meta-analysis of 218,166 adults and 19,268 children. PLoS Med. 8, e1001116 (2011).

9. 9.

Franks, P. W., Pearson, E. & Florez, J. C. Gene-environment and gene-treatment interactions in type 2 diabetes: progress, pitfalls, and prospects. Diabetes Care 36, 1413–1421 (2013).

10. 10.

Aguirre-Gamboa, R. et al. Differential effects of environmental and genetic factors on T and B cell immune traits. Cell Rep. 17, 2474–2487 (2016).

11. 11.

Ter Horst, R. et al. Host and environmental factors influencing individual human cytokine responses. Cell 167, 1111–1124 (2016).

12. 12.

Gref, A. et al. Genome-wide interaction analysis of air pollution exposure and childhood asthma with functional follow-up. Am. J. Respir. Crit. Care Med. 195, 1373–1383 (2017).

13. 13.

Ward-Caviness, C. K. et al. Genetic variants in the bone morphogenic protein gene family modify the association between residential exposure to traffic and peripheral arterial disease. PLoS ONE 11, e0152670 (2016).

14. 14.

Ward-Caviness, C. K. et al. A genome-wide trans-ethnic interaction study links the PIGR-FCAMR locus to coronary atherosclerosis via interactions between genetic variants and residential exposure to traffic. PLoS ONE 12, e0173880 (2017).

15. 15.

Idaghdour, Y. & Awadalla, P. Exploiting gene expression variation to capture gene-environment interactions for disease. Front. Genet. 3, 228 (2012).

16. 16.

Marigorta, U. M. & Gibson, G. A simulation study of gene-by-environment interactions in GWAS implies ample hidden effects. Front. Genet. 5, 225 (2014).

17. 17.

Patel, C. J. & Ioannidis, J. P. A. Studying the elusive environment in large scale. JAMA 311, 2173–2174 (2014).

18. 18.

Idaghdour, Y. et al. Evidence for additive and interaction effects of host genotype and infection in malaria. Proc. Natl Acad. Sci. USA 109, 16786–16793 (2012).

19. 19.

Awadalla, P. et al. Cohort profile of the CARTaGENE study: Quebec’s population-based biobank for public health and personalized genomics. Int. J. Epidemiol. 42, 1285–1299 (2013).

20. 20.

Hussin, J. G. et al. Recombination affects accumulation of damaging and disease-associated mutations in human populations. Nat. Genet. 47, 400–404 (2015).

21. 21.

Hodgkinson, A. et al. High-resolution genomic analysis of human mitochondrial RNA sequence variation. Science 344, 413–415 (2014).

22. 22.

Labelle, R., Brand, A., Buteau, S. & Smargiassi, A. Hospitalizations for respiratory problems and exposure to industrial emissions in children. Environ. Pollut. 4, 77 (2015).

23. 23.

Doyon, B., Bélanger, D. & Gosselin, P. The potential impact of climate change on annual and seasonal mortality for three cities in Quebec, Canada. Int. J. Health Geogr. 7, 1 (2008).

24. 24.

Charbonneau, H., Desjardins, B., Légaré, J. & Denis, H. The population of the St-Lawrence Valley, 1608–1760. A Population History of North America 99–142 (Cambridge, Cambridge University Press, 2000).

25. 25.

Gauvreau, D., Guérin, M. & Hamel, M. De Charlevoix au Saguenay: mesure et caractéristiques du mouvement migratoire avant 1911. Histoire d’un génome: Population et génétique dans l’est du Québec145–159 (Presses de l’Université du Québec, Sillery, 1991).

26. 26.

Jette, R., Gauvreau, D. & Guérin, M. Aux origines d’une région: le peuplement fondateur de Charlevoix avant 1850. Histoire d’un génome: Population et génétique dans l’est du Québec 75–106 (Presses de l’Université du Québec, Sillery, 1991).

27. 27.

Casals, F. et al. Whole-exome sequencing reveals a rapid change in the frequency of rare functional variants in a founding population of humans. PLoS Genet. 9, e1003815 (2013).

28. 28.

Roy-Gagnon, M.-H. et al. Genomic and genealogical investigation of the French Canadian founder population structure. Hum. Genet. 129, 521–531 (2011).

29. 29.

Leek, J. T. & Storey, J. D. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genet. 3, 1724–1735 (2007).

30. 30.

Wood, J. Canadian Environmental Indicators - Air Quality (Fraser Institute, 2012), Vancouver, British Columbia, Canada.

31. 31.

Dray, S., Chessel, D. & Thioulouse, J. Co-inertia analysis amd the linking of ecological tables. Ecology 84, 3078–3089 (2003).

32. 32.

Kelly, F. J. & Fussell, J. C. Linking ambient particulate matter pollution effects with oxidative biology and immune responses. Ann. N. Y. Acad. Sci. 1340, 84–94 (2015).

33. 33.

Yang, W. & Omaye, S. T. Air pollutants, oxidative stress and human health. Mutat. Res. 674, 45–54 (2009).

34. 34.

Nourshargh, S. & Alon, R. Leukocyte migration into inflamed tissues. Immunity 41, 694–707 (2014).

35. 35.

Campen, M. J., Lund, A. & Rosenfeld, M. Mechanisms linking traffic-related air pollution and atherosclerosis. Curr. Opin. Pulm. Med. 18, 155–160 (2012).

36. 36.

Chuang, K.-J., Chan, C.-C., Su, T.-C., Lee, C.-T. & Tang, C.-S. The effect of urban air pollution on inflammation, oxidative stress, coagulation, and autonomic dysfunction in young adults. Am. J. Respir. Crit. Care Med. 176, 370–376 (2007).

37. 37.

Brook, R. D. et al. Particulate matter air pollution and cardiovascular disease: an update to the scientific statement from the American Heart Association. Circulation 121, 2331–2378 (2010).

38. 38.

Dominici, F. et al. Fine particulate air pollution and hospital admission for cardiovascular and respiratory diseases. JAMA 295, 1127–1134 (2006).

39. 39.

Zhou, Z. et al. Transcriptomic analyses of the biological effects of airborne PM2.5 exposure on human bronchial epithelial cells. PLoS ONE 10, e0138267 (2015).

40. 40.

Paolicchi, A. et al. Human atherosclerotic plaques contain gamma-glutamyl transpeptidase enzyme activity. Circulation 109, 1440–1440 (2004).

41. 41.

Lee, D.-H. & Jacobs, D. R. Jr. Is serum gamma-glutamyltransferase a marker of exposure to various environmental pollutants? Free Radic. Res. 43, 533–537 (2009).

42. 42.

Markevych, I. et al. Air pollution and liver enzymes. Epidemiology 24, 934–935 (2013).

43. 43.

Lee, D. S. et al. Gamma glutamyl transferase and metabolic syndrome, cardiovascular disease, and mortality risk: the Framingham Heart Study. Arterioscler. Thromb. Vasc. Biol. 27, 127–133 (2007).

44. 44.

Manolio, T. A. et al. Finding the missing heritability of complex diseases. Nature 461, 747–753 (2009).

45. 45.

West-Eberhard, M. J. Developmental Plasticity and Evolution, Oxford University Press, Oxford, UK (2003).

46. 46.

Storey, J. D. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann. Stat. 31, 2013–2035 (2003).

47. 47.

Arrowsmith, C. H., Bountra, C., Fish, P. V., Lee, K. & Schapira, M. Epigenetic protein families: a new frontier for drug discovery. Nat. Rev. Drug Discov. 11, 384–400 (2012).

48. 48.

Poncet-Montange, G. et al. Observed bromodomain flexibility reveals histone peptide- and small molecule ligand-compatible forms of ATAD2. Biochem. J. 466, 337–346 (2015).

49. 49.

Ciró, M. et al. ATAD2 is a novel cofactor for MYC, overexpressed and amplified in aggressive tumors. Cancer Res. 69, 8491–8498 (2009).

50. 50.

Wu, G. et al. Epigenetic high regulation of ATAD2 regulates the Hh pathway in human hepatocellular carcinoma. Int. J. Oncol. 45, 351–361 (2014).

51. 51.

Battle, A. et al. Characterizing the genetic basis of transcriptome diversity through RNA-sequencing of 922 individuals. Genome Res. 24, 14–24 (2014).

52. 52.

Zhao, J. et al. A burden of rare variants associated with extremes of gene expression in human peripheral blood. Am. J. Hum. Genet. 98, 299–309 (2016).

53. 53.

Gibson, G. & Wagner, G. Canalization in evolutionary genetics: a stabilizing theory? Bioessays 22, 372–380 (2000).

54. 54.

Paaby, A. B. & Rockman, M. V. Cryptic genetic variation: evolution’s hidden substrate. Nat. Rev. Genet. 15, 247–258 (2014).

55. 55.

10 Years of Data from the National Air Pollution Surveillance (NAPS), Analysis and Air Quality Section, Air Quality Research Division, Science and Technology Branch, Environment Canada. Goverment of Canada, Ottawa, Canada (2013).

56. 56.

Ardlie, K. G. et al. The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science 348, 648–660 (2015).

57. 57.

Dobin, A. et al STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).

58. 58.

Lappalainen, T. et al. Transcriptome and genome sequencing uncovers functional variation in humans. Nature 501, 506–511 (2013).

59. 59.

Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).

60. 60.

Anders, S., Pyl, P. T. & Huber, W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169 (2014).

61. 61.

Lawson, D. J., Hellenthal, G., Myers, S. & Falush, D. Inference of population structure using dense haplotype data. PLoS Genet. 8, e1002453 (2012).

62. 62.

Delaneau, O., Zagury, J.-F. & Marchini, J. Improved whole-chromosome phasing for disease and population genetic studies. Nat. Methods 10, 5–6 (2013).

63. 63.

Bherer, C. et al. Admixed ancestry and stratification of Quebec regional populations. Am. J. Phys. Anthropol. 144, 432–441 (2011).

64. 64.

Howie, B., Fuchsberger, C., Stephens, M., Marchini, J. & Abecasis, G. R. Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. Nat. Genet. 44, 955–959 (2012).

65. 65.

Li, S. et al. Detecting and correcting systematic variation in large-scale RNA sequencing data. Nat. Biotechnol. 32, 888–895 (2014).

66. 66.

Kukurba, K. R. et al. Impact of the X chromosome and sex on regulatory variation. Genome Res. 26, 768–777 (2016).

67. 67.

Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).

68. 68.

Hystad, P. et al. Creating national air pollution models for population exposure assessement in Canada. Environ. Health Perspect. 119, 1123–1129 (2011).

69. 69.

van Donkelaar, A. et al. Global estimates of ambient fine particulate matter concentrations from satellite-based aerosol optical depth: development and application. Environ. Health Perspect. 118, 847–855 (2010).

70. 70.

Boys, B. L. et al. Fifteen-year global time series of satellite-derived fine particulate matter. Environ. Sci. Technol. 48, 11109–11118 (2014).

71. 71.

Dray, S., Dufour, A. B. & Chessel, D. The ade4 package-II: two-table and K-table methods. R. News 7, 47–52 (2007).

72. 72.

Reimand, J., Arak, T. & Vilo, J. G: Profiler—a web server for functional interpretation of gene lists (2011 update). Nucleic Acids Res 39, W307–W315 (2011).

73. 73.

Hartigan, J. A. & Wong, M. A. Algorithm AS 136: a K-means clustering algorithm. J. R. Stat. Soc. Ser. C Appl. Stat. 28, 100–108 (1979).

## Acknowledgements

We thank the Awadalla Lab for comments on the manuscript, as well as Paul C. Boutros and Veronica Y. Sabelnykova from OICR. We acknowledge financial support from Fond de Recherche du Québec—Santé (FRQS), Genome Quebec, Fonds de Recherche du Québec—Nature et Technologies (FQRNT), Canadian Foundation of Innovation, Ontario Ministry of Research and Innovation Principal Investigator Award (P.A.) and a Canadian Institute of Health Research award (#EC3-144623) to P.A. M.-J.F. is a CIHR Neuroinflammation Postdoctoral Fellow. F.C.L. is a FRQS Postdoctoral Fellow. A.H. is a FRQS Postdoctoral Fellow and currently holds an eMedLab Career Development Fellowship as part of the Medical Bioinformatics Initiative funded by the Medical Research Council, UK (grant number MR/L016311/1). Requests for data published here should be submitted to access@cartagene.qc.ca citing this study.

## Author information

### Affiliations

1. #### Ontario Institute for Cancer Research, Toronto, ON, M5G 0A3, Canada

• Marie-Julie Favé
• , Fabien C. Lamaze
• , David Soave
• , Vanessa Bruat
• , Jean-Christophe Grenier
• , Elias Gbeha
2. #### Sainte-Justine Research Center, Faculty of Medicine, University of Montreal, Montreal, QC, H3T 1C5, Canada

• Marie-Julie Favé
• , Alan Hodgkinson
• , Héloïse Gauvin
• , Vanessa Bruat
• , Jean-Christophe Grenier
3. #### Department of Medical and Molecular Genetics, Guy’s Hospital, King’s College London, London, WC2R 2LS, UK

• Alan Hodgkinson

• Héloïse Gauvin
5. #### Department of Environmental Health and Occupational Health, University of Montreal, Montreal, QC, H3N 1X9, Canada

• Audrey Smargiassi
6. #### Health Canada, Air Health Science Division, Ottawa, ON, K1A 0K9, Canada

• Markey Johnson
7. #### NYU Abu Dhabi, Abu Dhabi, UAE

• Youssef Idaghdour

### Contributions

M.-J.F., Y.I., and P.A. conceptualized the study. E.G. and Y.I. performed the experimental procedures for sequencing and genotyping. M.-J.F., A.H., J.-C.G., M.J., A.S., and V.B. prepared the data and performed quality control. M.-J.F., F.C.L., D.S., V.B., J.-C.G., and H.G. performed bioinformatics and statistical analyses. M.-J.F., F.C.L., D.S., K.S., V.B., and P.A. wrote and revised the manuscript.

### Competing interests

The authors declare no competing financial interests.

## Electronic supplementary material

### DOI

https://doi.org/10.1038/s41467-018-03202-2

• 1.
• Angela Mo
• , Urko M. Marigorta
• , Dalia Arafat
• , Lai Hin Kimi Chan
• , Lori Ponder
• , Se Ryeong Jang
• , Jarod Prince
• , Subra Kugathasan


text:
# Lesson 11

Skills for Modeling with Mathematics

These materials, when encountered before Algebra 1, Unit 5, Lesson 11 support success in that lesson.

## 11.1: Which One Doesn’t Belong: Four Graphs (5 minutes)

### Warm-up

This warm-up prompts students to compare four graphs. It gives students a reason to use language precisely (MP6). It gives the teacher an opportunity to hear how students use terminology and talk about characteristics of the items in comparison to one another.

### Launch

Arrange students in groups of 2–4. Display the graphs for all to see. Give students 1 minute of quiet think time and then time to share their thinking with their small group. In their small groups, ask each student to share their reasoning why a particular item does not belong, and together find at least one reason each item doesn’t belong.

### Student Facing

Which one doesn’t belong?

### Activity Synthesis

Ask each group to share one reason why a particular item does not belong. Record and display the responses for all to see. After each response, ask the class if they agree or disagree. Since there is no single correct answer to the question of which one does not belong, attend to students’ explanations and ensure the reasons given are correct. During the discussion, ask students to explain the meaning of any terminology they use, such as increasing, decreasing, linear, nonlinear, or quadrant. Also, press students on unsubstantiated claims.

## 11.2: Which Model? (20 minutes)

### Activity

The purpose of this activity is to learn or remember a straightforward way to check if a relationship is linear or exponential (subtracting and dividing successive terms). Also, in the synthesis, there is an opportunity for students to rehearse creating a scatter plot, and recall how to recognize that a graph of a relationship might indicate that it is linear or exponential.

### Student Facing

Consider each situation:

1. A person starts with $24,000 in a savings account. Each month, she deposits an additional$2,000 in the account.
2. A 30-year old puts $24,000 in a retirement account that increases by 10% each year. 3. The value of a car depreciates by a factor of $$\frac45$$ of the car’s value every year. The car initially cost$24,000.
4. A farmer has stored 24,000 pounds of grain. His cows eat 4,800 pounds of grain per month.

Match each situation to one of these tables:

A.

$$x$$ $$y$$
0 24000
1 19200
2 15360
3 12288

B.

$$x$$ $$y$$
0 24000
1 26000
2 28000
3 30000

C.

$$x$$ $$y$$
0 24000
1 19200
2 14400
3 9600

D.

$$x$$ $$y$$
0 24000
1 26400
2 29040
3 31944

### Activity Synthesis

Questions for discussion:

• “How can you tell if the relationship is linear or exponential?” (See if there is a constant rate of change by subtracting successive terms and seeing if the difference is constant. See if there is a constant growth factor by dividing successive terms and seeing if the quotient is constant.)
• Look more closely at the table for the depreciating car. How could you compute that the growth factor is $$\frac45$$? (If you divide any two successive terms, like $$19,\!200 \div 2,\!400$$ or $$15,\!360 \div 19,\!200$$, you get 0.8 which is equal to $$\frac45$$.

Demonstrate how to create a scatter plot using technology using the table of values about the farmer’s grain. Then, ask students to use technology to create a scatter plot using the table of values about the car depreciation. Display these graphs side by side and observe ways in which the graphs suggest that one would be better modeled with a linear function and the other with an exponential function.

If students have access to the digital version of the materials, the “Graphing Calculator” tool under Math Tools is recommended. Add a new table, enter the values, and then adjust the graphing window (under the wrench button) so that the scatter plot is visible. Here is an example https://www.desmos.com/calculator/gdxthrlgcz

## 11.3: Growth of a Small Business (20 minutes)

### Activity

This activity is an opportunity to practice deciding between a linear and exponential model based on data given in a table. It is also an opportunity to practice creating a scatter plot using technology. The data is not perfectly linear or exponential, which is addressed in the launch.

### Launch

Tell students that in this activity, they have an opportunity to decide whether a linear or exponential model is better for a given set of data and create a scatter plot and a model. Draw students’ attention to the word “approximately” in the task statement. In this case, the data can not be fit perfectly with a model. So, when they are checking for constant differences or quotients, they’re looking for a number that best encapsulates the differences or quotients they find.

### Student Facing

Here are two sets of data representing the annual revenue of two different small businesses for the past ten years. One of them had growth that was approximately linear, and one of them had growth that was approximately exponential. The revenue is expressed in thousands of dollars.

 year revenue 0 1 2 3 4 5 6 7 8 9 61.2 68.4 74.9 83.1 88.5 96.4 104.1 109.9 117 125.2

 year revenue 0 1 2 3 4 5 6 7 8 9 40 47.9 57 70.1 82.4 99.5 118.9 144.1 172 205.8
1. Which company’s growth could be modeled by a linear function, and which by an exponential function?
2. For the company with exponential growth:
1. What was the growth factor?
2. What is an equation that represents the relationship between year and revenue?
3. Use technology to make a scatter plot of the data and also graph your equation. If your equation does not look like a good model for the data, adjust it until you have a good model.
3. For the company with linear growth:
1. What was the rate of change?
2. What is an equation that represents the relationship between year and revenue?
3. Use technology to make a scatter plot of the data and also graph your equation. If your equation does not look like a good model for the data, adjust it until you have a good model.

format: plain


text:
# Learning Scala as a Python Programmer: Functional Programming Features for Control Flow

## Recap

In my previous post on key principles of functional programming, I explained how the functional programming paradigm differs from imperative programming, and discussed how the the concepts of idempotency and avoidance of side effects are linked to the property of referential transparency that enables equational reasoning in functional programming.

Before we dive into some of the features of functional programming, let’s start with a personal anecdote during my first 3 months of writing Scala code.

## There is no “If-Else” in Functional Code

I was writing a pure Scala function for a custom Spark UDF which computes revenue adjustments based on a custom tiered adjustment expressed in JSON string. While attempting to express the business logic in pure functional code (since that is the team’s coding style), I got pretty frustrated with my perceived drop in productivity to the point whereby I introduced “if-else” logic into my code in a bid to “get the job done”.

Let’s just say that I learnt a pretty tough lesson during code review for that particular merge request.

“No if-else in functional code, this is not imperative programming… No ifs, no elses.

Without “if-else”, how do we write “control flow” in functional programming?

The long answer: A combination of function composition and functional data structures.

As a deep-dive on each functional design pattern can be pretty lengthy, the focus of this post is to provide an overview of function composition and how it enables a more intuitive approach to designing data pipelines.

## A Brief Intro to Function Composition

Function composition

In mathematics, function composition is an operation that takes two functions f and g in sequence and forms a composite function h such that h(x) = g(f(x)) - function g is applied to the result of applying the function f to a generic input x. Mathematically, this operation can be expressed as:

$f:X\rightarrow&space;Y,&space;g:Y\rightarrow&space;Z\Rightarrow&space;g&space;\circ&space;f:X\rightarrow&space;Z$

where $g&space;\circ&space;f$ is a composite function.

Intuitively, the composite function maps x in X to g(f(x)) in domain Z for all values in domain X.

A useful analogy to illustrate the concept of function composition is making butter toast in an oven with a slice of bread and cold butter. There are two possible operations:

1. Toasting in the oven (operation f)
2. Spreading butter over the widest surface (operation g)

If we toast the bread in the oven first and spread cold butter over the widest surface of what comes out of the oven, we get a slice of toasted bread with cold butter spread ($g&space;\circ&space;f$).

If we spread cold butter over the widest surface of the bread first and toast the bread with cold butter spread in the oven, we get a slice of toasted bread with warm butter spread ($f&space;\circ&space;g$). And we know that “cold butter spread” != “warm butter spread”.

From these examples, we can intuitively infer that the order of function application matters in function composition. ($g&space;\circ&space;f&space;\neq&space;f&space;\circ&space;g$)

Similarly in designing data pipelines, we often write data transformations by applying functions to results of other functions. The ability to compose functions encourages refactoring of repeated code segments into functions for maintainability and reusability.

## Functions as First-Class Objects

The core idea in functional programming is: functions are values.

This feature implies that a function can be [2,3]:

1. assigned to a variable
2. passed as a parameter to other functions
3. returned as a value from other functions

For this to work, functions must be first-class objects (and stored in data structures) in the runtime environment - just like numbers, strings and arrays. First-class functions are supported in all functional languages including Scala, as well as some interpreted languages such as Python.

## Higher-Order Functions

A key implication resulting from the concept of functions as first-class objects is that function composition can be naturally expressed as a higher-order function.

A higher-order function has at least one of the following properties:

1. Accepts functions as parameters
2. Returns a function as a value

An example of a higher-order function is map.

When we look at the documentation for the Python built-in function map, it is stated that the map function takes in another function and an iterable as input parameters and returns an iterator that yields the results [4].

In Scala, each of the collection classes in package scala.collections and its subsets contain the map method that is defined by the following function signatures on ScalaDoc [5]:

def map[B](f: (A) => B): Iterable[B]    // for collection classes
def map[B](f: (A) => B): Iterator[B]    // for iterators that access elements of a collection

What the function signatures mean is that map takes a function input parameter f, and f transforms a generic input of type A to a resulting value of type B.

To square each value in a collection of integers, the iterative approach is to traverse each element in the collection, square the element and append the result to a collection of results that expands in length with each iteration.

• In Python:

  def square(x):
return x * x

def main(args):

collection = [1,2,3,4,5]
# initialize list to hold results
squared_collection = []
# loop till the end of the collection
for num in collection:
# square the current number
squared = square(num)
# add the result to list
squared_collection.append(squared)

print(squared_collection)

In the iterative approach, two state changes occur at each iteration within the loop:

1. The squared variable holding the result returned from the square function; and
2. The collection holding the results of the square function.

To perform the same operation using a functional approach (i.e. without using mutable variables), the map function can be used to “map” each element in the collection to a new collection with the same number of elements as the input collection - by applying the square operation to each element and collecting the results into the new collection.

• In Python:

  def square(x):
return x * x

def main(args):

collection = [1,2,3,4,5]
squared = list(map(square, collection))
print(squared)

• In Scala:

  object MapSquare {

def square(x: Int): Int = {
x * x
}

def main(args: Array[String]) {

val collection = List[1,2,3,4,5]
val squared = collection.map(square)
println(squared)
}
}

In both implementations, the map function accepts an input function that is applied to each element in a collection of values and returns a new collection containing the results. As map has the property of accepting another function as a parameter, it is a higher-order function.

A few quick side-notes on differences between Python and Scala implementations:

• Python map vs Scala map: An iterable function such as list is needed to convert the iterator returned from the Python map function into an iterable. In Scala, there is no need for explicit conversion of the result from the map function to an iterable, as all methods in the Iterable trait are defined in terms of an abstract method, iterator, which returns an instance of the Iterator trait that yields the collection’s elements one by one [6].
• How values are returned from a function: While the return keyword is used in Python to return a function result, the return keyword is rarely used in Scala. Instead, the last line within a function declaration is evaluated and the resultant value is returned when defining a function in Scala. In fact, using the return keyword in Scala is not good practice for functional programming as it abandons the current computation and is not referentially transparent [7-8].

### Anonymous Functions

When using higher-order functions, it is often convenient to be able to call input function parameters with function literals or anonymous functions without having to define them as named function objects before they can be used within the higher-order function.

In Python, anonymous functions are also known as lambda expressions due to their roots in lambda calculus. An anonymous function is created with the lambda keyword and wraps a single expression without using def or return keywords. For example, the square function in the previous example in Python can be expressed as an anonymous function in the map function, where the lambda expression lambda x: x * x is used as a function input parameter to map:

def main(args):

collection = [1,2,3,4,5]
squared = map(lambda x: x * x, collection)
print(squared)

In Scala, an anonymous function is defined in-line with the => notation - where the function arguments are defined to the left of the => arrow and the function expression is defined to the right of the => arrow. For example, the square function in the previous example in Scala can be expressed as an anonymous function with the (x: Int) => x * x syntax and used as a function input parameter to map:

object MapSquareAnonymous {

def main(args: Array[String]) {
val collection = List[1,2,3,4,5]
val squared = collection.map((x: Int) => x * x)
println(squared)
}
}

A key benefit of using anonymous functions in higher-order functions is that single-use single-expression functions need not be wrapped explicitly within a named function definition, hence optimizing lines of code and improving code maintainability.

## Recursion as a form of “functional iteration”

Recursion is a form of self-referential function composition - a recursive function takes the results of (smaller instances of) itself and uses them as inputs to another instance of itself. To prevent an infinite loop of recursive calls, a base case is required as a terminating condition to return a result without using recursion.

A classic example of recursion is the factorial function, which is defined as the product of all positive integers less than or equal to an integer n:

$n!=n\cdot&space;(n-1)\cdot&space;(n-2)\cdots&space;3\cdot&space;2\cdot&space;1$

There are two possible iterative approaches to implementing a factorial function: using for loop, and using while loop.

• In Python:

  def factorial_for(n):
# initialize variable to hold factorial
fact = 1
# loop from n to 1 in decrements of 1
for num in range(n, 1, -1):
# multiply current number with the current product
fact = fact * num
return fact

def factorial_while(n):
# initialize variable to hold factorial
fact = 1
# loop till n reaches 1
while n >= 1:
# multiply current number with the current product
fact = fact * n
# subtract the number by 1
n = n - 1
return fact

In both iterative implementations of the factorial function, two state changes occur at each iteration within the loop:

1. The factorial variable storing the current product; and
2. The number being multiplied.

To implement the factorial function using a functional approach, recursion is useful in dividing the problem into subproblems of the same type - in this case, the product of n and *(n-1)!).

The basic recursive approach for the factorial function looks like this:

• In Python:
  def factorial(n):
# base case to return value
if n <= 0: return 1
# recursive function call with another set of inputs
return n * factorial(n-1)

• In Scala:
  def factorial(n: Int): Long = {
if (n <= 0) 1 else n * factorial(n-1)
}

For the basic recursive approach, the factorial of 5 is evaluated in the following manner:

factorial(5)
if (5 <= 0) 1 else 5 * factorial(5 - 1)
5 * factorial(4)    // factorial(5) is added to call stack
5 * (4 * factorial(3))  // factorial(4) is added to call stack
5 * (4 * (3 * factorial(2)))    // factorial(3) is added to call stack
5 * (4 * (3 * (2 * factorial(1))))  // factorial(2) is added to call stack
5 * (4 * (3 * (2 * (1 * factorial(0)))))    // factorial(1) is added to call stack
5 * (4 * (3 * (2 * (1 * 1)))) // factorial(0) returns 1 to factorial(1)
5 * (4 * (3 * (2 * 1))) // factorial(1) return 1 * factorial(0) = 1 to factorial(2)
5 * (4 * (3 * 2))   // factorial(2) return 2 * factorial(1) = 2 to factorial(3)
5 * (4 * 6) // factorial(3) return 3 * factorial(2) = 6 to factorial(4)
5 * 24  // factorial(4) returns 4 * factorial(3) = 24 to factorial(5)
120 // factorial(5) returns 5 * factorial(4) = 120 to global execution context

For n = 5, the evaluation of the factorial function involves 6 recursive calls to the factorial function including the base case.

While the basic recursive approach expresses the factorial function more closely with its definition (and more naturally) compared with the iterative approach, it also uses more memory as each function call is pushed to the call stack as a stack frame and popped off the call stack when the function call returns a value.

For larger values of n, the recursion gets deeper with more function calls to itself and more space has to be allocated to the call stack. When the space needed to store the function calls exceeds the capacity for the call stack, a stack overflow occurs!

### Tail Recursion and Tail-call Optimization

To prevent infinite recursion from causing stack overflow and crashing the program, some optimizations have to be made to the recursive function in order to reduce consumption of stack frames in the call stack. A possible approach in optimizing the recursive function is by rewriting it as a tail recursive function.

A tail recursive function calls itself recursively and does not perform any computation after the recursive call returns. A function call is a tail call when it does nothing other than returning the value of the function call.

In functional programming languages such as Scala, tail-call optimization is typically included in the compiler to identify tail calls and compile the recursion to iterative loops that do not consume stack frames for each iteration. In fact, the stack frame can be reused for both the recursion function and the function being called within the recursion function [1].

With this optimization, the space performance for the recursion function can be reduced from O(N) to O(1) - from one stack frame per call to one stack frame for all calls [8]. In a way, a tail recursive function is a form of “functional iteration” with comparative performance to a loop.

For example, the factorial function can be expressed in the form of a tail recursion in Scala:

def factorialTailRec(n: Int): Long = {
def fact(n: Int, product: Long): Long = {
if (n <= 0) product
else fact(n-1, n * product)
}

fact(n, 1)
}

While tail-call optimization is automatically performed during compilation in Scala, it is not the case for Python. Moreover, there is a recursion limit in Python (the default value is 1000) as a prevention measure against an overflow of the C call stack for the CPython implementation.

## What’s next: Higher-Order Functions

In this post, we learn about:

1. Function composition
2. Higher-Order Functions as a key implication of functional programming
3. Recursion as a form of “functional iteration”

Have we found a replacement for “if-else” yet? Not entirely, but we now know how to write “loops” in functional programming using Higher-Order Functions and tail recursion.

In the next post, I will explore more on Higher-Order Functions and how they can be used in designing functional data pipelines.

## References

1. Functional Programming in Scala by Paul Chiusano and Rúnar Bjarnason
2. Functional Programming Simplified by Alvin Alexander
3. Functional Python Programming by Steven F. Lott, 2nd Edition
4. Built-in Functions - Python 3.9.6 Documentation
5. Scala Standard Library 2.13.6 - scala.collections.Iterable
6.  [Trait Iterable Collections Scala Documentation](https://docs.scala-lang.org/overviews/collections-2.13/trait-iterable.html)
7. tpolecat - The Point of No Return
8. Don’t Use Return in Scala? - Question - Scala Users
9. Tail Recursion - Functional Programming in OCaml
Written on July 4, 2021