dataesr/scientific-paragraphs-categorization huggingface dataset:
SNP analysis To estimate genetic divergence between the five accessions, SNPs were called from the RNA-Seq data in a twostep process. First, a set of high-confidence SNPs was identified within each accession separately. Second, these sets were evaluated in all five accessions to distinguish SNPs that are shared between all accessions (which reflect within-individual allelic polymorphisms that predate the most recent common ancestor of the accessions within the apomictic lineage) from SNPs that differentiate between accessions. To improve nucleotide coverage at all loci, trimmed and quality-filtered reads were concatenated per accession, pooling together biological replicates and both experiments. Using the BWA-backtrack program and default value parameters, resulting sequences were mapped to the Trinity reference transcriptome assembly  [51].
Mapping statistics were calculated with the FlagStat utility from SAM tools  [52]. Alignments were further processed using Picard tools (http://broadinstitute. github. io/picard) to add read group information for each accession, tag duplicates and to sort and index the final BAM files. The Freebayes analysis pipeline was used for SNP calling (Garrison and Marth, unpublished data) with the ploidy parameter set to three. Quality was filtered using the vcflib tool (downloaded from https://github. com/ ekg/vcflib). Read coverage for each SNP position in all accessions was checked with Bedtools  [53]. Due to the large amount of SNPs in the transcriptome, reflecting high intra-individual heterozygosity, we applied strict filters (Quality >200 and Read coverage >50 in each accession) to obtain a finalized set of high-confidence SNPs that are not shared between all five accessions.
This subset of high-confidence SNPs should accurately capture divergence of the accessions under asexual reproduction within the apomictic lineage. Using custom R scripts  [54] , a binary SNP matrix was built and processed to calculate pairwise Hamming distances between accessions  [55]  using the package APE  [56]. The hclust function generated the hierarchical clustering tree (from the fastcluster package  [57] ). Finally, in order to assess correlations between geographic, genomic and transcriptomic distances, Mantel tests were performed on the respective distance matrices using ZT software  [58].
95 Mb. Assembled transcripts were first annotated using protein databases (SwissProt and KEGG) and then using a nucleotide database of coding sequences from A. thaliana. Out of the 123,232 transcripts, 42,592 showed significant matches to SwissProt, 8,476 to KEGG and 39,685 to A. thaliana coding sequence database. Approximately 35 % of unique consensus sequences were successfully annotated in this study. This ratio is relatively low compared to other de novo transcriptome sequencing studies for non-model species due to short length of assembled transcripts. Lacking a genomic reference in addition with low-coverage RNA-Seq limits our ability to build full-length transcripts. For all libraries, ERCC spike-in reads were retrieved and mapped back to their known sequence reference. Approximately 30,000 reads (range ±12,000) were mapped for each library, representing less than 0.
Gene ontology and functional enrichment of DEG In DEG lists, 36. 5 % and 35. 5 % found a blast hit against A. thaliana and the SwissProt Database, respectively for the field and greenhouse experiments (Additional files 3 and 4). A large number of DEGs with descriptive functional annotations were retrieved as protein domains of transposable elements (34. 3 % of the annotated transcripts, Additional file 6). Most of these assembled transcripts were too short to have the full transposable element sequence. However, we detected protein domains of some Long Terminal Repeat retrotransposons (Class I transposable elements Copia and Gypsy) represented by GAG and POL genes (with Protease, Reverse Transcriptase, RNase H, Integrase and Endonuclease domains). Non-LTR retrotransposons were also well-represented with ORF2 and Reverse Transcriptase domains from putative LINE-1 and Tnt-1 elements.
Acknowledgments This work was supported by the  Dutch Organization for Scientific Research (NWO VIDI  grant  864. 10. 008  to KJFV). This is publication number 6032 of the  Netherlands Institute of Ecology (NIOO-KNAW ). We thank  Ingo Uhlemann ,  Bob Trávníček ,  Juhani Räsänen  and  Jan Kirschner  for proving seeds for this study. We also thank 3 anonymous reviewers for commenting on previous version of the manuscript and  Kelly Ramirez  for text editing.
In order to prioritise modules relevant to epilepsy we: (1) integrated modules with a large database of rare de novo mutations ascertained from patients with monogenic EE (and neurodevelopmental disease more generally); and (2) tested modules for enrichment of association to common forms of epilepsy using GWAS data. Utilising this integrative approach, we report a single network of 320 co-expressed genes genetically associated with monogenic and polygenic forms of epilepsy. Further analysis revealed functional disruption and/or downregulation of this network as a common mechanism regulating susceptibility to genetic and acquired epilepsy broadly, suggesting the network itself might be targeted as a novel therapeutic approach. As proof of concept for this, we show that among the many drugs capable of inducing transcriptional changes in cells  [28] , valproic acid (VPA), a widely used AED with a broad spectrum of clinical efficacy, is the one that is most significantly predicted to restore the expression of the network in epilepsy toward health.
Co-expression modules were inferred from these datasets using two approaches: (1) consensus Weighted Gene Coexpression Network Analysis (WGCNA) to identify coexpression modules common to all nine brain regions ('consensus modules'); and (2) DiffCoEx  [30]  to construct modules of differentially co-expressed genes in one or several brain regions compared to the other regions ('differential modules'). This identified 34 consensus modules and 13 differential modules, respectively (see Additional file 1: Table  S1 ). We then investigated whether these modules were reproducible in separate gene expression datasets. To this end, we used two unrelated publicly available brain expression datasets (see Additional file 2: Table  S2 ) and assessed module conservation using the Z summary value  [31]  (see 'Methods').
First, we used micro-array gene expression data spanning fetal and post-natal human brain development ('Brainspan' , GSE25219)  [32, 33] ; we observed that 32 of 47 (68%) modules were at least moderately preserved (Z summary >5) in at least one brain region during fetal development, and that these modules were generally more significantly (Z summary >10) preserved in cortical samples from post-natal brain (see Additional file 3: Figure  S2 ). We then analysed the preservation of co-expression modules using messenger RNA (mRNA)-sequencing (RNA-seq) expression data across different brain regions as well as in non-brain tissues (latter to assess specificity) from the Genotype-Tissue Expression (GTEx) project  [34, 35]. We observed that the brain co-expression modules from UKBEC data were generally poorly preserved in non-cortical tissue (e.
g. lymphocytes, fibroblasts) (see Additional file 3: Figure  S2 ), but well preserved using GTEx data relating to cortical brain regions; 20 out of 34 (59%) consensus modules and 8 out of the 13 (61%) differential modules were highly significantly preserved (Z summary >10) in at least one brain region from both the Brainspan and GTEx datasets. Finally, we assessed module preservation across-species using RNA-seq data from 100 healthy adult mouse hippocampus samples  [24] ; 20 out of 34 (59%) human consensus modules and 8 out the 14 (57%) differential modules were at least moderately preserved in the healthy mouse hippocampus (Z summary >5). In summary, taking into account the overlap of modules preserved according to these comparative network analyses, (See figure on previous page. ) Fig. 1  Schematic overview of study design.
We hypothesised that gene regulatory networks in the human healthy brain disrupted by de novo mutations ascertained from patients with epileptic encephalopathy could be informative for molecular processes involved in different types of epilepsy. We used post-mortem human brain samples ascertained from individuals with no history of psychiatric or neurological illness from the UK Brain Expression Consortium (UKBEC) dataset to build gene co-expression networks (modules) that are expressed across the whole brain ('consensus' modules, using the WCGNA consensus method), or differential to one brain region or another ('differential' modules, using DiffCoEx method). In order to prioritise modules relevant to epilepsy, we integrated modules with whole-exome-sequencing (WES) studies data of rare de novo mutations ascertained from patients with epileptic encephalopathy (EE) and neurodevelopmental disease more generally (ASD autism spectrum disorder, SCZ schizophrenia, ID intellectual disabilities, DDD broad developmental disease from the Deciphering Developmental Disorders study).
A single module was selected in this way: M30. Our hypothesis was that disruption of this gene network might lead to different types of epilepsy. M30 was therefore tested for enrichment of association to common forms of epilepsy using GWAS data from the International League against Epilepsy (ILAE) meta-analysis. Analysis of network genes' expression in disease in three epilepsies suggested functional disruption and/or downregulation of the network as a common mechanism regulating susceptibility to epilepsy broadly and therefore that the network itself might be targeted as a novel therapeutic strategy. As a proof of concept, we show that among the drugs capable of inducing transcriptional changes in cells of the Connectivity map (Cmap) dataset, VPA, a widely used AED with a broad spectrum of clinical efficacy, is the one that is most significantly predicted to restore the expression of M30 in epilepsy toward health.
Integration with rare de novo epilepsy variants De novo single nucleotide variants (de novo mutations, DNMs) were collated from published whole-exome sequencing (WES) studies of probands with monogenic EE and their unaffected parents (n = 356 trios)  [13, 14]. In view of the established genetic overlap between neurodevelopmental disorders  [25, 36, 37] , we extended the analysis to include WES trio datasets for autism spectrum disorder (ASD, n = 4186), schizophrenia (SCZ, n = 1004), ID (n = 192), and broad developmental disease from the Deciphering Developmental Disorders study (DDD, n = 1,133) (see 'Methods' for cohort references). For each co-expression module, we tested for the enrichment of DNMs, focusing on non-synonymous DNMs (consisting of all missense, nonsenses, and splice-site mutations) and including tests of enrichment for synonymous DNMs as a negative control.
Module M30 was significantly enriched in marker genes for interneurons (FET P value = 1. 65 × 10 -8 , OR = 4. 37) and pyramidal neurons (FET P value = 7. 71 × 10 -4 , OR = 3. 07). This mixed cell-type enrichment is consistent with the expression profile of M30 genes in the Allen Brain Institute single-cell RNA-seq brain dataset  [42]  (Fig. 3a ). We then explored the expression of M30 genes in different stages of human brain development utilising data from Kang et al. [32]  (see 'Methods') and observed a clear developmental gradient of expression of M30 genes beginning in early mid-fetal development (i. e. post-conception weeks  [16] [17] [18] [19] , maximal by birth and then persisting through all post-natal periods (Fig. 3b ). This clear developmental trajectory of expression of M30 led us to explore its transcriptional regulation.
Using the WebGestalt toolkit  [43]  to test for enrichment of transcription factor binding sites (TFBS), we found M30 was highly enriched for NRSF/REST (repressor element 1silencing transcription factor) targets (BH P = 2. 55 × 10 -8 , r = 8. 50). REST is a repressor of neuron-specific genes in early fetal development whose activity is downregulated in mature neurons  [44]  and dysregulated in a large array of brain pathologies including epilepsy  [45] [46] [47] [48]. The strong enrichment for REST-target genes in M30 is consistent with the tightly regulated developmental trajectory of expression of this module during brain development and supports the hypothesis that M30 genes non-randomly share a common regulation. We then investigated the physical interactions between the protein products of M30 genes.
First, selecting sources of protein-protein interaction (PPI) from the STRING database  [49] , we found highly significant enrichment for PPI among M30 genes (expected number of PPI = 120, observed number of PPI = 217, enrichment P = 1. 11 × 10 -15 ) (see 'Methods'). Then, we tested M30 for enrichment in PPI using the DAPPLE module  [50, 51]  of the GenePattern platform  [52]  and confirmed significant enrichment (expected Direct Edges Count (DEC): 59. 33, observed DEC: 180, empirical P = 0. 001). In addition, we collated non-redundant physical PPI data from three different sources: GeneMANIA  [53] , Hippie  [54]  and iRefWeb  [55]  (see 'Methods'). Using experimentally derived physical interactions that could be mapped to gene names using HUGO Gene Nomenclature Committee database  [56] , our final database contained information for 272,348 non-redundant protein interactions for 17,235 genes, 13,489 of which are expressed in the UKBEC dataset.
0044, FDR < 5%, Z-score = 2. 62) (see Table  1 , Additional file 5: Table  S4 ). None of the negative control GWAS datasets showed significant enrichment of association (see Additional file 5: Table  S4 ). A constituent gene within M30 is SCN1A, which is a known genome-wide significant susceptibility gene for epilepsy  [15, 54, 55]. We therefore re-tested the enrichment of association between M30 and epilepsy after excluding SCN1A from the list of M30 genes. After removing SCN1A from the set of M30 genes, the enrichment of association between M30 and epilepsy remained significant for both focal (FE: enrichment P = 0. 0019, Z-score = 2. 88) and generalised (GGE: enrichment P = 0. 0097, Z-score = 2. 34) epilepsy. In parallel, M30's enrichment for nonsynonymous rare variant DNM ascertained from patients Fig.
Relationship of M30 expression to broad forms of epilepsy We used genome-wide gene expression data from disease hippocampus tissues to investigate the association between expression of M30 genes and epilepsy. Differential expression analysis was carried out using gene expression data from: (1) the human epileptic hippocampus (24 TLE patients/23 controls) assayed by microarray; (2) the mouse epileptic hippocampus from the pilocarpine post-SE model of chronic TLE (100 cases/100 controls); and (3) the mouse hippocampus from a model of Dravet (SCN1A) disease (4 cases/ 4 controls) assayed by RNA-seq (see 'Methods'). GSEA  [57]  was applied genome-wide to the ranked list of gene scores (reflecting both the significance and the magnitude of expression changes in epilepsy, see 'Methods') to test for enrichment of differentially expressed genes in the set of M30 network genes.
Valproic acid reverses the downregulation of M30 We hypothesised that M30 represents a potential candidate network to target as a therapeutic strategy. We therefore aimed to identify small molecules whose effect on gene expression could restore the downregulation of M30 genes observed in epilepsy toward health. One approach to identifying drugs with this property is to computationally 'screen' drugs according to the extent and specificity of the overlap between genes differentially expressed by a drug and the component genes in the disease network. To date, the largest database reporting empirical changes in cellular gene expression following drug exposure is Connectivity Map (CMap), which reports transcriptional signatures for 1300 compounds  [28]. In order to screen drugs in CMap in terms of their predicted ability to restore M30 expression in epilepsy toward health, we first carried out differential expression analysis to identify genes with significant (FDR < 10%) drug-induced differential expression (DE) (see 'Methods').
Discussion The systematic integration of diverse genomic datasets, including brain gene expression data, de novo mutations for rare monogenic epilepsies and GWAS data from common epilepsy reported here, identified a convergence of genetic susceptibility for epilepsy on a generegulatory network (M30). The M30 network consists of 320 genes co-expressed across the human brain under tight developmental control. We observed that rare de novo mutations for monogenic epilepsy preferentially impact the mostly highly connected genes in the network. These results provide evidence for a shared pathogenesis between the rare and polygenic forms of epilepsy, although whether the enrichment of association of M30 with polygenic epilepsy reflects tagged rare variant contributions to the disease or true enrichment of common variants remains to be determined.
These findings suggested targeting the expression of M30 as a potential new therapeutic strategy. At the level of protein-protein interactions, several studies based on whole-exome sequencing trio analysis of neurodevelopmental disease have reported that genes with de novo loss-of-function mutations code for highly interconnected protein networks for autism  [63] , schizophrenia  [61]  and EE  [14]. In keeping with this, we found that the protein products of M30 genes are also significantly enriched for direct protein-protein interactions. In protein networks, the nodes with high degree (i. e. highly connected proteins), also termed hubs, are often functionally important. The most striking examples of this were initially discovered in worm, fly, yeast and E. coli proteomes, leading to the 'centrality-lethality' hypothesis which proposes that nodes with higher centrality in a network are more likely to produce disease phenotypes when removed compared to nodes with lower centrality  [64] [65] [66].
The analysis of differential expression of M30 genes in epileptic brain tissue versus controls showed a consistent enrichment in the M30 module for downregulated genes in epilepsy. Using quantitative data of seizure frequency in epileptic mice, we observed a significant enrichment in M30 for genes whose expression in the epileptic hippocampus was anti-correlated to the frequency of seizures (higher M30 expression correlated with fewer seizures). These analyses, together with genetic association of M30 to both monogenic and polygenic forms of epilepsy, suggest M30 may play a role in the mechanistic pathways leading to epilepsy. This led us to hypothesise that targeted increases in M30 expression could be exploited as a novel therapeutic approach in epilepsy. Using the CMap dataset, we prioritised several drugs and in particular VPA, a widely used AED, as having a preferential and therapeutic effect on M30 expression.
No other established AED had a significant effect on M30 gene expression in the CMap database and the identification of VPA as a candidate therapy for epilepsy was entirely independent of its known role in epilepsy or mechanism of action  [59, 71]  and based solely on an unsupervised analysis of disease and drug-related differential gene expression. Although VPA has a rapid onset of anticonvulsant activity, suggesting a direct effect on neuronal electrophysiology, the anticonvulsant activity of VPA has also been noted to increase during prolonged treatment  [72] , suggesting a potential transcriptional mechanism mediating its antiepileptic activity. These unsupervised analyses therefore suggest a systematic screen of drug-like molecules based on their measured effect on M30 gene expression might represent an efficient therapeutic strategy for new drug discovery in epilepsy.
Weighted gene co-expression network analysis of nine brain regions post-mortem human samples Brain region expression datasets The UK Brain Expression Consortium (UKBEC) released genome-wide gene expression data from several brain regions of individuals of European descent (GSE60862)  [80]. We selected samples for which data were available for the same 88 neuropathologically normal individuals in all the nine following brain regions: cerebellar cortex, frontal cortex, occipital cortex, temporal cortex, hippocampus, putamen, thalamus, medulla and white matter. Raw expression profiles from the Affymetrix Human Exon 1. 0 ST Array were processed to transcript-level expression with Affymetrix Power Tools (APT) (http://www. affymetrix. com/partners_programs/programs/ developer/tools/powertools. affx) using probe logarithmic intensity error (PLIER) normalisation  [81]  with probe GC-content correction.
Only the most reliable 'core' set of probes was used to generate transcript level expression profiles as defined in the Affymetrix chip definition file. Probes were retained if more than 50% of the samples in at least one brain region that had detection background P values < 0. 01. Gene-level expression was obtained by taking the median of the expression values of multiple exons mapping to the same gene. These filtering steps defined a final dataset of 20,388 probes, representing 15,199 protein coding unique genes (hg19/ GRCh37/Ensembl version 75), which were then used after quantile normalisation for network analysis and as the 'background' gene set for enrichment analyses. Before inferring gene co-expression networks, we used the probabilistic estimation of expression residuals (PEER) package  [82]  to determine hidden factors that explain much of the expression variability in the dataset.
We used principal component analysis (PCA) to calculate summary variables describing the variation in the microarray expression of the 15,199 genes and estimate the potential effects of the hidden factors identified by PEER and of the known clinical covariates on global gene expression variability. The plotting of hidden factors across samples allowed us to select the hidden factors F1 and F4 as potential batch effect without high relation to the variability explained by brain regions (see Additional file 3: Figure  S1 ). Gene expression levels were therefore adjusted to remove the effect of the selected hidden factors and covariates (F1, F4, age, gender, post-mortem interval, cause of death and brain-bank ID) by fitting linear models on gene expression using the lm function in R. The residuals from the linear model were then used as expression values (for 15,199 genes in 88 samples of each brain region dataset).
Co-expression network construction Genes were then grouped into modules using coexpression network approaches on the nine brain regions datasets: consensus method implemented in WGCNA to pinpoint co-expression modules common to the nine brain regions and DiffCoEx  [30]  to construct modules of differentially co-expressed genes in one or several regions compared to the others. The consensus co-expression network analysis was carried out using the blockwiseConsensusModules function in the WGCNA R package as previously described  [83, 84] , with the following parameters: β = 7 (chosen based on the scale free topology criterion r 2 > 0. 8), minModuleSize = 40, mergeCutHeight = 0. 25, minBlockSize = 20,000, cor-Type = 'bicor' (to compute robust pairwise correlations of gene expression using Tukey's biweight method  [85] ).