dataesr/scientific-paragraphs-categorization huggingface dataset:
Quantitative digital image analysis of immunohistochemical staining Quantitative comparison of immunohistochemical stainings were measured by the ImageJ digital image analysis public domain software (W. S. , Rasband, Image J, U. S. National Institutes of Health, Bethesda, MD; http://rsb. info. nih. gov/ij/, 1997-2011). The IHC profiler software plugin was required. The immunopositivity was expressed as a percentage of the total software-classified areas and the data obtained were plotted as histograms.
This study was supported by the  University of L'Aquila, L'Aquila, Italy ; the  Inserm   U995 , Lille, France;  Nogra Pharma Ltd, Dublin, Ireland ; the  Intestinal Biotech Development of Lille, France ; the  Digest Science foundation  and the association  F. Aupetit, Lille, SpA, Milano, Italy , has received research grants from  Giuliani Italy, Alfa Wassermann Italy and Sofar Italy. The authors also received funding from  Nogra Pharma Ltd  and the  Intestinal Biotech Development of Lille, France.
Another possibility is that these evolution studies have mainly focused on the overall level of mtDNA mutations, but not on their level in specific regions or individual positions. In contrast, functional studies have mostly focused on individual mutations. To shed more light on the subject, we aimed to explore whether site-specific natural selection might affect tumor somatic mtDNA mutations and, more importantly, their functions. We performed a systematic and comprehensive exploration of mtDNA mutations in various cancer types, and we developed the largest collection of mtDNA mutations in cancer to date using next generation sequencing (NGS) data generated in various sources, including The Cancer Genome Atlas (TCGA), International Cancer Genome Consortium (ICGC), other publications, as well as our own lab.
Data of somatic tumor mtDNA mutations from public resources and our lab We identified 5423 somatic mtDNA mutations from whole genome sequencing (WGS) and whole exome sequencing (WES) data reported in literature  [17, 18, 21] , as well as in TCGA (https://portal. gdc. cancer. gov/exploration?searchTableTab=mutations&ssmsTable_size=100) and ICGC (https://dcc. icgc. org/repositories) datasets. For patients who were present in two or more datasets and sequenced by the same method (WES or WGS), the mutations were intersected; for patients who were present in two or more datasets and sequenced by different methods, the data from WGS were retained. Seven patients with an extremely large number of somatic mutations (≥13) were excluded from further analyses. In addition, we used mtDNA capture sequencing approach to identify 104 somatic mutations from 54 colon cancer patients, 357 somatic mutations from 143 hepatocellular carcinoma (HCC) patients, and 36 somatic mutations from 15 kidney cancer patients.
All these mutations derived from tumor tissue in our lab were verified by parallel sequencing of the non-tumor tissue of the same patient. Considering the error rate of the Illumina sequencing platform, only the somatic mutations with ≥2% of variant allele fractions were counted. Finally, we generated a combined dataset containing 26 tumor types, 3277 patients, and 5920 mutations, the most comprehensive tumor mtDNA mutation dataset to date (Additional file 1: Table  S1  and Additional file 2: Table  S2 ). Among them, there were 2465 mutation records with heteroplasmic level data from three publications. The heteroplasmic levels of 497 mutations from our lab data were calculated as the ratio of the mutant reads in all reads that cover the mutation positions. All mutations were annotated using Annovar software  [22].
Germline variation data in natural human population Germline variation data containing 4322 mtDNA single base substitutions in the natural human population was extracted from phylotree (Build 17, 18 Feb 2016)  [23]. The htm file of mtDNA tree Build 17 containing all mtDNA variations among populations to date was downloaded (http:// www. phylotree. org/builds/mtDNA_tree_Build_17. zip) and the base substitution information was extracted using a homemade python script. (Additional file 3: Table  S3 ).
Simulation of mtDNA mutations We now describe a process to simulate mtDNA mutations under neutrality as a reference. Since it is impossible to know the mutation rate of each base under neutrality, we assumed that the vast majority of the mtDNA mutations in tumors were derived under neutrality, and selection plays a negligible role in shaping the overall mutation pattern. Under this assumption, we utilized the mutations in the large-scale dataset we collected to calibrate neutral mutation rates. Specially, we calculated mutation frequency of each of the possible nucleotide base substitutions (12 types in total) based on the observed base substitutions in this dataset. That is, for each of the 12 base substitutions, we counted the number of occurrences observed in this dataset, to represent the mutant frequency under neutrality; such estimates are expected to be robust given the large size of the collected mutation dataset.
In simulating a mutation dataset under neutrality, we first stored all the calibrated base substitution frequencies in a file (Additional file 4: Table  S4 ), and then used a homemade python script to randomly simulate 12 types of base substitutions, according to the mutation frequencies stored in this file (Fig. 1b ). The data thus simulated have the same distribution of mutations observed in the large-scale mutation dataset we collected, representing mutations derived under neutrality. We repeated this process N times to generate random replicates of mutations, and used the simulated datasets as a reference of neutrality to compare against the observed mutations in a dataset. Simulating scripts are available from https://github. com/gadbee/Bioinformatics.
Statistical analysis All statistical tests were performed using Graphpad Prism software (GraphPad Software, Inc. ). Comparisons between groups with continuous variables and categorical variables were performed using Student's t test, and then ANOVA and Chi Square test. Pearson's correlation coefficient was used to evaluate the correlations among related variables. All p-values were two-tailed. Residual was calculated by the difference between simulated data and observed data of a mutant codon or an amino acid. Standardized residual was calculated by formula: δ Ã ¼ δ-δ À Á =σ, in which δ* stands for standardized residual, δ for residual of a certain data point, δ for the average of all residuals, and σ for the standard deviation of all residuals. Standardized residuals greater than 2 were considered significant deviation of observed values from simulated values.
Characteristics of 5920 Single Base mtDNA mutations A total of 5920 single base mutations from 3277 patients of 26 tumor types (Additional file 1: Figure  S1 ) were collected from two databases (TCGA and ICGC), three previous publications, and our mtDNA sequencing data (see Materials and Methods and Additional file 2: Table  S2 ). Similar to previous report  [17] , a high consistency among WGS, WES, and mtDNA capture data was observed by comparing the proportion of synonymous and nonsynonymous mutations (Additional file 1: Figure  S1B ). In addition, we also analyzed the correlation of base substitution numbers or the amino acid substitution number between WGS and WES data (Additional file 1: Figure  S1C, D ) or between WGS + WES and mtDNA capture data (Additional file 1: Figure  S1E, F ) and found a high consistency of substitution pattern among the data derived from different sequencing methods.
1d ). The average mutation number per patient in different tumor types was 1. 82, ranging from 2. 29 (bile duct cancer) to 1. 22 (central nervous system tumors) (Fig. 1e , left). Interestingly, the average mutation number in sarcomas was significantly lower than those in carcinomas (Fig. 1e, right ). The mutation levels, available for 2962 mutations, ranged from 2% heteroplasmy to 100% homoplasmy. Of the 5920 somatic mutations, more than half (2951, 51. 10%) were not registered in the mitomap database updated on 29 Dec 2016  [25].
The relaxation of negative selection on tumor somatic mtDNA Single Base mutations To evaluate the effects of natural selection on single base mutations, mutation patterns were analyzed in general population, tumor, and simulation data. Compared with general population data, tumor data showed a higher proportion of nonsynonymous mutations in mRNA-coding regions, more deleterious mutations in tRNA coding regions, higher impact on protein function, and more codons that cannot be completely complementary with anticodon (Fig. 2a-d , left panel). These phenomena were considered as a result of positive selection in previous studies  [17]. Alternatively, another plausible model has been suggested that these tumor mutation patterns were shaped by mutation rate bias and relaxation of negative selection  [13, 17].
To further investigate whether the observed mutation patterns in the datasets that we developed are consistent with the relaxation of selection, we generated simulation data of mtDNA mutations under neutrality based on the mutation rate bias shown in Fig. 1b  (see Methods). No significant differences existed between the observed mutations and the simulated data, supporting the effects of the relaxation of selection on mtDNA mutations at the genome scale (left panel of Fig. 2 a-d ). To rule out the possibility that mutations with different heteroplasmic levels are under different selective constraints, we stratified the analyses by heteroplasmic levels (≤0. 05, 0. 05-0. 9, ≥0. 9) and found that the same conclusion holds regardless of the heteroplasmy levels (Fig. 2 a-d , middle panel). Furthermore, we found that heteroplasmic levels were not significantly different between synonymous and nonsynonymous mutations, between deleterious and benign mutations in tRNA coding regions, among mutations with different damage grades in mRNA coding regions, and among mutations resulting in different codon substitutions (Fig.
Discussion Somatic mtDNA mutations have attracted significant attention because of their high frequencies in cancers  [3, 4]. Specific mtDNA mutations have been implicated in tumor development and metastasis  [5, 6]. However, an increasing number of recent studies have shown that tumor mtDNA mutations are under relaxed negative selection  [13, 17, 18]. While the differing findings of these studies has raised controversy over the existence of specific functional mtDNA mutations, the indications of the relaxation of negative selection on the majority of tumor mtDNA mutations does not completely exclude the possibility of site-specific selection on individual mutations. To resolve these controversies, we developed a comprehensive dataset of tumor mtDNA mutations and systemically analyzed the effects of site-specific selection on these mutations.
Genotype-first combined with reverse phenotyping has shown to be a powerful tool in human genetics, especially in the era of next generation sequencing. This combines the identification of individuals with mutations in the same gene and linking these to consistent (endo)phenotypes to establish disease causality. We have performed a MIP (molecular inversion probe)based targeted re-sequencing study in 3,275 individuals with intellectual disability (ID) to facilitate a genotype-first approach for 24 genes previously implicated in ID. Combining our data with data from a publicly available database, we confirmed 11 of these 24 genes to be relevant for ID. Amongst these, PHIP was shown to have an enrichment of disruptive mutations in the individuals with ID (5 out of 3,275). Through international collaboration, we identified a total of 23 individuals with PHIP mutations and elucidated the associated phenotype.
Materials and methods We performed targeted re-sequencing of the coding sequence of 24 candidate ID genes (Table  1 ) as described previously  [25, 26]. A set of 1,490 MIPs was designed and DNA from 3,275 individuals with ID was collected through an international collaboration (Nijmegen, The Netherlands; Troina, Italy, Seattle, United States). MIP libraries of a maximum of 192 samples each were pooled and sequenced on one lane of a HiSeq2000 instrument (Illumina), resulting in an average coverage of >200 fold per individual. Sequence reads were mapped to UCSC human reference genome assembly hg19 and variants were called as described previously  [25]. Variant annotation was performed with an in-house pipeline  [3] , providing information of the variant effect at cDNA and protein level, evolutionary conservation scores, and frequency information based on > 15,000 in-house exomes and public databases such as ExAC and GoNL  [32] [33] [34].
*ExAC coverage data are provided in Supplementary Fig. S4. Most individuals had mild-severe ID whereas the individuals who did not have ID, had speech problems, global developmental delay in early childhood or learning difficulties. Behavioral problems were observed in 18 individuals (78%) and consisted of ADHD, features of autism, problems with impulse control, aggressive behavior, mood disorder, and anxiety. In addition to the DD/ID, the other most striking and overlapping feature in these individuals was overweight, which was observed in 17 individuals (74%), of whom 9 were obese. Shared dysmorphic facial features comprised high forehead (67%), full eyebrows/synophrys (59%) with 15 individuals having a characteristic upturned nose with thick alae nasi (68%), a long philtrum (45%) and thin lips (36%) (Fig.
Acknowledgements We are grateful to all individuals and their parents for participating in this study. We are also grateful to the  Radboud Genomics Technology Center  for performing Sanger sequencing. We also want to thank  Choli Lee  for technical assistance with HiSeq sequencing. We thank  Boris Keren ,  AP-HP, Groupe Hospitalier Pitié-Salpêtrière, Département de Génétique, Paris, France , for interpretation of the exome data of Individual 15. We thank  Ineke van der Burgt ,  Human Genetics ,  Radboudumc ,  Nijmegen ,  The Netherlands  and  Christina Ringmann Fagerberg ,  Department of Clinical Genetics, Odense University Hospital, Denmark  for clinical support. We are grateful to all of the families at the participating Simons Simplex Collection (SSC) sites, as well as the principal investigators ( A.
A global transcriptomic analysis by RNA sequencing (RNA-seq) was conducted on roots collected 3 and 24 h posttreatment (hpt). The results showed that for both rootstocks, the N-related gene expression pattern was extensively regulated in response to a local N supply. Interestingly, CS/RGM presented many more differentially expressed genes (DEGs) than CS/1103P according to the N treatment, and a weighted gene co-expression network analysis (WGCNA) approach allowed the identification of N-related modules specific to the CS/RGM combination.
RNA extraction, library preparation, and sequencing Total RNAs were isolated from 30 samples (three biological replicates per condition; the conditions are illustrated in Fig. 1C ), as described by  Cookson et al. (2013). Total RNAs were sent to IGBMC Microarray and Sequencing Platform (Strasbourg, France) and their quality was further determined using an Agilent 2100 Bioanalyzer (Agilent Technologies). The preparation of the libraries was performed by the IGBMC Microarray and Sequencing Platform using the TruSeq RNA Sample Prep Kit (Illumina Inc. ). Libraries were sequenced with an Illumina HiSeq 2000 sequencer (Illumina Inc. ) and 50 bp single reads were generated. The data have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (http://www. ncbi.
Pre-processing of RNA-seq data Quality checks were done with fastQC on the Galaxy website hosted on the Abims Bioinformatics Platform (Roscoff, France) (http:// abims. sb-roscoff. fr/resources/galaxy, last accessed 19 June 2017). The raw reads were then filtered by the Genotoul Bioinformatics Platform (Toulouse, France) with trimming of the low-quality reads with an average Phred score <28 (see Supplementary Table  S1  at JXB online).
De novo assembly De novo assembly was performed by the Genotoul Bioinformatics Platform using DRAP v1. 3  (Cabau et al. , 2017). This meta-assembler first performed an assembly with Trinity  (Grabherr et al. , 2011)  or Oases  (Schulz et al. , 2012)  using the multi-kmers approach. After assembly compaction, poorly supported and misassembled transcripts were filtered out using a reads per kilobase per million mapped reads (RPKM) threshold of 1, 3, or 10 (Supplementary Table  S2 ). Considering the results of the different assemblies (Supplementary Table  S2 ), the assembly with Oases multi-kmers and an RPKM threshold of 1 was selected.
Functional annotation of de novo assembled transcriptome Functional annotation was performed by the Genotoul Bioinformatics Platform with series of BLASTX against different protein databases from Ensembl hosted on the Genotoul server. Contigs were post-filtered on read alignment and the best-hit annotation was selected for each of them. Gene Ontology (GO) terms were predicted using Interproscan 5. All the data were available on a web interface RNAbrowse (http://ngspipelines2. toulouse. inra. fr:9005/ngspipelines/#!/NGSpipelines/Vitis%20 rupestris%20-%20RGMand1103P, last accessed 19 June 2017), which is described as an RNA-seq de novo assembly results browser  (Mariette et al. , 2014). The best annotation file was downloaded and a BIN code assignment was performed using Mercator  (Lohse et al.
Quality evaluation To assess the reliability of the de novo assembled transcriptome, raw reads were aligned against it using the aligner BWA MEM  (Li and Durbin, 2009). The completeness of the assembly was determined through the use of the Core Eukaryotic Genes Mapping Approach (CEGMA) (http://korflab. ucdavis. edu/Datasets/cegma/, last accessed 20 July 2016)  (Parra et al. , 2009) , which allows comparison between the transcriptome and a dataset containing highly conserved annotated core proteins.
Network analysis A co-expression gene network was constructed using the WGCNA software package (v1. 51) in R  (Langfelder and Horvath, 2008 ) using all the libraries except those corresponding to time 0 hpt (24 libraries). In order to remove low-expressed contigs (reflecting noise), a filter was applied to keep only contigs that had at least 20 counts in 70% of the libraries. A total of 45 358 contigs satisfying the above threshold were obtained, and their counts data were transformed using the function varianceStabilizingTransformation of the package DESeq2  (Love et al. , 2014). The resulting set of counts was used for network construction and module detection using the function blockwiseModules. Briefly, an adjacency matrix was created by calculating the biweight mid-correlation raised to a power β of 8 (softthreshold estimated with the pickSoftThreshold function) and the maxPoutliers parameter set to 0.
Therefore, a de novo assembly strategy was chosen in order to build the most comprehensive transcriptome allowing differential expression analysis between the two studied genotypes. Briefly, a total of 1. 19 billion cleaned reads from the 30 cDNA libraries were assembled into contigs using the Oases assembler, with an RPKM threshold of 1, which allowed a higher percentage of contigs showing protein assignments (Supplementary Table  S2 ). A total of 51 250 contigs (47 519 774 bp) were generated, with an N50 at 1481 (Supplementary Table  S2 ). The assembly generated a high number of transcripts enriched for smaller-sized contigs (51. 56% of contigs were in the size range 201-599 bp; Supplementary Fig. S1 ). To assign putative gene functions and annotations to all the assembled contigs, BLASTX alignment was performed against protein databases such as Swissprot, Refseq, and several Genotoul bioinformatics platform internal databases of different species, including Medicago truncatula, Theobroma cacao, Populus trichocarpa, A.
In total, 56 268 annotations were generated, corresponding to 51 012 contigs covering almost all the transcriptome (Fig. 3 ). In order to examine the representation of the transcriptome of the RNA-seq reads, reads were mapped back to the de novo assembly, and showed more than 90% of mapping for all the samples (data not shown). Finally, the completeness of the de novo assembled transcriptome was evaluated by comparison of its 51 250 contigs with a set of highly conserved core proteins using the CEGMA pipeline  (Parra et al. , 2009) ,  allowing the identification of 98. 39% of the 248 core eukaryotic genes contained in this core dataset.
13, NPF3. 1, and NPF4. 5). The expression profiles of seven genes with contrasting patterns were verified by qPCR. NR,  NIR, NPF6. 3, NRT2. 4a, NRT2. 4b, NRT3, and GS2  were differentially expressed in both combinations, following the same pattern of expression ratio between the HN and LN conditions (Fig. 7 ). The ratio of transcript abundance between HN and LN roots was 0. 82-to 6. 73-fold depending on the gene, and was higher in CS/RGM (average ratio=3. 04) than in CS/1103P (average ratio=2. 12). Pairwise Pearson's correlations between LFC values from edgeR analysis on RNA-seq data and those from qPCR data showed a strong correlation, indicating that the two methods gave very similar profiles (Fig. 7 ). Interestingly, the pronounced response observed for CS/RGM was found to be correlated with a higher difference of nitrate content between each side of the root system.
Weighted gene co-expression network analysis to identify hub genes in response to N availability A gene co-expression network analysis was conducted using the WGCNA R package  (Langfelder and Horvath, 2008). This approach was used to define clusters of highly correlated genes (modules). Including the grey module (which contained the unconnected genes), 26 modules were found and analysed for their association with each experimental trait (correlation with genotype, N condition, or time) (Fig. 9 ). Then, a kME value (module eigengene-based connectivity) was calculated for each transcript to every module (Supplementary Table  S5 ). For each module eigengene, highly correlated transcripts were filtered with a correlation coefficient >0. 80 and a P-value <0. 01. The first 100 transcripts were selected to perform an enrichment analysis using Mefisto; the resulting enriched categories are summarized in Supplementary Table  S6.
Heterogeneous nitrate supply modulates the root transcriptome in the two scion/rootstock combinations in different ways The differential analysis of transcript levels between HN and LN roots at 3 and 24 hpt highlighted that gene expression was more profoundly impacted by the difference in nitrate availability in CS/RGM than in CS/1103P. Nitrate measurements confirmed these results and suggested a close regulation of nitrate uptake, depending on its availability, for CS/RGM. In addition, a WGCNA analysis defined modules containing genes for which the expression profiles were associated to the different experimental traits and identified genes that could act as a hub in more than one module. A common set of genes exhibiting the same transcriptional regulation in both genotypes, that is, up-or down-regulation in HN roots compared with LN roots, was identified.
Transcriptional regulation of the ethylene and strigolactone pathway genes is rootstock dependent Genes belonging to the ethylene-related functional category were significantly over-represented in the DEGs between LN and HN roots in both genotypes, but with a time-dependent expression pattern. Most of these genes encode members of the apetala2/ethylene response factor (AP2/ERF-TF) super family. A recent study showed that the expression of these genes was regulated by N availability  (Zhao et al. , 2015). The interaction between ethylene and N affects several physiological processes, including root architecture  (Khan et al. , 2015). In addition, in Arabidopsis, ethylene regulates some nitrate transporters, such as NPF6. 3 and NRT2. 1  (Tian et al. , 2009). The WGCNA approach confirmed that ERF genes were found in modules associated with the LN treatment.