dataesr/scientific-paragraphs-categorization huggingface dataset:
Conclusions The accumulation of genome-wide studies has provided the opportunity for the analysis of the commonalities and idiosyncrasies in gene expression in host response to pathogen attack. Here, we synthesised 19 transcriptome datasets from experimentally infected honey bees and developed a new bioinformatics method, the directed rank product, to analyse gene expression profiles in order to identify the host specific responses to a diverse set of pathogens, and build a robust co-expression network. Although this method does not account for the amplitude of gene expression changes, using a rank product-based analysis has the advantage of enabling use of data from different platforms (e. g. microarray and RNA sequencing) in a single statistical analysis. Importantly, our analyses revealed a core set of genes involved in a common host response to phylogenetically distinct pathogens, yet also enabled identification of genes involved in pathogen specific host immune responses.
Dataset selection for meta-analysis We restricted our analysis to microarrays and RNA-seq datasets obtained from experimentally infected honey bee workers (Table  2 ). In total, we collected 19 transcriptome datasets obtained from nine experiments, reporting the differential expression of transcripts between control bees and samples parasitized by Nosema spp. , RNA virus and/or V. destructor and in which pathogen infection was a formal component of the experimental design (i. e. studies in which transcriptomes were generated for control and treatment groups). These 19 datasets were either from unpublished studies generated by the co-authors or recently published (and therefore publicly available) studies at the start of our work. Microarray probes and gene identifiers were converted or updated to the latest version of the honey bee genome assembly Amel_4.
5 and its annotation from NCBI  [63]. Differential gene expression data (treatment vs. control) were provided by authors of studies in terms of log 2 fold changes. All datasets were generated from worker honey bees experimentally infected by Varroa mites, RNA viruses and/or Nosema spp. , and for which gene expression was compared with uninfected control samples. Note that Varroa parasitism was also associated with high viral titers and therefore represented a 'Varroa plus virus' treatment. BQCV black queen cell virus, IAPV Israeli acute paralysis virus, DWV deformed wing virus. Categories (Cat. ) are N for Nosema, N/V for Nosema and virus co-infection, V for virus alone and M for Varroa mite ('Varroa plus virus'), as used across this study. Age and Days p. i. gives the age (i. e. days post-eclosion) and the number of days post infection when bees were collected for transcriptome analysis a Studies where honey bees were co-infected with two pathogens b This study used the model Sindbis virus expressing enhanced green fluorescent protein (SINV-GFP) c Transcripts from DWV (4 to 15 × 10 5 tags) and Varroa destructor virus (21 to 25 × 10 6 tags) present in brain transcriptomes from Varroa infested bees d Average proportion of reads attributed to DWV = 37.
6% (±14. 8 sem) e Average proportion of reads attributed to DWV = 47. 7% (±17. 7 sem) The use of one dataset (#3 in Table  2 ) required the reprocessing of the original raw data. We retrieved the pre-processed tiling array expression data (GSE25455) from NCBI GEO as described by Dussaubat et al. [26]. We then re-annotated the probe sequences of the tiling array by alignment to Apis mellifera transcripts extracted from Amel_4. 5 annotation as in Poeschl et al. [64]. We used the re-annotated probes to create sets of probes to measure the abundance of each transcript. We extracted the already computed log fold changes from the data files and applied quantile normalization. We used the new probe annotation to compute the median log 2 fold change of all probes assigned to represent a transcript.
Gene annotation Genes were annotated with GO terms using Blast2GO  [65]. The first step of sequence alignments was done inhouse using BLAST  [66]. The sequences of transcripts associated with gene identifiers from the honey bee genome assembly were recovered and compared to those in the non-redundant database  [67]  [downloaded on 2014/ 03/06, containing 35,149,712 sequences] using Blastx (parameters: e-value cutoff of 1E-6 and maximum number of alignments 20). Alignments were uploaded to the Blast2GO server and all following steps were done according to the Blast2GO pipeline using default settings. GoSlim-terms were chosen for annotation to reduce redundancy among overrepresented GO terms  [68].
Selection of genes for inclusion in the analysis We combined log 2 fold expression values of 11,165 genes from the 19 transcriptome datasets in one full synthesised dataset (fsd) file. Due to the diversity of expression detection platforms and changes in gene annotation within recent years, 56% of the genes contained missing values in at least one transcriptome dataset; hence, only 44% of the genes with complete observation across the 19 datasets would remain for the analysis. Therefore, to increase the number of genes to include in our synthesis, we constructed a restricted synthesised dataset (rsd) with a subset of the fsd containing log 2 fold expression values of 7,077 genes that had no more than three missing values (NAs) across the 19 transcriptome datasets. This enabled us to evaluate the expression changes of 63% of the annotated genes from the current genome assembly (see Additional file 1: Figure  S7 ) and constituted the dataset for further analyses.
Statistical analysis of gene expression required complete observations of a gene across all 19 transcriptome datasets. To overcome missing data for genes with incomplete observations, we ordered gene expression values by their log 2 fold change values and gave each a relative rank in each of the 19 transcriptome datasets. Relative ranks ranged between 1 for up-regulated genes and 0 for down-regulated genes. Missing values (5,015 of 134,463 in total) were then replaced by the average of non-missing relative ranks for the same gene from other datasets. We refer to this relative ranked dataset as ranked rsd, which is publicly available with the fsd at https://idata. idiv. de/DDM/Data/ShowData/35.
Rank product analysis Gene expression measurements using either microarrays or whole transcriptome sequencing (RNA-seq) vary greatly in methodology  [70] , resulting in substantial differences in the data produced. RNA-seq is generally more sensitive, producing gene expression levels spanning a greater dynamic range of values and resulting in a broad range of differential gene expression levels between control and infected samples (Additional file 1: Figure  S8 ). In contrast, microarrays typically report lower differential gene expression. To overcome this issue and compare a gene's level of expression across different datasets regardless of its differential expression range, we performed a rank product analysis, a non-parametric statistic used for detecting differentially expressed genes based on (log) fold changes.
The rank product analysis identifies genes that are consistently highly ranked in a number of datasets, and is a powerful approach to detect biologically relevant gene expression changes from heterogeneous datasets  [41]. For this we used the RankProd R package  [71] , which accepts pre-processed expression datasets produced from different platforms, and thus was appropriate for our ranked rsd. In short, all 7,077 genes were ordered based on their relative ranks and were given ranks from 1 to 7,077 in each transcriptome dataset. From these new rank values the rank product was computed for each gene across the 19 transcriptome datasets and p-values were assigned to genes using a permutation test (N = 10,000) to test for differential expression. We then corrected for multiple testing using the Benjamini-Hochberg procedure from the multtest R package  [72].
Using this approach, we identified genes that were (i) significantly up-regulated (genes with corrected rank product p-value < 0. 05), (ii) genes significantly downregulated (corrected p-value < 0. 05), and (iii) genes significantly differentially-regulated regardless of the orientation (up-regulated in some studies, down in the others, and here termed differentially-regulated; corrected p-value < 0. 05), across the 19 transcriptome datasets. For (i) and (ii) the ranked rsd was used to perform the analyses; note that a gene may be statistically significantly up-regulated across all 19 datasets by the rank product analysis but still be down-regulated in one or more datasets (and vice-versa for significantly downregulated genes). For (iii) we transformed the data as follows: the log 2 fold changes of each dataset were shifted by their median to obtain an equal amount of up-and down-regulated genes, to avoid biased analysis towards strongly up or down-regulated genes.
Signs of the log 2 fold changes were omitted and genes from group (iii) were first processed as described for the rsd, and then used to perform the analysis. Group (iii) therefore also contained genes of groups (i) and (ii). This resulted in three lists of genes and their ranks in the 19 transcriptome datasets, their rank product value, p-value as well as corrected p-value. To test whether a specific GO term was significantly over-and underrepresented in a specific gene set compared to a background set, we performed a two-sided Fisher's exact test using the stats R package  [69].
5  Methodological workflow of the directed rank product analysis (DiRank). This new method aims to identify genes with similar expression profile to a theoretical or observed profile of another gene. Gene expression values and profiles (geps) (shown in blue) and custom profile (cp) (shown in red), consisting of relative rank values, serve as input (yellow boxes). In rectangular matrices, gene expression values are reported in rows, while columns represent the transcriptome datasets. A custom profile can either be a user-defined profile or an existing gene expression profile. The directed rank product analysis aims to identify genes with a similar expression profile to the custom profile and to assign associated p-values. The custom profile is subtracted from each of the gene expression profiles and each difference (gep -cp) is transformed by 1 -| gep -cp|.
Specific response to pathogens To detect the specific response of honey bees to pathogens, we developed a novel method which applies the same rank product analysis as above, but on ranks that reflect the degree of similarity between a gene's expression profile and a custom expression profile across the 19 datasets. In other words, we identified genes whose expression profile is most similar to a theoretical expression profile corresponding to selected parameters. We called this method DiRank, for directed rank product analysis (R code available at https://idata. idiv. de/DDM/ Data/ShowData/35). For the current study, we identified genes that were specifically expressed in (i) abdominal tissues (i. e. gut, fat body or whole abdomens) after Nosema spp. infections (datasets #3-9, see Table  2 ) or (ii) in all tissues after RNA virus infection or Varroa parasitism, as mites transmit RNA viruses (datasets #12-19).
We thereby designed theoretical gene expression profiles across the 19 datasets. As an example, to identify genes specifically regulated in abdominal tissues after Nosema infection, we designed two custom profiles, one profile of increased expression in abdominal tissues of Nosema infected bees, with no changes of expression levels in other datasets, and a second profile of decreased expression in abdominal tissues of Nosema infected bees, with no changes of expression levels in other datasets. In the ranked rsd file, differential expression thresholds was defined, such as genes with only slight expression changes, with relative ranks ranging from 0. 3 to 0. 7, were considered as non-differentially expressed and assigned a relative rank of 0. 5. These thresholds reduced the impact of small variation in expression levels on the analysis.
All directed rank product analyses resulted in a list of genes annotated with their rank product values, their ranks according to the rank product value, p-value, corrected p-value, and ranks across the 19 transcriptome datasets. We used a cut-off uncorrected p-value < 0. 01 to identify genes specifically regulated under the selected parameters (Nosema in abdominal tissues or Varroa/virus) and to identify the associated overrepresented GO terms.
Gene co-expression networks Following an iterative process, we defined the expression profile of each gene as a custom profile for our DiRank method, so as to retrieve genes with similar expression across datasets, and the inverse of each gene profile (i. e. inverted) to identify genes with an opposite expression profile (see Additional file 1: Figure  S9 ). For each of the 7,077 genes from the ranked rsd file we obtained a list of other genes that showed a similar or an opposite expression profile. We then reconstructed the gene co-expression network using inter-gene connection falling under a Benjamini-Hochberg corrected p-value cut-off of 0. 05, after permutation test (N = 1,000). We visualized the network using Cytoscape  [73] , with genes as the 'nodes' of the network and gene interactions as the 'edges' between nodes, while a 'module' is a subset of interconnect nodes.
Additional files Additional file 1: Figure  S1 -S9. This file includes supplementary figures documenting our multidimensional scaling analysis results, a heat map of the differential expression of the 7,077 genes across the 19 datasets, a Venn diagram of differentially expressed genes, the expression profile of the gene coding for hymenoptaecin, the distribution of genes according to their number of inter-gene connections, the degree of connectivity of differentially expressed genes, the process of gene selection for this study, the distribution of genes' differential expression across datasets, and a diagram illustrating our new bioinformatics approach. (PDF 683 kb) Additional file 2: Table  ST1. List of 7,077 genes ordered by their rank product after the rank product analysis looking for up-regulated genes.
Funding This article is a joint effort of the working group TRANSBEE and an outcome of two workshops kindly supported by sDiv, the  Synthesis Centre for Biodiversity Sciences  within the  German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig , funded by the  German Science Foundation  ( FZT 118 ). New datasets were performed thanks to the  Insect Pollinators Initiative  ( IPI grant   BB/I000100/1  and  BB/I000151/1 ), with participation of the  UK-USA  exchange funded by the  BBSRC   BB/I025220/1  (datasets # 4 ,  11  and  14 ). The IPI is funded jointly by the  Biotechnology and Biological Sciences Research Council , the  Department for Environment, Food and Rural Affairs , the  Natural Environment Research Council , the  Scottish Government  and the  Wellcome Trust , under the  Living with Environmental Change Partnership.
Availability of data and materials Microarray expression data have been deposited on ArrayExpress (Archive of Functional Genomics Data from the European Bioinformatics Institute) under accession numbers E-MTAB-5508 (datasets #4, 11 and 14), and E-MTAB-5070 (#13). RNA-sequencing data (datasets #18 and 19) have been deposited to the Gene Expression Omnibus repository, under the accession number GSE94855.
Samples were representative of the mineralized, altered, and unaltered rocks from the Huincul Formation in the Tordillos area. Samples were analysed with a polarization-fluorescence Nikon Eclipse 50i microscope in order to detect organic compounds in diagenetic minerals. Modal mineralogy was determined by point counting (500 points) on each thin section of representative samples of unaltered and altered rocks. 2D image analyses were also used to obtain modal compositions and to quantify porosity using the software JMicroVision v. 1. 2. 7  (Roduit, 2008). Clay mineralogy was determined by X-ray diffraction (XRD), infrared spectroscopy (IR) and scanning electron microscopy (SEM/EDX) at the IC2MP Laboratory at the Université de Poitiers (France). Due to the poorly consolidated state of the samples, sandstones were gently disaggregated in water, and the clay fraction was taken into suspension by ultrasonic treatment.
Clay-size fractions (b2 and b 4 μm) were extracted by sedimentation for oriented and randomly oriented powder mounts. Only some problematic samples were treated with Ca saturation. All clay preparations were analyzed on a Bruker D8 Advance diffractometer. Diffracted beam radiation (CuKα1 + 2; 40 kV, 40 mA) was collected by a Lynx Eye detector. Relative humidity was not controlled during data acquisition. Quantitative clay analyses were carried out using the mineral intensity factor (MIF) method  (Moore and Reynolds, 1997). Infrared analyses were carried out on selected samples as a complementary technique for the determination of clay mineralogy. Middle-infrared (MIR) spectra (400 to 4000 cm -1 ) were acquired from KBr pellets using a Nicolet 760 (Fourier Transformation Infrared Spectroscopy -FTIR-) spectrometer equipped with a potassium bromide (KBr) beam splitter and DTGS-KBr detector.
Acknowledgements We are grateful to  Drs. Alex Brown  and  Franco Pirajno  for their constructive reviews, to  Ariel Testi  for his collaboration during field work, to  Dr. Laura Maydagán  for her collaboration in some of the sulfides analyses, to  Dra. Nora Cesaretti  for her contribution in the organic fluid inclusion analyses, to  Dr. Valentin Robin  due his support in the sample preparations for the XRD and MIRS analyses. This contribution forms part of a project financed by  CONICET  (PIP no. 6043  and PIP no. 1083 ), FONCyT, pict-2010-2608; Pict-prh  2008-00093 , Pict-2008-1120, Acmelab-IoStipend-2009,  PME-33 ,  Universidad Nacional del Comahue  (nos. 04/I167  and  04/I002 ) and a support from  European Union  through an  Erasmus Mundus mobility for EU and non-EU scholars-IMACS program  ( 2012-214  edition).
However, the common genetic variants explain only 3. 9% of the inherited component of EOC risk  15  and additional susceptibility loci are likely to exist, particularly for the less common, non-serous histotypes. We designed a custom Illumina array named the 'OncoArray', in order to identify new cancer susceptibility loci  18. The OncoArray includes ~533,000 variants (of which 260,660 formed a GWAS backbone) and has been used to genotype over 500,000 samples, including EOC case-control studies of the Ovarian Cancer Association Consortium (OCAC) and BRCA1 and BRCA2 mutation carriers of the Consortium of Investigators of Modifiers of BRCA1/2 (CIMBA). These data were combined with genotype data from the Collaborative Oncological Gene-environment Study (COGS) project  14, 19  and three EOC GWAS  8, 9.
Association analyses Genetic association analyses were performed using genotype data from 25,509 populationbased EOC cases and 40,941 controls from OCAC and meta-analysis of these data with 19,036 BRCA1 and 12,412 BRCA2 mutation carriers from CIMBA, of whom 2,933 and 954, respectively, were affected with EOC. The numbers of participants by study for OCAC and CIMBA are shown in Supplementary table 1 and Supplementary table 2, respectively. We used data from the 1000 Genomes Project  20  reference panel to impute genotypes for 11,403,952 common variants (MAF>1%) and evaluated the associations of these SNPs with EOC risk. In OCAC alone, nine histotypes were investigated (all invasive, serous invasive, HGSOC, LGSOC, serous borderline, LGSOC and serous borderline combined, ENOC, CCOC and MOC). Association analyses revealed six novel loci associated with serous EOC histotypes at genome-wide significance (p<5×10  -8  ): rs9870207 at 3q28, rs13113999 at 4q32.
SNPs were also intersected with ENCODE transcription factor ChIPseq data. All tissue types were evaluated for all risk loci. The SNP-biofeature intersections by tissue type are illustrated in Figure  2  and Supplementary table 6. Nine (1. 6%) of the 571 candidate causal SNPs lie in protein coding sequences. Five of these are synonymous and four are non-synonymous but predicted to be benign by Polyphen-2 (Supplementary table  6 ). Four SNPs lie within untranslated regions of protein-coding genes and so could affect mRNA stability: rs1051149 and rs1051150 in the 3' UTR of LAMA3 and rs12327412 in the 5' UTR of TTC39C, all at the 18q11. 2 locus; and rs1018128 in the 3' UTR of GMNC at 3q28. The majority of biofeature-SNP intersections (n=166, 29% of all candidate causal SNPs and 97% of candidate causal SNPs overlapping a biofeature) were for SNPs lying within active chromatin, and/or open chromatin.
Therefore we tested for enrichment of SNP-H3K27ac intersections at each locus because H3K27ac was the most comprehensively profiled regulatory feature across different tissue types (Supplementary table  7 ). At the 12q24. 31  ). At 10q24. 33 (which is associated with LGSOC and borderline SOC) we identified a significant enrichment of acetylated H3K27 in normal fallopian cells (FT33 P=1. 09×10 -4 , FT246 P=4. 29×10 -3 ), HGSOC ovarian cancer cells (UWB1. 289 P=6. 23×10 -3 ), MOC cells (GTFR230 P=5. 16×10 -3 ) as well as, somewhat surprisingly, colorectal cancer cells (HCT116 P=2. 64×10 -4 ) and cervical cancer cells (HeLa-S3 P=9. 60×10  -12  ). This locus contains several clusters of H3K27ac activity and TF binding in ovarian and ENCODE datasets, and these highly active regions showed extensive overlap with candidate causal alleles (Figure  3 ).
01. Assuming a polygenic variance of 1. 45  29  the 27 loci account for approximately 6. 4% of the polygenic risk in the population. Incorporating common EOC susceptibility variants into risk assessment tools will improve risk prediction and may be particularly useful for refining risk estimates in BRCA1 and BRCA2 mutation carriers. Some strata in the OCAC analyses pooled data from several studies from the same country. This might increase the potential for bias because of population stratification, but we expect any inflation due to population stratification to be effectively removed by adjusting for the principal components. In order to evaluate the possible magnitude of such a bias we compared the inflation of the median test statistic for the analysis of the OncoArray data stratified by study with an analysis in which all the cases and controls were combined into a single stratum.
The lack of heterogeneity between studies of varying designs, carried out in different populations, and the high levels of statistical significance, with confirmation of known EOC susceptibility loci, provide evidence that these are robust associations. Molecular analyses of genes and the tissue specific regulatory architecture at these loci, which combined publicly available datasets with systematic, large-scale genome wide profiling experiments, point to a small number of non-coding biofeatures and target genes that may play a histotype-specific role in EOC initiation and development. Detailed functional studies will be required to define the underlying biology of SNP-regulatory interactions to identify the causal SNP(s) at each locus, and to confirm which candidate susceptibility genes represent the targets of these risk SNPs.