dataesr/scientific-paragraphs-categorization huggingface dataset:
DISCUSSION In this manuscript, we describe BrainVar, a unique resource of paired genome (WGS) and transcriptome (bulk tissue RNAseq) data derived from 176 human DLPFC samples across prenatal and postnatal development (Figure  1 ). We identify 23,782 genes expressed during human cortical development, gene lists relating to developmental trajectories and coexpression, and common variants that alter gene expression (eQTLs). Our analyses show how these datasets relate to each other and to gene expression in cell types derived from single-cell RNA-seq data and to CNS traits and disorders derived from genomic analyses (exome sequencing and GWAS). In addition to developing a resource with utility for future studies of human development, neurobiology, and neuropsychiatric disorders, we also describe key biological insights, including the nature of the late-fetal transition in gene expression (Figures  2  and 3 ), an early-fetal transition (Figures  2  and 3 ), developmental processes and cell types implicated in CNS traits and disorders (Figure  4 ), eQTLs split by effect size across development (Figure  5 ), differing characteristics of genes with constant versus temporal-predominant eQTLs (Figure  5 ), and the application of this dataset to implicate specific genes at GWAS loci (Figure  6 ).
Principal component analysis identifies developmental age as the most important factor underlying the variance in gene expression in this dataset. The majority of this temporal variance occurs in two transitional phases (Figure  2 ), the early-fetal and late-fetal transitions. The early-fetal transition is a coordinated decrease in expression of multiple genes in early development (epoch 0; periods 1-2; 6-10 post-conception weeks) that coincides with the establishment of regional identity across the brain. Concordant with a possible role of the early-fetal transition in this process, the expression of several genes associated with noncortical tissues (e. g. , UNCX and DMBX1) is decreased during this period. In addition, we found that the early-fetal transition is captured in the M8 red module (Figure  3 ), which is enriched for lncRNA transcripts and Gene Ontology terms related to morphogenesis and cell fate.
Higher-resolution datasets across development, including single cells, additional brain regions, and larger sample sizes, along with complementary analyses of brains of individuals with neuropsychiatric disorders and rare genetic disorders, are likely to provide additional insights. The combination of genomic and transcriptomic data across development allows us to interrogate human cortical development from a molecular perspective at a higher resolution than before. Understanding patterns of temporal and cell type specificity, along with eQTL colocalization to resolve GWAS loci, has already provided insights into the pathology underlying neuropsychiatric disorders. Further delineation of these patterns is likely to be critical for a detailed understanding of etiology as a foundation for therapeutic development.
QUANTIFICATION AND STATISTICAL ANALYSIS WGS variant calling-Using the pipeline from the Centers for the Common Disease Genomics project  (Regier et al. , 2018) , FASTQ reads were aligned to the GRCh38 reference from the 1000 Genomes Project using BWA-MEM version 0. 7. 15  (Li and Durbin, 2009). Reads were sorted and duplicates were removed with Picard, version 2. 17. 5 (https:// github. com/broadinstitute/picard/); base quality score recalibration was then performed with the Genome Analysis Toolkit (GATK), v3. 8-0-ge9d806836  (McKenna et al. , 2010). Variant calling and joint genotyping were done with the Haplotyper and Genotyper tools from Sentieon v201711. 01, a toolkit containing modules that are mathematically equivalent to their counterparts in the GATK  (Freed et al. , 2017). SNP and indel recalibration were performed on the joint genotyped VCF file.
RNA-seq alignment and gene-level read count quantification-RNA-seq reads were aligned to the human genome (hg38/GRCh38) using STAR aligner  (Dobin et al. , 2013)  and gene-level read counts were calculated using HTSeq  (Anders et al. , 2015)  based on GENCODE v21 annotation  (Harrow et al. , 2012). Read counts per gene were then converted to counts per million (CPM), which were logarithmically scaled to base 2 (log 2 CPM). Of the 60,155 genes assessed, 23,782 were defined as being cortically expressed, based on CPM ≥ 1 in at least 50% of samples of either sex in at least one of the 12 developmental periods (Table  S2 ). For these 23,782 genes, the median log 2 CPM ranged from -5. 9 to 12. 3, with a median of 2. 3.
DCX is involved in neuron migration, and is expressed most strongly during prenatal development, with distinctly decreased postnatal expression. All 176 samples showed the expected DCX expression levels given samples' reported age. Similar results across all expressed genes were obtained using principal component analysis (below). Ancestry estimation-Ancestry was estimated using principal component analysis of common SNPs and indels in the WGS data, run alongside 3,804 additional individuals of known ancestry with WGS data (parents from the Simons Simplex Collection;  An et al. , 2018). From 10,688,106 SNPs and indels with allele frequency ≥ 5% in either this dataset, one of the three batches of Simons Simplex Collection data, or GnomAD genomes, variants were pruned for independence with linkage disequilibrium r 2 < 0.
Comparison between BrainVar and BrainSpan-While the BrainVar bulk tissue RNA-seq dataset catalogs gene expression in the DLPFC of 176 donors across development (6 post-conception weeks to 20 years of age), BrainSpan catalogs gene expression across 16 brain regions, including DLPFC, in 40 brains ranging from 8 post-conception weeks to 40 postnatal years  (Li et al. , 2018). Of note, data were generated for 14 brains in both BrainVar and BrainSpan. To compare the data from the 14 samples profiled in both BrainSpan  (Li et al. , 2018)  and BrainVar, we reprocessed the entire BrainSpan dataset using the BrainVar RNA-seq analysis pipeline, transforming raw read counts to log 2 CPM values. We then filtered to genes that passed BrainVar minimum expression criteria and calculated sampleto-sample and geneto-gene Pearson correlation coefficients.
We performed principal component analysis on a covariance matrix of the 14 commonly sequenced samples and the 23,782 BrainVar-expressed genes for both datasets. Despite differences in library preparation, with BrainSpan using poly-A priming compared to TruSeq random priming in BrainVar, gene expression was highly correlated per-sample and per-gene. BrainSpan also generated ChIP-seq datasets of histone 3 lysine 27 acetylation (H3K27ac), a marker of active genes, in DLPFC samples  (Li et al. , 2018). These H3K27ac data mirror the expression profiles in BrainVar, with Falling genes enriched for fetal-biased H3K27ac peaks and depleted for adult-biased H3K27ac peaks, while Rising genes show the opposite pattern of enrichment (Figure  S2 ). Transcriptome temporal trajectory estimation-Statistically, the temporal dynamics of the expression of genes can be modeled as a mixture of K distinct trajectories, each with Gaussian noise  (Jones et al.
Gene ontology functional enrichment for temporal trajectories-For functional enrichment, we characterized genes sets for each trajectory using the R package, gProfiler  (Reimand et al. , 2011). The pathway enrichment test was performed using Gene Ontology Biological Process terms, which contain between 10 and 2,000 genes, and all 23,782 cortically expressed genes were used as background. Enrichment tests were subject to the "moderate" hierarchical filtering parameter, and the FDR multiple correction in the gProfiler (Figure  2F ). Assessing enrichment in tissue-specific genes from GTEx-To assess enrichment across tissues, 27,546 transcripts with an RPKM ≥ 0. 5 in 80% of samples one or more tissues in GTEx (https://gtexportal. org/home/) were defined and log-transformed (log 2 [RPKM+1]). For each gene, expression between each tissue and all other tissues was assessed using a moderated t-test (R package limma), with models adjusted for age, RIN, gender, and surrogate variables.
The Benjamini and Hochberg method was used to estimate false discovery rate (FDR) and tissue enriched genes were defined as: log fold-change > 0. 5 and FDR < 0. 05. The enrichment of the Falling, Non-transitional, and Rising genes was assessed with using Fisher's exact test with 23,782 cortically expressed genes as a background (Figure  2 ). Identifying genes enriched in cell types from single cell data-To identify cell type-enriched genes we calculated a tau metric (Kryuchkova-Mostacci and Robinson-Rechavi, 2017) from the per gene log 2 TPM of log 2 UMI values for genes within cell type clusters in prenatal forebrain  (Nowakowski et al. , 2017)  and postnatal cortex  (Li et al. , 2018). Each gene was also ranked across all clusters within the dataset on the basis of TPM/ UMI, so that the cluster with the highest TPM/UMI was ranked as "1," while the cluster with the second highest TPM/UMI for that gene was ranked as "2," etc.
The genes were sorted by TPM/UMI rank and visually inspected for cell type specificity or substantial enrichment in single cell RNA-seq data from the prenatal and postnatal human cortex (https://cells. ucsc. edu/?ds=cortex-dev#). For each cell type, the top ten genes that showed clear specificity/enrichment were selected. Enrichment of gene trajectories in temporal putative cis-regulatory elements -H3K27ac peaks present in more than two samples of fetal or adult dorsal frontal cortex in BrainSpan  (Li et al. , 2018)  were tested for fetal versus adult temporal bias using DESeq2  (Love et al. , 2014). Temporally biased genes were defined as adjusted p < 0. 01 and fold change ≥ 2. A category of "non-temporal" H3K27ac peaks was generated with peaks showing p > 0. 05. All peaks were annotated using the gene with the closest transcription start site in Gencode v21.
Genes were classified as only-fetal or only-adult if they were associated with fetal or adult-only H3K27ac peaks, respectively. Enrichments of each category of H3K27ac-genes in each category of eGenes were tested by means of a Fisher Exact's test and p-values were adjusted using Benjamini-Hochberg, using genes associated to non-variant H3K27ac peaks as a reference background. WGCNA network construction and module definition-To assess the functional topology in cortical samples, we applied Weighted Gene Co-Expression Network Analysis (WGCNA)  (Langfelder and Horvath, 2008)  to 23,782 cortically expressed transcripts. Network analysis was performed with WGCNA (version 1. 63) using a signed network, choosing a soft-threshold power, the mean connectivity less than 50, and scale-free topology greater than 0.
, 2011) , as described above in the analysis of temporal trajectory genes. To identify WGCNA module genes that are regulatory targets, we searched for transcription factor binding targets using the ChEA  (Lachmann et al. , 2010; Satterstrom et al. , 2020) , and TRANSFAC  (Matys et al. , 2003) , and microRNA using the mirTarbase database  (Chou et al. , 2018). WGCNA module preservation-To assess whether 19 co-expression modules in our samples were preserved in other, independent DLPFC or frontal cortex expression datasets, we compared our dataset with non-overlapping samples from the BrainSpan dataset  (Li et al. , 2018)  and applied the module preservation function from the WGCNA R package  (Langfelder et al. , 2011). From the BrainSpan dataset, we selected DLPFC (n = 30), 10 non-DLPFC neocortical regions (n = 317), and subcortical regions excluding the cerebellum (n = 140).
Given the original co-expression network constructed above, we used modulePreservation to calculate module preservation statistics from 100 permutations (Table  S3 ). Clustering analysis in protein-protein interaction network-To examine functional association of a group of genes/proteins, we performed a clustering analysis of a proteinprotein interaction network implemented in SANTA R package  (Cornish and Markowetz, 2014 ). Ripley's K function provides a measure of whether points are clustered together or randomly dispersed (homogeneous) in a network and the SANTA R package reformulated Ripley's K function for a protein-protein interaction network. Clustering can be indicative of functional association between the genes/proteins under consideration. To test departure from homogeneity of a given gene set, we drew an empirical null distribution of clustering from 1,000 random samples of matching sized gene sets from the BioGRID protein-protein interaction network data  (Stark et al.
, 2006; Winter et al. , 2011) (v3. 4. 132). We reported departure from null distribution as a Z-score (i. e. , Z > 0: level of clustering of a gene set greater than null expectation and Z < 0: level of clustering of a gene set less than null expectation). Cis-eQTL detection and classification-Cis-eQTLs were identified for all high quality, common variants (N = 6,573,196) within 1 Mb of a gene boundary using the linreg function in Hail 0. 1, with period, sex, and the first five principal components of common variant ancestry as covariates. This analysis was run on three cuts of the BrainVar dataset: complete sample (N = 176, periods 1-12), prenatal-only (N = 112, periods 1-6), and postnatal-only (N = 60, periods 8-12). Separately for the results of each analysis, false discovery rate (FDR) was calculated for all gene-variant pairs using the Benjamini-Hochberg procedure.
We then classified all gene-variant pairs with FDR ≤ 0. 05 from at least one analysis into groups defined by the temporal specificity of their eQTL effects. To do this, we first identified one variant per gene with the smallest, FDR-significant p value, from any of the three analyses. We then used a Z-test to compare the regression coefficients for these variant-gene pairs from the prenatal and postnatal analyses: Z = β Pre -β P ost SE Pre 2 + SE P ost 2 Using the results from the eQTL analyses and from this prenatal-postnatal comparison, we then classified each of these top (one per gene) gene-variant pairs and their corresponding target gene (eGene) into one of five groups: 1. Constant eQTLs/eGenes, characterized by consistent effects across this developmental dataset: FDR ≤ 0. 05 in the complete sample analysis, same direction of effect and unadjusted p ≤ 0.
Using variant (GRCh38 position, reference, and alternate alleles) and gene (Ensembl gene IDs) identifiers, we matched significant variant-gene pairs separately from GTEx frontal cortex and prenatal brain to the variant-gene pairs meeting FDR significance (%0. 05) in the BrainVar analyses. For all significant eQTLs in BrainVar that overlapped with the reference datasets, we then compared the effect of the variant on the expression of its associated gene to determine the percentage of overlapping eQTLs with concordant direction of effect, as well as the Pearson correlation between the eQTL effects (beta from BrainVar, slope from GTEx or prenatal whole brain). The significance of this correlation was evaluated using the cor. test function in R. Distance between eQTLs and transcription start site-The distance between each significant eQTL and the transcription start site (TSS) of its associated eGene was calculated by comparing the variant position to the TSS position and strand of the gene according to Gencode v21.
At a Benjamini-Hochberg-adjusted p value ≤ 0. 05, we identified 8,094 genes with greater prenatal variance, and 1,752 genes with greater postnatal variance. Gene sets associated with CNS traits and disorders-To compare with 23,782 cortically expressed genes, we created the list of gene sets for previous trait and disorder association and functional properties. Gene identifiers were converted between studies based on the complete HUGO Gene Nomenclature Committee dataset. Autism spectrum disorder (ASD) risk genes were obtained from  Satterstrom et al. (2020) , an exome sequencing based gene discovery refining to high-confidence genes (n = 99) at a false discovery rate (FDR) ≤ 0. 1. Genes associated with developmental delay (n = 93) were selected from the exome sequencing analysis of the Deciphering Developmental Disorders project (Deciphering Developmental Disorders Study, 2017).
, 2015). We input the eGene lists from Prenatal-predominant, Postnatal-predominant, Constant, Prenataltrending, and Postnatal-trending (Table  S5 ) and summary statistics from published GWAS of schizophrenia (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014), autism spectrum disorder  (Grove et al. , 2019) , educational attainment  (Lee et al. , 2018) , multiple sclerosis  (Beecham et al. , 2013) , Alzheimers disease  (Lambert et al. , 2013) , triglycerides  (Willer et al. , 2013) , and height  (Wood et al. , 2014). First, SNPs from the GWAS summary statistics files were annotated to NCBI protein-coding genes that passed RNA-seq QC in our DLPFC expression data with a 10kb window on either side of the gene boundaries. Next, a gene-level analysis was done to determine the strength of association for each gene with phenotype of interest.
To assess whether genes in the eGene gene-sets are more strongly associated with the phenotype of interest than other genes, gene-based zscores are regressed on a gene-set indicator variable and MAGMA default covariates (gene size, gene density, sample size, 1/MAC, and the log of each of these). The beta coefficient for the gene-set indicator variable is tested for significance H A : β 1 > 0. Results from this analysis are reported in Table  S6. We did not observe significant enrichment of GWAS signal from any of the six phenotypes tested in any of the five temporally assigned eGene gene sets. However, we repeated a similar test, annotating GWAS summary statistics SNPs to NCBI protein-coding genes with 10kb flanking region that passed RNA-seq QC in our DLPFC expression data and had an assigned pLI score, and found that a gene-set defined as pLI score > = 0.
DATA AND CODE AVAILABILITY Open source scripts used in this study are referenced throughout. The pipeline for wholegenome analysis is available online at: https://github. com/sanderslab/psychcore-compute-platform. The accession number for the raw RNA-seq and WGS data reported in this paper, along with processed files, is PsychENCODE Knowledge Portal: syn21557948 on Synapse. org (https://www. synapse. org/#!Synapse:syn4921369).
ACKNOWLEDGMENTS Data were generated as part of the  PsychENCODE  Consortium, supported by  U01MH103339 ,  U01MH103365 ,  U01MH103392 ,  U01MH103340 ,  U01MH103346 ,  R01MH105472 ,  R01MH094714 ,  R01MH105898 ,  R21MH102791 ,  R21MH105881 ,  R21MH103877 , and  P50MH106934  awarded to  Schahram Akbarian  ( Icahn School of Medicine at Mount Sinai ),  Gregory Crawford  (Duke),  Stella Dracheva  ( Icahn School of Medicine at Mount Sinai ),  Peggy Farnham (USC) ,  Mark Gerstein (Yale) ,  Daniel Geschwind  ( UCLA ),  Thomas M. Hyde  (LIBD),  Andrew Jaffe  (LIBD),  James A. Knowles (USC) ,  Chunyu Liu  (UIC),  Dalila Pinto  ( Icahn School of Medicine at Mount Sinai ),  Nenad Sestan (Yale) ,  Pamela Sklar  ( Icahn School of Medicine at Mount Sinai ),  Matthew State  ( UCSF ),  Patrick Sullivan (UNC) ,  Flora Vaccarino (Yale) ,  Sherman Weissman (Yale) ,  Kevin White  ( UChicago ), and  Peter Zandi  ( JHU ).
This work was supported by funding provided by  Autism Science Foundation postdoctoral fellowships  (to D. W. and J. -Y. A. ) and a research award (to S. J. S. );  Simons Foundation Autism Research Initiative (SFARI)  grants  574598  (to S. J. S. ),  402281  (to S. J. S. , M. W. S. , B. D. , and K. R. ), and  573206  (to M. E. T. );  National Institute for Mental Health (NIMH)  grants  R01 MH109901  (to S. J. S. and M. W. S. ),  R01 MH110928  (to S. J. S. and M. W. S. ),  U01 MH103339  (to M. W. S. ),  R01 MH111662  (to S. J. S. and M. W. S. ),  U01 MH105575  (to M. W. S. ),  U01 MH106874  (to N. S. ),  P50 MH106934  (to N. S. ),  R01 MH109904  (to N. S. ),  R01 MH110926  (to N. S. ),  U01 MH116488  (to N. S. ),  R37 MH057881  (to B. D. ), and  R01 MH115957  (to M. E. T. );  National Institute of Child Health and Human Development (NICHD)  grants  R24 HD000836  (to I.
Rosat et al. (2009)  already performed a comparison of GNSS and SG measurements with hydrological model at the Strasbourg (France) station, but they have computed neither the spectral coherency nor the gravity-to-height changes ratio. In the present paper, using collocated gravity and displacement measurements, we will estimate the gravity (g) to height (u) change ratio dg/du that we compare with theoretical predictions from the MERRA2  (Gelaro et al. , 2017)  model. In the following, we first explain the methods used to compute the transfer function and the gravity-to-height change ratio. We then remind quickly about the data processing. Finally, we show the results for the hydrological surface loading predictions and for GNSS and SG observations at nine IGETS sites.
GNSS data processing We compute daily displacements from a global set of 117 worldwide stations using the GAMIT/GLOBK (v10. 6) software  (Herring et al. 2015). We use the latest tropospheric mapping function (VMF1;  Böhm et al. 2006)  with a priori values of zenith hydrostatic delay derived from the ECMWF meteorological reanalysis fields and residual Zenith Wet Delays (ZTDs) estimated at 2-hour intervals with two gradients per day. Standard solid Earth tide, ocean tidal loading (using FES2014a;  Carrère et al. 2016 ) and pole tide corrections are applied to follow standards recommended by the International Earth Rotation and Reference Systems Service convention (IERS2010)  (Petit and Luzum, 2010). We slightly modify the GAMIT software to include in our processing 3-hourly atmospheric and non-tidal oceanic loading effects at the observation level using ECMWF surface pressure field and TUGO-m  (Carrère and Lyard, 2003)  barotropic ocean model forced by air pressure and winds (see  Boy and Lyard, 2008; Gegout et al.
Summary We have verified the spectral coherency at seasonal time-scales between surface deformation and surface gravity changes recorded at nine collocated sites. The ratio of surface gravity perturbations over vertical height changes strongly depends on the local masses fluctuations. We have shown that the seasonal signal is coherent at eight sites over the nine studied here. The obtained dg/du ratios at the annual period are different between real observations and hydrological loading predictions. Influence of data pre-processing and local hydrological masses should be further investigated. Interpreting the obtained dg/du ratios in terms of rheological properties of the Earth would be the next step.
Data A structured questionnaire was used to obtain data on biosecurity measures on 34 and 93 farms from Catalonia and Galicia, Spain, respectively, that voluntarily participated in the project. The biosecurity questionnaire (available in Spanish upon request) was structured in 4 parts: (1) general data of the farm; (2) animal movements (e. g. , origin of the animals, frequency of introductions, test, quarantine facilities, external rearing farms, cattle fairs, pasture); (3) vehicles (e. g. , vehicles entering inside farm perimeter, vehicles coming in with other animals); and (4) visitors and staff (e. g. , external workers, frequency of visitors, use of protective clothing). The Autonomous Governments of Catalonia and Galicia provided records of cattle movements (national and international) for 2017.
Model Development To estimate the probability of introducing BVDV and BoHV-1 into a dairy cattle farm in a 12-mo period, a stochastic risk assessment model was developed using the mc2d package  (Pouillot and Delignette-Muller, 2010)  implemented in R (https: / / www. R -project. org/ ). Monte Carlo simulations (10,000 iterations) were performed, and all nonfixed input parameters were included as uncertain parameters.