dataesr/scientific-paragraphs-categorization huggingface dataset:
Hence, RPS31 may be involved in chemical stress resistance by regulating protein metabolism. Avt1 and Spe4 phosphorylations were correlated with multiple morphological traits. Avt1 is involved in the import of large neutral amino acids into the vacuole  (Russnak et al, 2001). The phosphoresiduals of Avt1 were correlated with traits related to general cell shape and size (Dataset EV5-S9), and all of those traits were targets of the same QTL within hotspot chrVIII:1 (GPA1/STE20 hotspot, discussed below). Spe4 phospho-residuals were correlated with traits related to actin (Dataset EV5-S9). Spe4 is required for the synthesis of spermine, and disruption of spermine synthesis was previously associated with abnormal distribution of actin  (Balasundaram et al, 1991). Formal mediation analysis suggested that the phosphoresiduals of Avt1 and Spe4 significantly mediate the effects of genetic variation on the respective cell morphology traits (all FDR < 0.
phQTL effects on regulatory network states The above findings suggest that protein phosphorylation provides a crucial link between genomic variability and cellular (stress resistance) traits. Since these effects were likely mediated via signaling processes, we next investigated QTL effects on the states of signaling networks. To study QTL effects on the regulome, we integrated the five types of molecular traits with regulatory network information. We found that phResQTL targets were often functionally related to the putative causal genes underlying a given hotspot and/ or targets of kinases that were affected by the same locus. The IRA2 hotspot (XV:1)  (Smith & Kruglyak, 2008)  is an example for such a case. This hotspot affected resistance to six chemical compounds and had targets among all five types of molecular traits (Dataset EV5-S2, Fig  3B ).
polymorphisms of this cross in the IRA2 coding sequence affect the activity of the Ras/Pka pathway  (Smith & Kruglyak, 2008). Several targets of phQTLs in the IRA2 hotspot were functionally connected to IRA2 including a local phQTL targeting Ira2 itself, a phQTL targeting Cdc25, which promotes the GTP-bound form of Ras2, and a phQTL targeting the Ras2 protein  (Jian et al, 2010). In addition, we detected phQTLs for numerous downstream targets of Ras/Pka signaling at the IRA2 hotspot (Dataset EV5-S10). Hence, our data reveal widespread differences in phosphorylation states of key signaling molecules in the Ras/Pka pathway upon genetic variation and therefore contribute to a better understanding of the molecular mechanisms through which this QTL hotspot affects cellular phenotypes. Another example illustrating the effect of genetic variants on cellular signaling is the HAP1 hotspot on chromosome 12 (chrXII:2)  (Brem et al, 2002; Albert et al, 2014) , which was previously shown to affect cellular growth in multiple conditions  (Bloom et al, 2013).
Hap1 regulates respiration in response to oxygen and iron deprivation  (Verdi ere et al, 1985). Again, our analysis revealed that this locus has much more prominent effects on the phospho-layer than on transcript or protein abundance  (Fig 4A). To determine a potential regulator through which the HAP1 hotspot affects the phosphorylation state of these proteins, we leveraged a previously published dataset of strains with deleted kinases and phosphatases (Bodenmiller et al, 2010) (see Materials and Methods for details). Among all tested regulators, targets of the HAP1 hotspot were most significantly enriched in proteins changing their phosphorylation upon perturbation of the protein kinase PSK2, a known regulator of carbohydrate metabolism  (Rutter et al, 2002). Strikingly, we found Psk2 itself to be a target of the HAP1 hotspot both on the transcript and on the phosphorylation layer.
Together, these findings suggest a signaling cascade from the HAP1 locus via altered Psk2 activity to Psk2 substrates. The added value of integrating protein phosphorylation with transcript and protein abundance was even more apparent from a hotspot on chromosome 8 (chrVIII:1). GPA1, a key factor in the pheromone response pathway, was earlier shown to harbor a variant that explains effects on some targets of this hotspot  (Yvert et al, 2003). The hotspot had a disproportionate effect on phospho-traits, harboring only 69 eQTLs, 8 pQTLs, and 1 ptQTL, but 41 phQTLs (Dataset EV5-S2, Fig  3B ). Many of the hotspot eQTL targets (41%) could be found downstream of the mating pheromone response pathway, which supports GPA1 as a causal gene. In addition, the hotspot region included the kinase STE20.
Ste20 is the key activator of multiple mitogen-activated protein kinase (MAPK) pathways, including not only the mating pheromone response but also the invasive growth regulation, regulation of sterol uptake, and osmotic stress response  (Klipp & Liebermeister, 2006 ) (Fig  6A ). STE20 contains variants that co-segregate with those in GPA1, making it a second candidate for a causal gene for this hotspot. Indeed, for up to 60% of the target genes of this hotspot we found a link to Ste20 and/or Gpa1 (Dataset EV5-S11 and Dataset EV5-S12). Further, most of the phosphopeptides and phRes traits targeted by the hotspot corresponded to proteins phosphorylated by components of MAPK pathways downstream of Ste20 (51%; Fig 6A, Dataset EV5-S11 and Dataset EV5-S12). GO enrichment analysis of targets of the hotspot at the transcript layer revealed enrichment for biological processes under the influence of STE20 (Dataset EV5-S13).
Discussion In this study, we present the first dataset that concurrently quantifies the effects of natural genomic variation on the transcript (eQTL), protein (pQTL), and protein phosphorylation (phQTL) layers in the same sample set. Whereas previous "omics" QTL studies have primarily focused on the genetic effects on the abundance of transcripts and proteins, our study integrates abundance with the phosphorylation states of proteins. This allows us to associate the consequences of genomic variation not just on transcript and protein abundance changes but also on phosphorylation-mediated activation of signaling pathways and regulatory networks. The detailed dissection of the three hotspots around IRA2, HAP1, and GPA1/ STE20 exemplifies this notion. Quantifying transcript and protein levels along with the activity state of the proteins indicated by phosphorylation from the same yeast culture maximized the comparability of the data, thus facilitating data integration.
Thus, in order to interpret genomic variants affecting complex traits (including loci identified in human GWAS), it would help to leverage comprehensive catalogues of functional phosphorylation sites  (Needham et al, 2019; Ochoa et al, 2020). In conclusion, our study is a first step toward a better mechanistic understanding of genomic effects on multi-layered cellular networks and physiological traits. This study captures some of the enormous complexity of how genetic variants influence signal transmission between different molecular layers. While the data here were obtained from just one standard growth condition, mapping phosphorylation states will be even more critical when aiming to understand the effects of genetic variation in the context of dynamic environmental cues. In addition, future work should address additional aspects of protein state changes, including other PTMs, changes in protein folding, protein complex formation, and proteinligand interactions.
Sample preparation All media were prepared in a single batch to limit experimental variability. The BYxRM yeast strain collection, which we obtained from Rachel Brem, was originally derived from a cross between the two parental strains, BY4716, an S288C derivative (MATα lys2Δ0), and RM11-1a (MATa leu2Δ0 ura3Δ0 ho::KAN)  (Brem et al, 2002). A subset of 129 strains were picked in random series of 16 (Dataset EV3), pre-cultured in in-house made synthetic dextrose medium (S. D. , containing per liter: 1. 7 g yeast nitrogen base without amino acids (Chemie Brunschwig), 5 g ammonium sulfate, 2% glucose (w/v), 0. 03 g isoleucine, 0. 15 g valine, 0. 04 g adenine, 0. 02 g arginine, 0. 02 g histidine, 0. 1 g leucine, 0. 03 g lysine, 0. 02 g methionine, 0. 05 g phenylalanine, 0. 2 g threonine, 0. 04 g tryptophan, 0.
RNA extraction Total RNA was isolated from deep frozen aliquots of yeast pellets using the RiboPure™ RNA Purification Kit, yeast (Ambion), which includes a DNase treatment to eliminate contamination. RNA quality was assessed using RNA ScreenTape assay (Agilent). All RNAs were of very high quality (Dataset EV3 median RIN 9. 8, minimal RIN 9. 1). RNA-seq library preparation cDNA libraries were prepared from poly(A) selected RNA applying the Illumina TruSeq protocol for mRNA using a total of 1 μg RNA per sample and 14 PCR cycles.
RNA-seq-based genotyping The BYxRM yeast cross has been widely used for QTL mapping. Microarray-based genotype information at 2,957 markers is available for the segregants  (Brem & Kruglyak, 2005). However, deep sequencing enables a more accurate genotyping of recombinant lines. To this end, we have exploited published resequencing data of the parental strains  (Bloom et al, 2013)  together with the RNA-seq data generated here to infer the genotype of the segregants of the BYxRM cross using a method that we previously developed (Cl ement-Ziza et al, 2014). Resequencing data of the parental strains Sequence variation information of the parental strains was obtained from http://genomics-pubs. princeton. edu/YeastCross_BYxRM/; only calls with a MQ ≥ 30 were considered, which represents 42,769 polymorphic sites.
Genotype inferring from RNA-seq data RNA-seq data were mapped to the S. cerevisiae reference genome (SaCer3) using TopHat 2 (Ref. Kim et al, 2013)  with the following options: --min-intron-length 10 --min-segment-intron 10 --b2-very-sensitive --max-multihits 1 --librarytype fr-secondstrand. Read group information was added, and BAM files were sorted using Picard utilities (http:// broadinstitute. github. io/picard/). RNA-seq data were further processed using the GATK pipeline (version 3. 4-46-gbc02625) following the best practice guide (Van der Auwera et al, 2013). First reads containing exon-exon junction were split using the following options: -T SplitNCigar-Reads -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS; then, the variants were called using the UnifiedGenotyper at the 42,769 polymorphic sites identified in the resequencing data of the parental strains.
Filtering and missing value imputing For every polymorphic site between the parental strains, we compared the polymorphisms in the segregants and the parental strains to infer which allele was inherited. We further excluded polymorphisms (i) that could not be called (or correctly called) in the parental strains based on RNA-seq data, (ii) that could be called in < 70% of the segregants, and (iii) with a lower allele frequency of < 20%. As previously discussed (Cl ement-Ziza et al, 2014), genotypes called differing from the two direct flanking markers in more than one segregant (93 cases) probably denote erroneous genotype calls; the corresponding polymorphisms were excluded from the analysis. This resulted in 25,590 polymorphisms that were considered as genetic markers. Finally, missing genotype values were inferred from the two neighboring polymorphisms if those were each within 20 kb of the polymorphism of interest and were inherited from the same parental strain (Cl ement-Ziza et al, 2014) (Dataset EV2).
5. 0a; Dobin et al, 2013)  with the following options: --alignIntronMin 10 --quantMode GeneCounts. The gene-specific read counts (strand-specific) generated by STAR were used to quantify gene expression. RNA-seq coverage was computed by dividing the sum of the length of all reads by the sum of the length of the coding regions of the quantified transcripts per sample. Raw read counts were normalized using the rlog method of DESeq2  (Love et al, 2014). Normalized data were further corrected for effects due to culture batches using the non-parametric empirical Bayes framework Com-Bat  (Johnson et al, 2007)  (Dataset EV1). Normalized and batchcorrected read counts were corrected for gene length as follows for gene i in sample j: c 0 ij ¼ log 2 2 cij Á 1000 l i ; where l i is the length of the coding region of gene i, excluding intronic regions.
LC-MS data acquisition The peptide concentration in all samples was measured on a Nano-Drop at OD280 and normalized to allow injection ~1 µg of material into the mass spectrometer. The samples were randomized and then either injected individually for SWATH-MS acquisition or pooled and injected in technical duplicates for shotgun acquisition. The LC-MS acquisitions were performed on an AB Sciex 5600 TripleTOF coupled to a NanoLC-2D Plus HPLC system (for the main QTL dataset), or on an AB Sciex 6600 TripleTOF coupled to a NanoLC-1D Plus HPLC system (for the STE20/GPA1 validation). The liquid chromatographic separation and mass spectrometric acquisition parameters were essentially as described earlier  (Selevsek et al, 2015). The peptide separation was performed on a 75-µm-diameter PicoTip/PicoFrit emitter packed with 20 cm of Magic C18 AQ 3 resin using a 2-35% buffer B at 300 nl/min (buffer A: 2% acetonitrile, 0.
, 1,174-1,200 m/z. The collision energy for each window was determined based on the collision energy for a putative doubly charged ion centered in the respective window using the same equation as above with a collision energy spread of 15 eV. All the MS data files were visually inspected and curated at this stage for low total ion chromatogram intensities, and the corresponding samples were reinjected when possible. For the main QTL dataset, this resulted in a final set of 179 SWATH data files for nonenriched and 179 SWATH data files for phospho-enriched samples (Dataset EV3) that were used for data extraction. Similarly, 40 DDA files for non-enriched and 30 DDA files for phospho-enriched samples were selected for database searching and library generation. The STE20 validation dataset additionally contained 16 phosphoand 16 non-phospho-data files.
LC-MS database searching The shotgun data were searched with Sorcerer-Sequest (TurboSequest v4. 0. 3rev11 running on a Sage-N Sorcerer v4. 0. 4) and Mascot (version 2. 3. 0) against the SGD database (release 03 Feb. 2011, containing 6,750 yeast protein entries, concatenated with 6,750 corresponding "tryptic peptide pseudo-reverse" decoy protein sequences). For the search, we allowed for semi-tryptic peptides and up to two missed cleavages per peptide. For the non-enriched samples, we used carbamidomethylation as a fixed modification on cysteine residues and oxidation as variable modification on methionine residues. For the phospho-enriched samples, we additionally allowed for phosphorylation as variable modification on serine, threonine, and tyrosine residues. The Sequest and Mascot search results were converted to pep.
Those "split peptide assays" were given different protein entry names labeled Subgroup_0_ProteinX to Subgroup_N_ProteinX, respectively. The fragment ion coordinates for the peptides contained the top 6 most intense (singly or doubly charged) y or b fragment ions for each spectrum, excluding those in the SWATH precursor isolation window for the corresponding peptide. The nonenriched assay library comprised assays for 19,473 peptides (thereof 18,074 proteotypic peptides matching a total of 3,119 unique proteins). The phospho-enriched assay library comprised assays for 14,339 peptides (thereof 12,969 phosphopeptide sequences) or assays for 13,786 proteotypic peptides (thereof 12,678 proteotypic phosphopeptides, matching a total of 1,676 unique phosphoproteins). The same phospho-and non-phospho-libraries were used to extract both the QTL and the STE20 validation datasets.
The SWATH-MS data extraction was performed separately for the main QTL and the STE20 validation datasets using the iPortal workflow manager  (Kunszt et al, 2015)  calling OpenSWATH (openMS v. 1. 10) (Ro ¨st et al, 2014) and pyProphet  (Teleman et al, 2015). The precursors were then realigned across runs using TRIC  (Ro ¨st et al, 2016). The two resulting SWATH identification result files contained a total of 18,273 identified peptides (thereof 16,922 proteotypic peptides matching a total of 2,940 proteins) for the nonenriched datasets; in the phospho-enriched datasets, there were 13,748 identified peptides (thereof 12,412 phosphopeptides) or 13,218 proteotypic peptides (thereof 12,139 proteotypic phosphopeptides matching a total of 2,247 unique phosphoproteins). After alignment, we used a set of in-house scripts to compare the chromatographic elution profiles of the various isobaric phosphopeptide isoforms matching a same delocalized peptide form (peptide sequence + number of phosphorylations) within each single run and to eventually group those co-eluting phosphopeptide assays into the proper corresponding number of phospho-peak clusters (labeled _ cluster0 to _clusterN, respectively).
First, the genotype matrix was normalized as described in  Patterson et al (2006)  (see equations (1) to (3) there). Then, we carried out a singular value decomposition on the normalized genotype matrix followed by an eigenvector decomposition. We then selected those eigenvectors corresponding to the top seven eigenvalues as covariates for the QTL mapping (additional predictors for growing the Random Forests). These first seven vectors explained more than 25% of the genotype variance (Dataset EV2). We adapted the approach described in Cl ement-Ziza et al (2014) using a different score describing the importance of each predictor in the Random Forest (variable importance measure, VIM). While previous work relied on the selection frequency as the VIM, i. e. , the number of times a predictor was used in a forest, to quantify the importance of each predictor, we combined two previously published VIMs, RSS and PI  (Liaw & Wiener, 2002) , to compute an informative combined score serving as the VIM in our study.
QTL hotspot detection To formally identify hotspots in our dataset, the genome was divided into 40-kb bins (293 bins, bins at chromosome extremities could be bigger). As proposed before  (Brem et al, 2002) , if the linkages were randomly distributed across the genome, the number of linkages in each bin would be expected to follow a Poisson distribution with the mean of N linkage =N bin. We use this distribution to estimate the highest number of linkages that a bin could contain at a probability lower than 0. 01. These numbers were 30 for eQTLs, 13 for pQTLs, 9 for ptQTLs, 11 for phQTLs, and 4 for phResQTLs. Bins with a sufficient number of QTLs at a given molecular layer were considered as QTL hotspots. Consecutive bins on the same chromosome with enough QTLs were combined into single hotspots.
The matrix containing the read counts of all samples where RNA extraction was possible (all except BY:GPA1-BY:STE20-3) was processed with the function DESeqDataSetFromMatrix to generate a DESeq dataset. The DESeq function was run on this object using the following design: ~batch + strain, where batch refers to the RNAextraction batch. For each pair of strains, we used the results function from DESeq2 to compute log2 fold changes and P-values for differential expression.
Acknowledgments JG and CS received funding from the  German Science Foundation  ( DFG CRC 680 & CRC 1310 ). MCZ and CS received funding from the  German Federal Ministry for Science and Education (BMBF)  as part of the e: Med program  (grant no. 01ZX1406 ). OTS received funding from the  German Academic Exchange Service (DAAD) , the  Swiss National Science Foundation , and the  Human Frontiers Science Program. IB was supported by the  Swiss National Science Foundation  (grant no. 31003A_166435 ). The group of RA was supported by the  European Research Council  ( Proteomics4D  ( AdvG grant   670821  and Proteomics v3. 0; AdvG-233226). We thank  Dr. Rachel Brem  for providing the strains of the BY/RM cross. We thank the  Regional Computing Center of the University of Cologne (RRZK)  for providing computing time on the DFG-funded HPC cluster  CHEOPS and support.
Data availability The datasets and computer code produced in this study are available in the following databases: • RNA-seq reads: ArrayExpress (E-MTAB-8146, https://www. ebi. ac. uk/arrayexpress/experiments/E-MTAB-8146/) • MS data: Pride (PXD010893, http://proteomecentral. proteomexchange. org/cgi/GetDataset?ID=PXD010893) • Analysis scripts: GitHub (https://github. com/Jan88/phosphoQTL) Expanded View for this article is available online.
All subjects had normal or corrected-to-normal visual acuity, as checked with the Freiburg Visual Acuity Test  57. One VGP has been removed from the EEG analysis due to excessive noise in the recorded signal. All participants were given a compensation of 45€ for their participation and provided a written informed consent to take part in the study. The study has been approved by the local ethics committee of the University of Strasbourg (i. e. Comité d'Éthique de Recherche). Experimental Protocol. The experiment used the Unity software (Unity technologies, v. 2019. 3. 9f1) to create the virtual environment. The HTC Vive Eye Pro (HTC Corp. ) headset and controllers were used to immerse the participants in VR. Participants wore both the EEG and VR headsets while sitting on a chair. The use of VR has several advantages.
Firstly, VR allows a better trade-off between fully-controlled experimental settings and ecological experience (e. g. 3D visual percepts) in comparison with 2D screen setups. Secondly, it increases the engagement of the participant to the task at hand. Thirdly, the embedded eyetracking system to the VR headset allows researchers to easily track and record the gaze in the 3D space. Eye-tracking was used to trigger visual stimuli. Participants were instructed to fixate the warning signal without moving their eyes, and it was only after a time interval free of saccades and eye-blinks that the warning signal was switched off. This procedure improves the data quality of EEG recordings. Each participant was immersed in a virtual room, facing four 3D robots, each with a light whose color and onset were manipulated (Figure  1A ).
The experiment was composed of two intertwined tasks: a variable foreperiod task and an asynchrony detection task. The asynchrony detection task consisted of discriminating whether the lights of two robots appeared synchronously or asynchronously (using delays of 11 ms, 33 ms, or 66 ms). Participants performed four blocks of trials for each task. The participant switched tasks every block to reduce boredom. At the end of each experimental block, a break was proposed to the participant to remove the VR headset. Given the research questions investigated in this article, only data collected from the variable foreperiod task are presented. Variable foreperiod task and stimuli. In case the button was pressed before the go-signal, a warning sound was delivered to the participant signaling the incorrect response.
Second, raw data (i. e. the 3D gaze position over time) of each trial were interpolated with a spline method to increase the temporal precision followed by the calculations of the derivations of the speed of vertical and horizontal movements. Saccades were detected using a modification of a published algorithm  60  based on gaze's velocities. A threshold criterion for saccades detection was determined in the 2D velocity space based on the horizontal and the vertical velocities of gaze movement. This 2D space represented a plane surface located at the stimulus area (i. e. the warning signal and target). This threshold was represented by a 2D ellipse. For each trial, we set the threshold to be six times the SD of the gaze velocity 37,39,61 using a median-based estimate of the SD. Small fixational saccades were defined when six or more consecutive velocity samples (i.