dataesr/scientific-paragraphs-categorization huggingface dataset:
39. 5. The ISA, proposed by Dufrêne and Legendre  28  , is based on the relative frequency of a specie (OTU in our case) within and inter treatments. It gives provides an indicator value ranging from 100 (perfect indicator) to 0 (no indicator). The perfect indicator species would be the one found in high abundance in all samples of a given treatment and absent in all other treatments. It also uses a randomization test to evaluate the significance of the specie distribution. Nucleotide sequence accession numbers. The data generated were deposited in the NCBI Sequence Read Archive (SRA) and are available under Bioproject accession numbers PRJNA471919. Gene quantification by qPCR. Genes were chosen based on their involvement in the soil nitrogen and methane cycles. Primers were selected according to the literature and assessed with the Primer blast tool of the National Center of Biotechnology Information (NCBI) (https://www.
Impact of land use change and Eucalyptus logging on microbial community structure. We assessed a total of 525,420 sequences (35,028 per sample) after applying quality filters and data normalization, clustered into 6,831 OTUs (3% dissimilarity threshold). Rarefaction curves show that our sequencing effort describes well the diversity of each sample (Supplemental Fig. 1 ). Bacterial richness (represented by numbers of OTUs) was significantly different among treatments, being highest in YE, followed by OE, and lastly NF (Table  4 ). Shannon diversity indices also differed significantly among treatments, being higher in YE and OE than in NF (Table  4 ). No significant difference was found between the two Eucalyptus plantations. Taxonomic assignments revealed that all treatment areas were dominated by the same phyla, but with significant differences in the abundances of some phyla among treatments (Fig.
Library preparation and data alignment We prepared scM&T-seq libraries 10 by isolating mRNA on magnetic beads and separating from the single-cell lysate as described  30  prior to reverse transcription and amplification using Smartseq2  31  but with 25 PCR cycles. We then processed the lysate containing genomic DNA according to the published single-cell bisulfite sequencing protocol  32. Single-cell RNA-seq libraries were aligned using HiSat2 with options --sp 1000,1000 --no-mixed --no-discordant  33. Single-cell bisulfite libraries were processed using Bismark 34 as described  10. Mapped RNA-seq data were quantitated using the RNA-seq quantitation pipeline in Seqmonk software (www. bioinformatics. babraham. ac. uk/projects/seqmonk/).
Data analysis RNA-seq Gene expression levels were estimated in terms of reads per million of mapped reads to the transcriptome. A score of variability per gene (named distance to the median) was calculated by fitting the squared coefficient of variation as a function of the mean normalized counts and then calculating the distance to a rolling average (window size=100) (Fig. S8 )  24. We included only genes with an average normalized read count of at least 10. The top 1000 most variable genes of the entire data set were used to perform principal component analyses (as log 2transformed and median-cantered values) (Fig. 1B , Table  S3 ). Single cell differential expression (SCDE) was used to calculate differential expression analysis between young and old cells (Table  S1 )  35. Cell-to-cell correlation analyses were performed using the top 500 most variable genes within each individual and using Spearman's correlation as the measure of similarity between cells (Fig.
DNA-methylome We discarded cells that had less than 1 million paired-end alignments or less than 500,000 CpG sites covered (Fig. S1 ). To avoid biases that might occur due to different sequencing depths or number of cells between individuals, we down-sampled the data to 1 million reads for each cell and randomly selected 35 cells from each individual (2 young and 2 old). ChIP-seq datasets for H3K4me3, H3K27me3, H3K36me3 in satellite cells and H3K27ac in myoblast were obtained from existing studies  23, 29. Bowtie2 and MACS2 were used for mapping and peak calling respectively.
Acknowledgments: We would like to thank the  Flow Cytometry Platform of the Center for Technological Resources and Research (Institut Pasteur)  and the  Wellcome Trust Sanger Institute sequencing facility  for assistance with Illumina sequencing. This project was supported by grants from  Institut Pasteur ,  Agence Nationale de la Recherche  ( Laboratoire d'Excellence Revive, Investissement d'Avenir ;  ANR-10-LABX-73 ),  Association Française contre les Myopathies  ( 21857 ),  CNRS , the  European Research Council  ( Advanced Research Grant   332893 ), and the  Biotechnology and Biological Sciences Research Council (BBSRC ,  CBBS/E/B/000C0425 ).
Results The promoter interactome of human islets. To create a genome-wide map of long-range interactions between gene promoters and distant regulatory elements in human pancreatic islets, we prepared Hi-C libraries from four human islet samples. We performed hybridization capture of 31,253 promoter-containing HindIII fragment baits and their ligated DNA fragments, which were then sequenced and processed with the CHiCAGO algorithm to define promoter-interacting DNA fragments  20, 21  (Figure  1a, b, ). This resulted in 175,784 high-confidence interactions (CHiCAGO score > 5) between annotated promoters and distal genomic regions. We independently validated pcHi-C landscapes at two loci by 4C-seq analysis in the EndoC-βH1 human β cell line  22  (Supplementary Figure  1j, k ). In parallel, we created new human islet ChIP-seq and ATAC-seq datasets, and refined human islet epigenome annotations to subclassify active enhancers according to Mediator, cohesin, and H3K27ac occupancy patterns (Figure  1b-d,  and  Extended Data 1.
These studies, therefore, identify functional target genes of an enhancer that harbors diabetes-associated variants, and points to functional connectivity of two enhancer hub promoters. To further explore the behavior of hubs as functional domains, we again examined islets exposed to moderately low vs. high glucose concentrations, using above-mentioned datasets. Glucose-induced enhancers and mRNAs were highly enriched in hubs, compared with non-hub genomic regions (Fisher's P = 1. 1 X 10 -7 and 2. 2 X 10 -16 , respectively). Of 297 promoters that showed increased H3K27ac at high glucose, 94 were contained in hubs, and 65% of these also showed a significant glucose-dependent mRNA increase (Supplementary Table  3 ). If hubs are indeed regulatory domains, we expect that hub enhancers connected to glucose-induced genes should show concordant glucose-dependent changes.
Common variants in islet hubs provide tissue-specific risk scores. The improved definitions of genomic regions relevant to islet function prompted us to examine if islet hub variants could be harnessed to identify individuals in whom variation in islet function, rather than other T2D-relevant cellular functions, plays a preponderant role in T2D susceptibility. Recent studies suggest that polygenic risk scores (PRS) that integrate the combined effects of very large number of variants, including many that lack genome-wide significant association, can be used to identify individuals at risk for polygenic diseases  34, [40] [41] [42] 45. The predictive power of PRS for T2D is likely to increase as larger GWAS datasets become available, although it is already possible to explore the potential of hub variants to define tissue-specific specific risk profiles.
We first calculated a PRS using all common variants from a recent BMIadjusted T2D GWAS meta-analysis  46  , and examined the ability of this genome-wide PRS to predict T2D in the UK Biobank population cohort  47, 48. The UK Biobank cohort was divided into a testing dataset used to train the PRS model, and a validation dataset of 236,236 individuals. This showed that individuals of the validation dataset with the 2. 5% highest PRS had a 7. 1-fold higher frequency of T2D than those in the lowest 2. 5%, and an incident T2D hazard ratio of 2. 5 relative to the entire population (Figure  6f , Supplementary Figure  7h ). Next, we built PRS models with islet hub enhancer and promoter variants. Although hub regulatory regions represent < 2% of genomic space, 2. 5% of individuals with highest hub PRS scores had 4-fold higher frequency of T2D than those in the lowest bin and an incident T2D hazard ratio of 2.
Genome-wide PRS Pseudo-hubs PRS (n=100 sets)  (f) Impact of islet hub polygenic risk scores (PRS) on T2D frequency. T2D frequency (y-axis) was calculated in 40 bins, each one representing 2. 5% of individuals in the UK Biobank validation dataset. PRS values were calculated with the entire set of common genetic variants (light blue dots), the fraction of islet non-hub variants (dark blue dots) and islet hub variants (orange dots), or variants from 100 sets of pseudo-hubs (boxplots with IQR). (g) Polygenic risk stratified by BMI. The T2D risk ratio was calculated as the frequency of T2D cases in the highest 2. 5% PRS bin divided by the frequency of T2D cases in the lowest 2. 5% risk bin. Boxplots show the risk ratio for PRS from 100 sets of pseudo-hubs, the dark blue dot corresponds to the islet non-hub PRS and the orange dot shows that of islet hub PRS.
Controls were stratified by age at recruitment. Statistical significance for islet hub OR was calculated using a Z-score defined as the number of standards deviations above the average OR of the pseudo-hub distribution. See also Supplementary Figure  6 , Supplementary Table  4  and Data   (g) Tissue-selectivity of islet pcHi-C interactions. pcHi-C datasets from four primary hematopoietic cell types (erythroblasts, macrophages, naïve CD4+ T cells and total B lymphocytes) were used to classify islet interactions as tissue-invariant, if they were also observed in at least four datasets, or islet-selective, if they were observed in islets and no other datasets. Intermediate selectivity was classified as others. (h) Islet-selectivity of interactions is concordant with islet-selective gene expression.
5 µl of 25 µM Forward/Reverse ATAC-seq index primers from Nextera kit (FC-121-1030) and 20 µl of the purified transposed DNA. We used the qPCR plot on ABI7900 (Applied Biosystems) to determine the number of additional cycles of PCR amplification that were required, by using the cycle number that corresponds to 1/3 of the maximum fluorescent intensity. ATAC libraries were purified using Qiagen MinElute PCR Purification Kit and qPCR of open chromatin sites and negative controls were used before sequencing (data not shown). ChIP-seq and ATAC-seq analysis. Related to Figure  1. Illumina TruSeq adapters were removed from ChIP-seq reads using cutadapt 1. 9. 1 (options: -m 20)  12  In ATAC-seq reads, low quality bases were trimmed using Trimgalore 0. 4. 1 (options -quality 15 --nextera), which also removes Nextera transposase adapters (https://github.
com/FelixKrueger/TrimGalore). Trimmed reads were aligned to hg19 genome build using bowtie2 2. 1. 0 (options: --no-unal) allowing no mismatches  13. Aligned reads were filtered to retain only uniquely mapped reads (MAPQ>=30) using samtools 1. 2  14  and duplicate reads were removed using picard 2. 6. 0  15. Reads mapping to blacklisted regions  16  were also removed using BEDTools 2. 13. 3  17. For ATAC-seq, reads mapping to mitochondrial genome were also removed. Mapped reads obtained at this stage are referred to as usable reads. The quality of ChIP data was assessed with SPP. R script from phantompeaktools  18. The list of ChIP-seq and ATAC-seq experiments from human islet samples, alignment statistics and accession numbers are shown in Supplementary Table  8. MACS2 was used to determine regions of statistically significant enrichment over corresponding input samples.
Aggregation plots for pcHi-C baits and interacting regions. Related to Figure 1d and Supplementary Figure  1d. A window of 25 kb around the center point of HindIII fragments containing pcHi-C baits and promoter-interacting regions was selected for histone marks and transcription factor density calculations. The distance density was estimated using Gaussian kernels (python 2. 7 function scipy. stats. gaussian_kde)  26. The expected distribution was computed using 10 randomizations of the epigenomic feature coordinates across the mappable nonblacklisted genome using shuffleBed from BEDTools  17  , retaining only nonoverlapping random regions. Density plots were computed using computeMatrix and plotProfile tools from Deeptools 2  27. Aggregation plots of Mediator and H3K27ac ChIP-seq signals from three human islet donor samples were plotted using computeMatrix and plotProfile tools from Deeptools 2  27.
Input DNA was used as a reference. To create control interactions, we used the method implemented in the CHiCAGO package to create 100 sets of distance-matched interactions  8. Briefly, promoter interacting regions were randomly re-organized across all 22,076 baits such that the distance distribution of the interactions was maintained. The overlaps were also defined through the CHiCAGO package as the number of promoter interacting HindIII fragments that overlap with the genomic regions of interest, using either the set of high-confidence interactions (CHiCAGO score>5), giving us the observed overlap, or the sets of distance-matched control interactions, giving us the expected overlap. For the expected overlaps, a 95% confidence interval of the mean was also obtained. 2k-l. Epigenomic features that showed enrichment in high-confidence interacting regions (e.
Expression selectivity analysis of islet-selective and tissue-invariant pcHi-C interactions. Related to Supplementary Figure 1h. To assess whether isletselective pcHi-C interactions were preferentially connected with genes showing isletselective expression, we selected baits that exclusively contained islet-selective or tissue-invariant interactions and computed the gene expression islet specificity Zscore of the genes contained in the baits. Statistical analysis was performed using Wilcoxon's signed ranked test. 4C-seq Library Preparation. Related to Supplementary Figure  1j, k. 4C-seq library preparation of human EndoC-βH1 β cells was performed following the protocol described in  9. 4C-seq datasets from native human pancreatic islets in MAFB and ISL1 loci shown in Supplementary Figure  1j ,k are taken from  9.
𝐷𝐼 = ! 𝐵 -𝐴 |𝐵 -𝐴| ! ! (𝐴 -𝐸) ! 𝐸 + (𝐵 -𝐸) ! 𝐸 ! 𝐸 = (𝐴 + 𝐵)/2 Identification of Promoter-Associated Territories (PATs). Related to Supplementary Figure  2f, g. PATs were created for imputing enhancer-promoter interactions. They were defined as the linear space covered by all the interactions originating from a pcHi-C bait, within the same islet TAD-like compartment. A total number of 16,030 PATs were defined in this manner (see also Supplementary Figure  3e ). ChromHMM segmentation of regulatory elements. Related to Supplementary Figure  2h-j. To segment the genome into chromatin states based on combinations of chromatin marks in pancreatic islets, we used ChromHMM v1. 11  31. Previously published human islet in-house datasets (H3K27me3, H3K36me3, H3K4me1, H3K4me3, H2AZ)  9  , externally generated H3K9ac and H3k9me3 datasets  32  , and newly generated datasets described in the current study (cohesin, Mediator, CTCF and ATAC) were used.
The specific parameters used were as follows: reads were shifted in the 5′to 3′direction by 100 bp. For each of the aligned ATAC-seq and ChIP-seq datasets, read counts were computed in non-overlapping 200-bp bins across the entire genome, which were binarized into 1 (enrichment) or 0 (no enrichment) by comparing the read counts within each bin against a background signal using a Poisson p value threshold of 1 × 10 -4. Models from 10, 15, 20 and 25 states were trained in parallel. We chose to focus on the 15-state model as it best summarized the pancreatic islet epigenome. ChromHMM enrichment in PATs. Related to Supplementary Figure  2h, i. For each PAT we considered the entire linear genomic space, and computed the fraction of genomic space occupied by every ChromHMM state segment. The enrichment of ChromHMM states in PATs was calculated as the fraction of PAT space occupied by a given ChromHMM state divided by the fraction of genomic space occupied by the same state.
For this analysis, we used 14 human pancreatic islet samples, including 7 samples exposed to 11mM glucose and 4mM glucose, and 51 tissues from Epigenome Roadmap Consortium consolidated epigenomes. Epigenome RoadMap aligned reads from H3K27ac ChIP-seq samples and corresponding inputs were downloaded from egg2. wustl. edu/roadmap/data/byFileType/alignments/consolidated/ and converted to BAM format using bamToBed from BEDTools 17. To avoid artificial bias in H3K27ac signal strength due to differences in sequencing depth, all datasets selected contained at least 15 million usable reads and were uniformly subsampled to a maximum depth of 30 million usable reads. Active islet enhancers were uniformly extended to +/-750 bp from the center of the peak. Active islet promoters were used without further modification.
The statistical significance was assessed performing a chi-square test (python 2. 7 function scipy. stats. chi2_contingency) comparing the frequency of each gene expression class among the two "assigned genes" and "control genes" lists. (c) Concordant glucose regulation of assigned enhancer-promoter pairs. Related to Figure  2c, d. H3K27ac ChIP-seq and RNA-seq datasets were analyzed in islets from 7 human donors that were cultured for three days in either 11 or 4 mM glucose, as described above. To assess glucose-regulation of enhancers, we defined H3K27ac-enriched regions, rather than using annotated enhancers, which typically contain nucleosomedepleted subregions that do not show H3K27ac enrichment. We thus defined consistent H3K27ac-enriched regions for each human donor and each glucose condition treatment using MACS2 with the parameters described above.
A total number of 90,814 narrow H3K27ac-enriched regions were interrogated for this analysis. The number of H3K27ac reads mapping to each peak was calculated using FeatureCounts v1. 5. 0 program (--ignoreDup -O --minOverlap 10). Then, paired DESeq2 (v1. 10. 1) analysis was used to assess differential signal strength. Peaks consortium. org/downloads. html)  34 and (b) the re-analysis of T2D publicly available GWAS datasets and imputation with 1000 Genomes Project (1000G) and UK10K  35. For loci that were examined in both studies, we considered the union of all variants. Additionally, a list of T2D and Fasting Glycemia-associated lead SNPs was downloaded from the NHGRI-EBI GWAS Catalog on 31/05/2016, from which we selected SNPs if they complied with the following criteria: (a) genome wide significance (p ≤ 5 x 10 -8 ) in a study that could be verified in a peer-reviewed publication; (b) the phenotype in the study was reported as Type 2 diabetes, or Type 2 diabetes (and other traits), for T2D-associated SNPs, "Fasting glucose-related traits" or "Fasting plasma glucose", and (c) they were not identified in or were in high linkage disequilibrium (LD) with SNP in the credible sets from  34  and/or  35.
Epigenomic features interrogated in logistic regression analysis: --H3K4me3 signal at TSS in human pancreatic islets --H3K27me3 signal at TSS in human pancreatic islets --H3K9me3 signal at TSS in human pancreatic islets Logistic regression (LR) was implemented using python 2. 7 and scikit-learn  48  with class_weight mode set as "balanced" to automatically adjust weights in an inversely proportional manner to the class frequencies in the input data. The machine learning classifier was trained using 70% of the full dataset with 250 iterations and random sampling. The remaining 30% of the dataset was used for validations. In each iteration, the machine learning classifier computed the logistic regression coefficients,  ! !∈ !;! , for a list of  epigenomic features listed above,  ! !∈ !;! , that fit a logit function.
Generation of 3D ensemble structures of enhancer hubs. Related to Figure  5b, c  and Supplementary Figure  5b, c. Normalized interaction matrices were modelled using TADdyn, a molecular dynamic-based protocol implemented in the TADbit software  52. Similarly to TADbit, TADdyn generates models using a restraint-based approach, in which experimental frequencies of interaction were transformed into a set of spatial restraints, as previously described  53  , but suited for sparse datasets such as pcHi-C. The underlying chromatin fiber was represented by a bead-spring polymer model  [54] [55] [56] [57] [58]  , where the size () of each bead was defined by the relationship 0. 01 nm/bp assuming the canonical 30 nm fiber  59, 60  , as implemented in TADbit. The system scoring function (H) was defined as:  =  !"# +  !"#$!! where H pol is the contribution of the polymer and H pcHi-C is defined by a set of distance restrains proportional to the experimental frequencies of interaction.
Furthermore, TADbit tools was used to calculate the mean distances particles containing genomic regions of interest in the model. Using the mean distances, a promoter-enhancer interaction network was built as a weighted undirected graph, in which each node represents either a promoter or an enhancer in a specific locus and each edge is weighted by its mean distances values. The promoter-enhancer interaction network was decomposed into communities using MCODE clustering algorithm as implemented in clusterMaker2  64  and visualized with Cytoscape 3. 5  65. Spatial distribution of chromatin marks in 3D models. Related to Figure  5d  and Supplementary Figure  5d. The 3D spatial distribution of ChIP-seq data from human islets was analyzed as follows: for a given central viewpoint, an initial sphere with a radius of 200 nm was constructed from the viewpoint.
Then, a set of spherical shells, that occupied a volume equal the initial sphere, were added. Each particle neighboring the viewpoint was assigned to a spherical shell based on their relative distance to the viewpoint. ChIP-Seq data from MACS2  66  was then binned at 5 kb resolution and the mean score value of their coverage score was calculated. The odd ratio was calculated for each shell using a Fisher's exact test for 2 × 2 contingency tables comparing bins with and without ChIP-Seq signal. A total of 13 different ChIP-seq datasets of human pancreatic islets (CTCF, FOXA2, H3K27me3, H3K27ac, H3K36me3, H3K4me1, H3K4me3, MAFB, NKX2. 2, NKX6. 1, PDX1, Mediator and cohesin) were analysed and the corresponding circular 2D plot was generated. SNP enrichment analyses. Related to Figure  6  and Supplementary Figure  6.
To assess class I enhancers inside hubs we generated Q-Q plots and λ measures for variants overlapping (i) class I enhancers inside hubs, (ii) open chromatin regulatory elements outside hubs and (iii) the rest of variants that did not overlap either of the previous categories. GWAS meta-analysis of Insulinogenic index from oral glucose tolerance test (OGTT) measures. Related to Figure  6. 7,807 individuals from four population studies were included in these analyses: the Inter99 study (ClinicalTrials. gov ID-no: NCT00289237) which is a population-based non-pharmacological intervention study for ischemic heart disease (n=5,305)  69  , the Health2008 cohort (n=605)  70  , the 1936 Birth Cohort (n=709)  71  and the ADDITION-Pro cohort (n=1,188)  72. All study participants gave informed consent and the studies were approved by the appropriate Ethical Committees and were performed in accordance with the scientific principles of the Helsinki Declaration II.