dataesr/scientific-paragraphs-categorization huggingface dataset:
Gene expression dataset We applied our method to a cancer diagnosis problem extracted from microarray data. The used data comes from a study of cross-experiment compiling the gene expression profile from about 40,000 publicly available Affymetrix HG-U133Plus2 arrays  [24]. Combining different expression datasets gives a global gene expression map that contains variability related to the type of tissues and the experimental protocols. This allows us to address new questions and to make original studies that may lead to new biological discoveries. After quality control and normalization, the dataset contains the expression of 54675 probes from 27887 tissues of which 9450 are healthy and 18437 are cancer. The dataset is accessible via the ArrayExpress database (accession number E-MTAB-3732). We divided this dataset into a training set (14750 cancer, 7560 non-cancer) and a test set (3687 cancer,1890 non-cancer).
3  Comparison of LRP and WM on a toy example values since the latter is independent from the input data. By definition WM values are computed only from the weights of the connections. These values could represent the importance of the inputs if the variables were independent and uniformly distributed. However, this is not the case in the context of gene expression data. The interpretation must depend on the distribution of the data that are provided to the network. We propose an experiment to compare the quality of biological interpretation produced by WM and the LRP on our gene expression dataset. In this experiment, we measure the consistency of the LRP and WM interpretations of our cancer prediction network with known biological knowledge. We use all the test examples predicted as cancer by our network to compute the average of the LRP relevance scores of the input neurons.
Discussion We point out that the goal of the interpretation is to explain how the model works and not how biology works. Sometimes, there is no obvious relation between the biological functions or metabolic pathways, returned by the interpretation, and the predicted phenotype. This does not necessarily mean that the predictions are not reliable. Recall that a model looks for correlations between the output and the input and not for causalities. When a biological function, that is not related to the phenotype, is identified in the important neurons, it is possible that this function either has an indirect correlation or is linked by an unknown causality relation with the phenotype. Moreover, let us not forget that the biological databases used are not a perfect description of biology, but just a representation of current knowledge of biology.
In these databases, a lot of information is missing because it is unknown and some available information may be incorrect. We finish by giving the different types of conclusions that we can draw from a biological interpretation of a model. We identify three cases based on the results of the interpretation. In the first case, the majority of the elements provided by the biological interpretation are related to the predicted phenotype. This means that the model bases its predictions on elements consistent with the biological knowledge. This should improve the trust in the model in addition to its prediction performance. The second case is the opposite. Most parts of the elements provided by the biological interpretation are known to be unrelated to the predicted phenotype. Since the predictions are based on elements inconsistent with the current biological knowledge, the reliability of the model must be questioned.
Conclusion In this paper, we propose an original approach for biological interpretation of deep learning models for phenotype prediction from gene expression data. Our main objective is to identify the neurons and inputs of the NN that contribute to the predictions and to link them to biological knowledge. The model is reduced to a sub-network containing the relevant connections and neurons involved in the prediction. These neurons are then associated with a list of genes and the corresponding biological knowledge (GO, KEGG, and DOLite). Our experiments, based on cancer prediction, show that (1) deep learning approach outperforms classic machine learning methods on large training datasets, (2) our approach produces interpretations more coherent with biology than the state-of-the-art based on WM approach, (3) we can provide a comprehensive explanation of the predictions for biologists and physicians.
Methods We present the deep neural network architecture used for gene expression data and our biological interpretation approach. The gradient method for neural network interpretation is the Layer-wise Relevance Propagation (LRP), which is adapted to identify the most important neurons that lead to the prediction as well as the identification of the set of genes that activate these important neurons. Finally, the important neurons and genes are linked to the Gene Ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Disease Ontology Annotation List (DOLite) in order to propose a biological interpretation of the neural network model.
size. We show that for small training set size ( n < 1000 ) the state-of-the-art algorithms provide higher accuracy than NN, for moderate training set size ( 2000 < n < 10000 ), the accuracy of NN is similar to the accuracy of the other methods, for large training set size ( n > 15000 ) the NN significantly outperforms the other methods. When the whole training dataset is used, the accuracy of the NN is estimated at 95. 7% on the test set. Among the 238 prediction errors, 115 are non-cancer patients predicted as cancer and 123 are cancer patients predicted as non-cancer. This model will be analyzed in detail in the rest of the paper. We note that the performances of SVM, LASSO, RF and XGB reach their maximum for n ∈ [10000, 15000] examples, adding new training examples does not improve the accuracy.
2a, b )  21. Furthermore, the two polyadenylation (PolyA) sites identified by 3′ RACE agree with Cancer Genome PolyA Site & Usage data (Supplementary Fig. 2a, b ). A closer examination of lincNORS revealed well-recognized features of a bona fide lncRNA. First, multiple coding-potential analyses produced consistently low scores (Supplementary Fig. 2c ), and the SORFS mass spectrometry database (http://sorfs. org/ database) did not list any peptides associated with this locus at the time of writing. Second, the locus is separated from the flanking ones. Third, as expected for a lncRNA, lincNORS is expressed in a tissue-specific pattern, with a higher abundance in hormone-responsive tissues and skeletal muscle (Supplementary Fig. 3 ). Furthermore, GTEx data sets indicate an overall abundance comparable to, or higher than, the well-established lncRNAs HOTAIR, CCAT1, and MIR210HG (Supplementary Fig.
2c ). This result is consistent with results reported by Zhang et al. listing MIR193BHG (lincNORS) among the HIF-2 targets  23. The preferential role of HIF-2 has been previously described for hypoxia lncRNAs, in particular NEAT1, MALAT1, and NLUCAT1  6, 11. A survey of the genomic region surrounding lincNORS using available ChIP-seq data sets (GSE28352) did not reveal functional HIF recruitment sites within the lincNORS proximal regulatory region, in contrast to standard hypoxia responders typically driven by HIF-1 from proximal promoters (e. g. , CA9 and miR210/MIR210HG) (Fig. 2d )  24. Based on the knowledge that HIFs (and in particular HIF-2) sometimes transactivate from distant enhancers, we identified major peaks at the neighboring locus LINC02130 (Fig. 2e  and Supplementary Fig. 5a )  24, 25.
This locus is a plausible distal regulatory element based on (1) ENCODE genome segmentation features (Supplementary Fig. 5b ); (2) chromatin interaction with the lincNORS proximal regulatory region (RNAPII ChIA-PET GSE33664) (Fig. 2e )  26  ; (3) LINC02130 transcript is much less abundant than lincNORS in tissues and cell lines (Fig. 2e -GTEx track summary and Supplementary Fig. 1b  for MCF7 cells); and (4) evidence of activation by hypoxia but the corresponding transcript has very low abundance (Supplementary Fig. 1b ). Interestingly, the HIF prolyl-hydroxylase inhibitor DMOG was much less effective in inducing lincNORS than low O 2 (Supplementary Fig. S4c, d ). Despite inducing the expected hypoxic responses (increased HIF-2 protein abundance; CA9 transcript upregulation), DMOG impact on lincNORS was modest in MDA-MB-231 and reproducibly absent in MCF7.
7b ). We next performed RNA-seq analysis in hypoxic MCF7 cells to determine the consequences of lincNORS downregulation (Supplementary Fig. 8a  and Supplementary Data 4). Differentially regulated transcripts were subjected to Gene Set Enrichment Analysis (GSEA) and Ingenuity Pathway Analysis (IPA), both of which identified the upregulation of cholesterol biosynthesis as the top transcriptomic signature (Fig. 4a  and Supplementary Fig. 8b ). As illustrated in Fig. 4b  and Supplementary Fig. 8c , the program responsible for the de novo synthesis of cholesterol and related steroids was coordinately upregulated by lincNORS siRNA, with significant upregulation of many individual genes. A similar effect was evidenced in ER-negative MDA-MB-468 cells using the AmpliSeq platform (Fig. 4c, d  and Supplementary Data 5).
8d ). We then performed gain-of-function studies and we stably expressed full-length lincNORS (NORS_OE) or vector-only control and performed RNA-seq on cells growing under ambient O 2 (Supplementary Fig. 9a  and Supplementary Data 6). Both GSEA and IPA identified cholesterol biosynthesis among the top suppressed signatures (Supplementary Fig. 9b, c ), with multiple components downregulated (Supplementary Fig. 9d ). We confirmed by qPCR the downregulation of MSMO1, SQLE, and HMGCR as a result of lincNORS transcript induction (Supplementary Fig. 9e ). As expected, this intervention did not significantly affect miR-193b or miR-365a levels (Supplementary Fig. 7b ). Based on RNA-seq data, lincNORS downregulation and overexpression did not significantly affect the expression of flanking genes MKL2/MRTFB and PARN, thus indicating a predominantly trans activity (Supplementary Fig.
Furthermore, an examination of the MKL2-PARN genomic region immediately revealed that these SNPs are more likely to impact MIR193BHG/ lincNORS than MKL2 (or PARN). First, the top GWAS-reported SNP for puberty onset (rs246185) is very close to the TSS of lincNORS and overlaps with histone acetylation/methylation peaks (Fig. 7a ). Second, expression quantitative trait loci (eQTL) analysis in the GTEx database shows that these genetic variations are correlated with lincNORS, but not MKL2 expression in several tissues (example shown in Fig. 7b ). Colocalization analysis showed that the GWAS signal for age at menarche and the eQTL have a high probability of representing the same underlying signal in skeletal muscle and artery (Supplementary Data 8)  52, 53. Finally, based on ARIES mQTL information, these variants also exhibit significant associations with the methylation of the MIR193BHG locus rather than MKL2.
A schematic diagram of the cholesterol synthesis pathway is shown on the right. Differentially expressed genes (FDR < 0. 05) are marked with red. c Gene set enrichment analysis of hallmark gene sets associated with lincNORS knockdown versus the control in MDA-MB-468 (FDR < 0. 1). Enrichment plot of cholesterol homeostasis gene set is shown on the right. d Heatmap depicting expression of genes in the cholesterol synthesis pathway in MDA-MB-468 cells transfected with lincNORS siRNA or siRNA control. The transcriptome was profiled by AmpliSeq. e qPCR analysis of selected cholesterol synthesis genes in MCF-7 cells transfected with two different lincNORS siRNA. Data represent mean ± SD from three/four biological replicates (*p < 0. 05, **p < 0. 01, two-way ANOVA with Dunnett's test). f qPCR analysis of selected cholesterol synthesis genes in PC-3 cells transfected with two different lincNORS siRNA.
Data represent mean ± SD from three biological replicates (*p < 0. 05, **p < 0. 01, ***p < 0. 001, twosided Student's t-test). rs181766, which is in high linkage disequilibrium with rs246185 and exhibits similarly significant phenotypic associations according to the UK Biobank database (Fig. 8a ). However, based on ChIP-seq and open chromatin by Formaldehyde Assisted Isolation of Regulatory Element (FAIRE-seq) data sets from ENCODE, rs181766 represents a potentially stronger causal effector SNP candidate (Supplementary Fig. 13b ). Sequences corresponding to the reference or alternative alleles of rs246185 (T/C) and rs181766 (T/C) were subcloned into pGL4 vectors and reporter activity was assayed in MCF7 cells. Although the Alt rs246185 (C) variant had no measurable impact, substitution of the rs181766-C allele for the T allele had a significant effect on reporter activity (Fig.
20  following the method formulated by Maud H. W. Starmans et al. 69. HIF ChIP-Seq data (GSE28352) were used to assess whether HIF was directly involved in lincNORS induction by hypoxia. Raw data of HIF-1α ChIP-Seq, HIF-2α ChIP-Seq, and pre-immune ChIP-seq control were downloaded and aligned to hg19 using Bowtie (version 1. Subcellular fractionation. Subcellular fractionation was performed as described  70. Briefly, cells were pelleted and resuspended in ice-cold HLB buffer (10 mM Tris pH 7. 5, 10 mM NaCl, 3 mM MgCl 2 , 0. 3% (v/v) NP-40, 10% (v/v) glycerol) supplemented with either RiboLock RNase Inhibitor (Thermo #EO0381) or Halt Protease Inhibitor (Thermo #78442) for RNA and protein samples respectively. After 10 min incubation on ice, samples were centrifuged at 1000 × g at 4 °C for 3 min and the supernatant was collected as the cytoplasmic fraction.
Reverse transcription was performed using TaqMan® Reverse Transcription Reagent (Applied Biosystems, CA) from 1 µg of RNA, and qPCR was performed using gene-specific primers (Supplementary Data 9) and SYBR green master mix (Applied Biosystem, CA). 18 S rRNA or B2M RNA level was used for normalization. For microRNA analysis, reverse transcription was performed with TaqMan® MicroRNA Reverse Transcription Kit using 10 ng RNA, and qPCR was performed using TaqMan® Gene  AmpliSeq and data analysis. Automated library preparation was performed according to the protocol available with the Ion Chef (MAN0013432) and using the 14,240,000 14,250,000 14,260,000 14,270,000 14,280,000 14,290,000 14,300,000 14,310,000 14,320,000 14,330,000 14,340,000 14,350,000 14,360,000 14,370,000 14,380,000 Ion AmpliSeq Kit for Chef DL8 (A29024).
Library was further diluted to 50 pM for sequencing. Template preparation sequencing were executed using the Ion PI Hi-Q Chef Kit (A27198, ThermoFisher Scientific) with template quantification using the IonSphere Quality Control Kit (4468656, ThermoFisher Scientific). The raw read counts are included in the Supplementary Data 5. Standard workflow for gene-level expression was conducted using Torrent Suite Version 4. 4. The AmpliSeq Transcriptome Plugin was used to generate standard expression files. Differential expression was performed in R using the DESeq2 package (version 1. 12. 3) in R/Bioconductor (R version 3. 3. 1). Differential expression with FDR < 0. 05 was considered significant. Pathway analysis was performed using IPA and GSEA. Gene set enrichment analysis. Gene set enrichment analysis was performed following GSEA User Guide (The Broad Institute)  74.
Briefly, each condition was considered as a group and gene list was ranked with GSEA default ranking metrics. Gene sets from the Molecular Signature Database were used in the analysis to identify the pathways significantly enriched in each group. Gene sets were permutated 1000 times to obtain empirical FDR corrected p-values. Western blot analysis. Cell lysates were prepared in RIPA buffer containing protease inhibitor cocktail (ThermoFisher Scientific). For HIF western blot, cells were fractionated using RLN buffer (10 mM Tris-HCl pH 8. 0, 140 mM NaCl, 1. 5 mM MgCl 2 , 0. 5% NP-40, proteinase inhibitor cocktail) and the nuclear fraction was collected for analysis. Proteins were separated by SDS-polyacrylamide electrophoresis (SDS-PAGE), then transferred onto PVDF membranes (Biorad). Primary antibodies used in this study included those that recognized HIF-1α (R&D AF1935), HIF-2α (R&D AF2886; Cell Signaling #7096), β-actin (Thermo MA5-15739; Santa Cruz sc-1615), β-tubulin (Santa Cruz sc-5274), acetyl-histone H3 (Millipore 06-599), Lamin-B1 (Abcam ab16048), SREBP2 (R&D AF7119), hnRNPC (Santa Cruz sc-32308), TRIM28 (Abcam ab10483), RALY (Abcam ab170105), and Enterobacterio Phage MS2 coat protein (Millipore ABE76).
For a complete list of antibodies used please see Supplementary Data 10. Total cholesterol quantification. Cellular cholesterol was measured using the Wako Cholesterol assay kit (WAKO, #439-17501) and liquid chromatographymass spectrometry. For the measurement using Wako Cholesterol assay kit, cellular lipids were extracted using 750 µL chloroform: methanol (2:1) solution. After extraction for 30 min on a rocker, samples were centrifuged for 10 min at 1500 × g at 4 °C, then the lower chloroform phase was collected and air-dried in a fume hood until all liquid evaporated. The lipids were resuspended in 10% Triton-X100 in isopropanol. Total cholesterol content was measured using the Wako Cholesterol assay kit (WAKO, #439-17501) according to the manufacturer's instructions. For the measurement using liquid chromatography-mass spectrometry, Oxysterols were extracted from cells using a modified sterol extraction method from McDonald et al.
5 mL of tetrahydrofuran (THF). 0. 15 mL of this mixture were added to each sample along with 0. 02 mL of triethylamine (TEA). The reaction proceeded at room temperature for 1 h, then the samples were dried and subsequently reconstituted in 0. 1 mL of acetonitrile just prior to LC/MS/MS analysis. The analysis was done with an Agilent 1200 liquid chromatography system coupled to an Agilent 6460 QQQ mass spectrometer (Santa Clara, CA) in Bindley Bioscience Center at Purdue University. Data were processed using Agilent Masshunter Quantitative analysis software (V. B. 06). Bromouridine (BrU)-labeled RNA Immunoprecipitation and Bru-Seq. BrU labeling of newly synthesized RNA and BrU-RNA isolation were performed following the method established by Ljungman's group  22    transfected with either control siRNA or lincNORS siRNA were cultured in hypoxic conditions.
Purified RNA was then converted to cDNA using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific), following the manufacturer's instructions. The qPCRBIO SyGreen Mix (PCR Biosystems) was used as a master mix for qRT-PCR. The reactions were performed on a CFX384 thermal cycler (Biorad) in technical duplicates for each target and for a total of four independent biological replicas. RALY RIP-seq data were previously published by Macchi lab and re-analyzed for lncRNA and cholesterol gene enrichment using the same method as previously described  44. Luciferase assay. Oligos of lincNORS promoter region containing rs246185 and rs181766 reference or alternative allele were ordered as gBlocks fragments from IDT and inserted into pGL4 vectors. MCF7 were seeded at 70-80% confluence in 24-well tissue culture plates 24 h before transfection.
53. Colocalization analysis was performed in R (version 3. 4. 2) using the "coloc" package  52. Data are presented as the posterior probability of the two signals colocalizing at this locus (PP4) and the posterior probability of the two signals colocalizing conditioned on there being a signal from both the eQTL and the GWAS data (PP4/PP3 + PP4). A higher value indicates a higher likelihood that the signal is shared. Statistical analysis. Statistical analysis was performed in Excel, R (version 3. 0. 1, http:///www. rproject. org/) or GraphPad Prism (version 6). Statistical significance of differences was determined by unpaired two-sided Student's t-test or ANOVA with post hoc comparison where appropriate unless stated otherwise in the figure legend. The difference was considered statistically significant when p-value was <0.
2. 1. 1). The peaks were visualized in UCSC genome browser. The other ChIP-Seq data presented in this study were obtained from ChIP-Altas (http://chip-atlas. org/) and visualized in IGV. Longrange chromatin interaction in MCF-7 was determined by ChIA-PET data published by Li et al. (GSE33664) and visualized in WashU EpiGenome Browser 26.
Acknowledgements We thank members of the  IU Simon Cancer Center In Vivo therapeutics (IVT) Core  for assistance in animal experiments. We also thank  Dr. Je-Hyun Yoon  for the technical advice on the MS2 RNA pulldown system, and  Dr. Mervin Yoder  for the critical comments and suggestions for the manuscript. This work was supported in part by  NIH   R01CA155332 ,  Pancreatic Cancer Signature Center at Indiana University, and Biomedical Research Grant at Indiana University. M. G. was supported by the  NIA IRP, NIH. C. M. N. was supported by the  Foundation for Cellular and Molecular Medicine. D. L. C. was supported by  1K99HD099330-01. J. R. is supported by  Genome Canada Genome Technology Platform  Grant  12505 ,  Canada Foundation for Innovation  Project  33408.
Variations in COVID-19 transmission rates across countries We assembled a composite dataset by combining a dated phylogeny of SARS-CoV-2 (Fig. 2a ), publicly available from Nextstrain  20  and built from the GISAID initiative data  21  , with a detailed timeline of non-pharmaceutical interventions available from the Oxford COVID-19 Government Response Tracker (OxCGRT)  22. Extended Data Fig. 1  shows a flowchart outlining the data sources, sample sizes and selection steps of the study. Phylogenetic and intervention data covered the early phase of the epidemic up to May 4, 2020. The 5,198 SARS-CoV-2 genomes used to reconstruct the dated phylogeny were collected from 74 countries. Detailed per-country data including sample sizes are shown in Extended Data Table  1. Among the 10,394 branches in the phylogeny, 2,162 branches (20.
8%) could not be assigned to a country with >95% confidence and were excluded, also reducing the number of represented countries from 74 to 59 (Extended Data Fig. 1 ; a comparison of included and excluded branches is shown in Extended Data Fig. 2 ). The remaining 4,025 internal branches had a mean time-to-event (delay between transmission events) of 4. 4 days (Fig. 2b ). These data were congruent with previous estimates of the mean serial interval of COVID-19 ranging from 3. 1 days to 7. 5 days  23. The 4,207 terminal branches had a mean time-to-censoring (delay from infection to detection) of 10. 6 days (Fig. 2a, b ). This pattern of longer terminal vs. internal branches is typical of a viral population in fast expansion  11. We compared the timing and dynamics of COVID-19 spread in countries represented in our dataset (Fig.
SARS-CoV-2 genome sequences have been continuously submitted to the Global Initiative on Sharing All Influenza Data (GISAID) by laboratories worldwide  21. To circumvent the computational limits of phylogeny reconstruction and time calibration techniques, the sequences of the GISAID database are subsampled before analysis by the Nextstrain initiative, using a balanced subsampling scheme through time and space  20, 43. Phylogenetic reconstruction uses maximum-likelihood phylogenetic inference based on IQ-TREE  44  and time-calibration uses TreeTime  45. See  46  for further details on the Nextstrain bioinformatics pipeline. A dated phylogeny of 5,211 SARS-CoV-2 genomes, along with sampling dates and locations, was retrieved from nextstrain. org/ncov on May 12, 2020. Genomes of non-human origin (n = 13) were discarded from analysis.
Acknowledgements: We thank  Philip Supply ,  François Vandenesch ,  Jean-Sébastien Casalegno ,  Vanessa Escuret , and  Christophe Ramière  for fruitful discussions and reviews of our work. We thank the  GISAID, Nextstrain and OxCGRT teams  for making their high-quality datasets available to the community. A list of authors and laboratories contributing SARS-CoV-2 genome sequences is shown in Extended Data Table 3.