dataesr/scientific-paragraphs-categorization huggingface dataset:
Given the practical difficulties of obtaining multiomics datasets from nonaccesible tissues, a question that remains unanswered is how deep molecular phenotyping in accessible tissues such as blood may aid in understanding the genetics of complex diseases. Here, we use data from the DIabetes REsearCh on patient straTification (DIRECT) consortium to investigate the genetic regulation of multiple molecular phenotypes. Using blood and plasma samples from DIRECT participants, we performed local (cis) and distal (trans) quantitative trait loci (QTLs) analysis in derived gene expression, protein and metabolite phenotypes. A deep analysis of these QTLs investigated allelic heterogeneity and pleotropic effects shared across molecular phenotypes with local and distal effects. Comparison with genetic associations across multiple tissues allowed us to evaluate the value of deep phenotyping in an accessible tissue to understand the genetics of complex traits.
Identification of local genetic regulation of molecular phenotypes The DIRECT dataset consists of 3029 individuals of European descent, with genotype information, gene expression quantified by RNA sequencing (RNA-seq), targeted proteomics (multiplexed immunoassays) and metabolites (targeted and untargeted mass spectrometry) in whole blood (Fig. 1A , Methods, Supplementary Table  1 ). To identify independent local (cis-) genetic effects that regulate mRNA expression (cis-eQTLs) and levels of circulating protein (cis-pQTLs), we performed QTL analyses followed by a stepwise regression conditional analysis for each phenotype. We identified independent significant associations for 94. 4% (15,305) of genes and 97. 3% (363) of proteins (Supplementary Table  2 , Methods, Supplementary Data 1-3), finding similarities in the regulation of both gene expression and proteins (Fig.
1B-E,  and  Fig. 1  | Multiomic QTL analysis identifies extensive allelic heterogeneity and pleotropic effects across molecular phenotypes. A The DIRECT consortium derived genetic, transcriptomic, proteomic and metabolite data from blood and plasma samples from 3029 individuals. Significant genetic associations (FDR < 0. 05) after linear regression between the molecular phenotypes and SNPs in cis (cis-QTLs) and trans (trans-QTLs or GWAS) were used to build a network of genetic perturbations affecting molecular phenotypes. Partially created with BioRender. com. B Location of the lead eSNP with respect to the TSS of the significantly associated genes (FDR < 0. 05) for cis-eQTL. The most associated eSNP per gene (primary) is shown in black (n = 15,305). Secondary cis-eQTL are shown in orange (n = 44,667).
A Using n = 3027 biologically independent samples, we detected a cis-pQTL for CCL16 in whole blood (P value = 9. 5e-243, n = 3029). The GTEx consortium reported a cis-eQTL, with the same SNP (rs10445391) affecting the expression of the gene in liver (n = 208). Violin plots show the median and first and last quartiles as defined by ggplot geom_violin function. Partially created with BioRender. com B Between 91. 2% (pancreatic islets) and 71. 6% (esophagus mucosa) of cis-eQTLs discovered by GTEx v8 were also active in whole blood DIRECT datasets (n = 3029) as shown by the π 1 values (y-axis). The number of P values per tissue used to calculate the π 1 estimates ranged from 334 in kidney to 14,920 in thyroid. C Comparison of the effect size of cis-eQTLs from pancreatic islets (InsPIRE) and whole blood (DIRECT).
336) and brain putamen basal ganglia (π1 = 0. 136) among others, while 15 tissues did not identify any significant replicated trans-eQTLs from DIRECT blood (π1 = 0, Supplementary Data 11). After multiple testing correction, we replicated 278 significant blood trans-eQTLs (237 unique genes, 159 unique SNPs), corresponding to 643 gene-SNP pairs across all tissues. In contrast, only 4 trans-pQTL were observed as trans-eQTLs in GTEx tissues after multiple testing correction: FCRL5-rs569841457 (adipose subcutaneous), MMP9-rs919377 adrenal gland, SPINK5-rs12462111 (nerve tibial) and OLR1-rs76604815 (uterus). Additionally, we evaluated the number of significant blood trans-eQTLs, trans-pQTLs and metabo-QTLs that replicated in other blood and plasma datasets. Using the eQTLGen dataset  26  , we were able to evaluate 514 gene-SNP pairs from DIRECT trans-eQTLs, of which 463 were also significant (Supplementary Data 11).
However, when we attempted to identify cis-eQTLs from other tissues, we could detect between 71. 6% and 91. 2% of these signals in our whole blood dataset as active in genetic regulation in other tissues. This suggests that a sufficiently large sample size in an accessible tissue can be informative of regulation in other tissues. However, we also reported up to a 20% decrease in the degree of sharing genetic signal across tissues when the sample size of the GTEx dataset increased with the transition from version 6 to version 8. Since this decrease was not uniform across tissues, researchers should evaluate carefully what accessible tissue may be most suitable for their research. The choice of tissue for the study of genetic regulation is even more critical for those aiming to understand the underlying Fig.
Cohort The DIRECT (Diabetes Research on Patient Stratification) consortium (Supplementary Data 19) includes pre-diabetic participants (target sample size 2,200-2,700) and patients with newly diagnosed type 2 diabetes (target sample size ~1,000) with detailed metabolic phenotyping. Characteristic of the cohort as well as inclusion/exclusion criteria have been described elsewhere  49. In short, fasting blood samples from venous blood were collected and DNA extractions and other biochemical analyses were carried out as detailed bellow. All measurements for each molecular data type were taken from different samples. Ethics approval for the study protocol was obtained by all the regional research ethics review boards (Lund, Sweden: 20130312105459927; Copenhagen, Denmark: H-1-2012-166 and H-1-2012-100; Amsterdam, Netherlands: NL40099.
Genotyping Genotyping was conducted using the Illumina HumanCore array (HCE24 v1. 0) and genotypes were called using Illumina's GenCall algorithm. Samples were excluded for any of the following reasons: call rate <97%; low or excess mean heterozygosity; gender discordance; and monozygosity. Genotyping quality control was then performed to provide high-quality genotype data for downstream analyses using the following criteria: call rate <99%; deviation from Hardy-Weinberg equilibrium (exact p < 0. 001); variants not mapped to human genome build GRCh37; and variants with duplicate chromosome positions. To identify possible ethnic outliers in the DIRECT data, we performed a principal component analysis (PCA) using the genotype data from our studied population (3102 samples; 547,644 markers) using the following cut-offs MAF > 0.
RNAseq data generation mRNA samples were processed and quality was assessed using the TapeStation Software (A. 01. 04) with an RNA Screen Tape from Agilent. Truseq Stranded mRNA from Illumina was used to generate libraries and their quality was evaluated using Qubit and TapeStation using DNA1000 Screen Tape. The samples were then sequenced on the Illumina HiSeq2000 platform using 49 bp paired-end reads. The 49-bp sequenced paired-end reads were mapped to the GRCh37 reference genome  52  with GEMTools 1. 7. 1  53. Exon quantifications were calculated for all elements annotated on GENCODE v19  54. Gene quantifications were calculated as FPKM values. This pipeline is fully described in Delaneau et al. 55  , as part of QTLtools. Samples with a total number of exonic reads lower than 5e + 06 reads or with a proportion of exonic reads over the total number of reads lower than 20% were considered of low quality and removed.
Proteins data Plasma proteins were measured using the Olink® Cardiometabolic, Cardiovascular II, Cardiovascular III, Development and Metabolism panels (Olink Proteomics AB, Uppsala, Sweden) according to the manufacturer's instructions. The Proximity Extension Assay (PEA) technology used for the Olink protocol has been well described  59. The obtained data was processed using Olink's NPX manager software version 0. 0. 85. 0. Internal and external controls were used for quality control and normalization of the data. Quality control included calculating the standard deviation for the detection control and the incubation/immuno controls and comparison of the results for the detection control and one of the incubation controls against the run median. All protein measurements were reported, but proteins with more than 50% samples below LOD levels were excluded, leaving a total of 373 proteins for further analysis from 3027 individuals.
Metabolites data Metabolites abundance was assessed using targeted and untargeted technologies. The assays, quality control and analyses were performed separately, and results combined for discussion. Plasma targeted metabolites were evaluated for 163 metabolites using a FIA-ESI-MS/MSbased targeted metabolomics approach with the AbsoluteIDQ TM p150 kit (BIOCRATES Life Sciences AG, Innsbruck, Austria) as described in Biocrates manual AS-P150. Mass spectrometric analyses were done on an API 4000 triple quadrupole system (Sciex Deutschland GmbH, Darmstadt, Germany) equipped with a 1200 Series HPLC (Agilent Technologies Deutschland GmbH, Böblingen, Germany) and a HTC PAL auto sampler (CTC Analytics, Zwingen, Switzerland) controlled by the software Analyst 1. 6. 2. Data evaluation for quantification of metabolite concentrations and quality assessment was performed with the software MultiQuant 3.
Controls included a pooled matrix sample generated by taking a small volume of each experimental sample serving as a technical replicate throughout the data set. Experimental samples were randomized across the platform run with QC samples spaced evenly among the injections. Metabolite concentration was assessed using Liquid Chromatography-Tandem Mass Spectrometry (LC-MS/MS). The LC-MS portion of the platform was based on a Waters ACQUITY ultra-performance liquid chromatography (UPLC) and a Thermo-Finnigan LTQ mass spectrometer operated at nominal mass resolution, which consisted of an electrospray ionization (ESI) source and linear ion-trap (LIT) mass analyzer. Raw data was extracted, peak-identified and QC processed using Metabolon's hardware and software. For studies spanning multiple days, a data normalization step was performed to correct variation resulting from instrument inter-day tuning differences.
Trans-QTL discovery Trans-QTL analysis was performed between molecular phenotypes with a genomic location (RPKM expression and proteins) and all SNPs further than 5MB from the TSS of the expressed gene or coding gene for the protein. For all associations, since phenotypes PCs may capture global trans effects removing true signals, we used residuals after removing known technical covariates with a linear mixed model implemented in the lme4 package in R 62. To control for false positives in trans-eQTLs  57  we removed: any gene with low mappability and SNPs located in repetitive regions and with mapping issues. We then performed a genome wide scan of SNP-gene associations using QTLtools  55  storing all P values < 1e-04. Pairs of SNP-genes with known crossmapability issues in a 100 kb window were then removed.
Multiple testing corrections for trans-eQTLs was done using 50 permutations. The best P values per phenotype were used to define a phenotype level threshold and calculate adjusted P values for each gene. To adjust for multiple testing across genes, we used the qvalue package to estimate the false discovery rate. For proteins, 100 permutations per protein were used. For the identification of secondary trans-QTL, we used a stepwise regression analysis after defining non-overlapping cis-windows around 1MB up-and down-stream of any significant trans-QTLs signal. This created 2021 genomic regions for expression and 2589 regions for proteins to be tested for any nearby signal. Then, we ran the standard conditional analysis using FastQTL to look for further signals.
Transcription factor enrichment To estimate the proportion of significant trans-SNPs which also affected local molecular phenotypes, we extracted all the cis-QTL P values for SNPs with a significant trans-SNP effect and calculated π 1 estimates using all significant associations and the qvalue package. The product of all P values per SNPs were used to calculate the probability of one gene to be associated to the trans-SNP, allowing us to calculate one P value per phenotype using π 1 enrichment. To investigate if cis-SNPs that also have an effect in trans were more likely to regulate a transcription factor (TF) in cis, we extracted all genes or proteins reported by Lambert et al. 22  as transcription factors (TFs). We then used a Fisher exact test to calculate the OR of a cis-eQTL affecting TFs, given the number of TFs in the whole dataset that have trans associations.
Tissue and phenotype specific associations, replications To identify the proportion of cis-eQTLs, cis-pQTLs and metabo-QTLs active in other tissues, we extracted the same SNP-gene in GTEx v6p  25  , GTEx v8  12  ,and InsPIRE 3 that were significantly associated in DIRECT. To identify the proportion of eQTLs that were identified in DIRECT as QTLs of any type, we extracted the SNP-gene pairs from the significant GTEx and Inspire significant associations. P values were then evaluated for enrichment using qvalue. For metabolites-QTLs with no TSS of reference, we used all cis-eQTLs in 45 studies and cis-pQTLs in DIRECT associated to the lead metabo-SNPs to calculate π 1 estimates. To identify the proportion of trans-eQTLs and trans-pQTLs active as trans-eQTLs in other tissues we extracted all gene-SNP and protein-SNPs pairs from the trans-QTL significant results in DIRECT-blood in GTEx v8 data and performed a linear regression.
We controlled for covariates reported by GTEx (PCs from genotypes 1 to 5, 15 to 60 inferred covariates depending on tissue, pcr, platform, and sex). We considered SNP-genes or SNP-proteins to be active in other tissues when a test for association in the GTEx v8 dataset was significant after multiple testing corrections (qvalue < 0. 05). Replication of distal QTLs was done by extracting P values for genetic association from publicly available summary statistics. For trans-eQTLs we used associations from the eQTLGen consortium (https://eqtlgen. org/)  26  , extracting all gene-SNP pairs. For trans-pQTLs we used associations from the INTERVAL consortium (http://www. phpc. cam. ac. uk/ceu/proteins/)  17. In both cases we evaluated the replication level using π1 and the qvalue package across all pairs commonly present between studies.
3 ): 1) Direction 1 model ( [1] and [2]): The genetic variant affects first phenotype 1, then phenotype 2; 2) Direction 2 model ( [3] and [4]): The genetic variant affects first phenotype 2, then phenotype 1; 3) Independent model ( [5] and [6]): The genetic variant affects both phenotypes 1 and 2 independently. The probabilities can be described in the following formulas, for one or two SNPs: 1. Direction We used normalized quantifications after removing any technical covariates. In addition, we removed from the phenotypes the effects of other QTLs by including SNPs as covariates, creating two phenotypes per trio. Computation to generate these pseudo-phenotypes were done using the -single-signal option from CaVEMaN  63. The phenotypes were used to compute the likelihood of the three possible BN topologies using the R/bnlearn package  64.
As cis-eQTL variants selected for theses analyses were previously found to be associated with the traits studied, we reported the probability that both GWAS and cis-eQTLs were shared as P(H4') = P(H4)/(P(H3) + P(H4)). H3 is the probability that both traits have different causal variants and H4 the probability of both traits sharing the same causal variant. For GWAS studies reporting different sample sizes per significant variant due to population differences, we calculated probabilities using the median number of samples across studies. To identify enrichment of GWAS associations among QTLs, we calculated and evaluated odd ratios using a fisher test. Since the GWAS catalog includes results from pQTL and metabolite-QTL studies, we removed any GWAS-SNP from blood proteins levels, metabolites levels and related phenotypes (Supplementary Data 18).
Figures Data and results figures were generated using R or Cytoscape 66 v3. 7. 2 for the networks. Figures  1A  and 3A  were partially created with BioRender. com. Figure  1B  was created using a network from Cytoscape and a lolliplot figure generated using the lolliplot() function from the trackViewer package  67. Figures were combined in panels using Inkscape v1. 0. 1. Tables and objects to load networks into Cytoscape are included in the Zenodo submission (https://zenodo. org/record/ 7521410).
Acknowledgements We are grateful to IMI-DIRECT study participants who volunteered for phenotyping, and clinical and technical staff across involved European study centres who contributed to recruitment and clinical assessment of study participants. The work leading to this publication has received support from the  Innovative Medicines Initiative Joint Undertaking  under grant agreement  115317  (DIRECT), resources of which are composed of financial contribution from the  European Union 's  Seventh Framework Programme  ( FP7/2007-2013 ) and  EFPIA  companies' in kind contribution. Information on the initiatives and activities of the IMI-DIRECT research consortium is available at https://directdiabetes. org.
Data availability The molecular and clinical raw data as well as the processed are available under restricted access due to the informed consent given by study participants, the various national ethical approvals for the present study, and the European General Data Protection Regulation (GDPR), individual-level clinical and molecular data cannot be transferred from the centralized IMI-DIRECT repository. Requests for access will be informed on how data can be accessed via the DIRECT secure analysis platform following submission of an appropriate application. The IMI-DIRECT data access policy is available at https:// directdiabetes. org. As described in the methods section we used the human genome build GRCh37 as a reference for genomic location of genotypes and transcriptomics data, and Gencode v19  54  for gene models and TSS information (https://www.
gencodegenes. org/ human/release_19. html). For functional enrichment analyses we used dataset from the Ensembl Variant Effect Predictor 14 (VEP) information v98 (https://grch37. ensembl. org/info/docs/tools/vep/ script/vep_download. html) and ChromHMM 13 models (https://egg2. wustl. edu/roadmap/data/byFileType/chromhmmSegmentations/ ChmmModels/imputed12marks/jointModel/final/). GTEx v6p and v8 summary statistics were accessed using the GTEx Portal (https:// www. gtexportal. org/home/). The 16 GWAS studies evaluated for colocalization and the links to their summary statistics are listed in Supplementary Data 15. Complete summary statistics including cis and trans genetic associations for gene expression, proteins and metabolites, as well as files to visualize networks on Cytoscape are freely available in the following link https://zenodo.
Instead of separate treatment of the vegetation canopy and ground, ISBA considers the composite soil-vegetation energy budget (Fig. 2B, C ). In the most detailed soil scheme ISBA-Diffusion (ISBA-DIF,  Boone et al. , 2000; Decharme et al. , 2011) , parameters may be defined by the user, or obtained from global or regional databases (e. g. Faroux et al. , 2013)  and pedotransfer functions  (Noilhan and Mahfouf, 1996; Peters-Lidard et al. , 1998). In the presence of full snow cover, the surface energy budget is solved by Crocus (Sect. 2. 2. 3). For partial snow cover, Crocus is used to solve the snow covered fraction while the energy balance of the snow-free fraction is computed by ISBA, and total surface energy fluxes are computed as weighted averages of the snow and snow-free fractions (Sect.
Model forcing Meteorological forcing consist of hourly observations of T a , V a , precipitation rate (P ), q a (T s ), R g and R A and atmospheric pressure. The available meteorological observations from the nearest meteorological stations were obtained from Finnish Me-   (Alekseychik et al. , 2022a; Mammarella et al. , 2019). At S-WET and S-FOR the shortwave and longwave radiation were obtained from the SMEAR database, while at N-WET and N-FOR data from FMI stations were used. The diffuse to total shortwave radiation ratio, r, was estimated as a function of the cosine of the sun zenith angle, µ. More specifically, a 3rd degree polynomial fit between r and µ was obtained using the atmospheric model SBDART  (Ricchiazzi et al. , 1998)  to simulate diffuse and total solar radiation in clear sky conditions.
Model evaluation data We use surface energy flux observations, height of snow (HS) and soil temperatures in model evaluation. The availability period of each variable is given in Appendix A (Table  A1 ). On all sites, reflected shortwave radiation (SWU) and outgoing longwave radiation (LWU) were measured using pyranometers and pyrgeometers, while ground heat flux (G) was measured using soil heat flux plates between 5 and 10 cm depths. The H and LE were measured by the eddy-covariance (EC) technique. The EC systems consist of USA-1 (METEK) threeaxis and Gill HS-50 sonic anemometers as well as closed-path LI-7000 and LI-7200 (Li-cor, Inc. ) CO 2 /H 2 O analysers  (Aurela et al. , 2015; Mammarella et al. , 2016; Alekseychik et al. , 2022b). The detailed procedure for obtaining the turbulent heat fluxes from raw eddy covariance data is detailed in  Aurela et al.
(2015) ;  Mammarella et al. (2016). In short, the sensible and latent heat fluxes were screened for instrument failure and data outliers, and data quality flags were made according to friction velocity (u * ) and flux stationarity (FST) criteria  (Foken et al. , 2005) : flag 2: all data (after screening of instrument failures and outliers) ). On N-WET the automated HS measurement is at 0. 7 m height, and therefore the exceeding snow depths were taken from biweekly manual measurements. To account for the spatial variability of snow depth in the forests, manual HS measurements from a snow course in the close proximity of the automated measurements were used  (Aalto et al. , 2022; Marttila et al. , 2021). Each site has different configuration of soil temperature sensors. At N-FOR and N-WET stations, soil temperatures are measured at 5, and 20 cm depths  (Aurela et al.
Insights on energy flux partitioning in boreal environments We used a novel dataset including all surface energy balance components from two peatland and two forest sites. The potential of such data has not yet been fully leveraged in LSM snow model evaluation. Observations showed that in winter the latent heat flux was minimal at all sites and the negative net radiation (down to -50 Wm -2 ) almost completely counterbalanced by the sensible heat flux (see Fig  13  and Fig. 4 ). The G had only small contribution to the winter energy budget on the studied peatlands, whereas it has been reported to have a rather important role for the open sites in Canadian Arctic  (Lackner et al. , 2021) , Siberia  (Langer et al. , 2011)  and Svalbard  (Langer et al. , 2011). This is due to the high heat capacity of the peatland, and its large water storage which progressively freezes from the top keeping the temperatures in the soil-snow surface nearly constant at minimum 0°C.