dataesr/scientific-paragraphs-categorization huggingface dataset:
4c  and Extended Data Fig. 3 ), indicating that the differences between exAMs and in vivo AMs disappeared upon transplantation. Comparison to tissue-resident macrophage populations from published datasets 4 using principal-component analysis (PCA) indicated that texAMs were indistinguishable from tAMs or resident host AMs, but different to all other resident macrophage populations (Fig. 4d ). As an indication of the high degree of their similarity, tAMs, texAMs and host AMs were closer to each other than to the published AM datasets from another laboratory  4  , generated with a different protocol (Fig. 4d ). In fact, only 56 genes were significantly differentially expressed with more than a twofold change (FDR < 0. 05), between tAMs and texAMs, compared to 3,547 genes that were differentially expressed between exAMs in culture and texAMs (Fig.
g, Flow cytometric analysis showing expression of cell-surface markers F4/80, CD64, CD11c, SiglecF and CD11b (colored) against isotype control (gray) on AMs and exAMs (1 month). h, RNA-seq gene expression from two replicates each pooled from three mice/cultures showing core AM upregulated signature genes  3  in AMs or exAMs (2 months) compared to peritoneal macrophages (PMs) and AMs from the Immgen consortium  48  (GSE122108). The box plot shows the median with the lower and upper quartiles representing a 25th to 75th percentile range with whiskers extending to 1 × IQR and outliers falling outside of this range. Pairwise comparison between average PMs and other samples by one-sided Wilcoxon signed-rank test Benjamini-Yekutieli adjusted P values <10  -18. i, Heat map of AM-specific transcription factor gene expression from two replicates each pooled from three mice/cultures of AMs or exAMs (2 months) compared to PMs and AMs (Immgen), showing rlog normalized values from minimum (blue) to maximum (red).
FACS data acquisition was conducted with a BD LSR Fortessa and analysis was performed using FlowJo software (v. 10 Treestar). Phagocytosis assay using fluorescent latex beads. AMs were incubated in growth medium containing fluorescent yellow-green, carboxylate-modified polystyrene latex beads (L4655, Sigma-Aldrich) at 1:400 dilution for 1 h at 37 °C. Supernatant medium containing beads was discarded post-centrifugation and the cells were washed three times for 5 min with PBS. After cytocentrifugation, cells were fixed using 4% PFA for 10 min and permeabilized with 0. 1% Triton-X in PBS for 15 min at room temperature. Cells were incubated in blocking buffer (5% donkey serum + 1% BSA in PBS) for 15 min at room temperature followed by Actin Red staining (two drops per ml in PBS). All the steps were carried out with washing steps in between (three times, 5 min each using PBS).
Libraries were sequenced with 75 paired-ends using NextSeq 500 instrument (Illumina) for an average of 2 × 10 7 reads per sample. Bioinformatic analysis. RNA-seq analysis. The following analyses were conducted with R (v. 3. 3). Samples were aligned to mouse genome version mm10 using Bowtie aligner (v. 1. 1. 2). A table of raw counts was generated using htseq-count (v. 0. 6. 1) and R (v. 3. 3). All samples were normalized to 1 M reads using the c. p. m. function from the edgeR package. To check reliability of samples replicates, we conducted PCA on normalized expression data on all expressed genes (excluding lowly expressed genes) and correlation tests on normalized expression data. For identification of DEGs, a filtering step was conducted depending on the comparison of interest. Genes having more than 4 c.
p. m. in at least one condition in a pairwise comparison were kept for the remaining analysis. For detection of DEGs, raw data were treated using DESeq package (v. 1. 20). GO analyses were conducted on DEGs after k-means clustering between exAMs versus host AMs with the help of the ChIPseeker package (v. 1. 26). R software (v. 3. 3) was used to generate heat maps. ATAC-seq analysis. Data were pooled from three independent experiments. Fastq file quality control was assessed by the fastqc tool. Adaptor trimming was performed using Cutadapt (v. 2. 6). Bowtie2 (v. 2. 3. 5. 1) was used to align and index sequences on the mm10 mouse genome and thus generated the first series of bam files. Then, peak detection was made on these bam files by Genrich (v. 0. 6. 1). Then, from the first series of bam files, a second series was generated as following: after using Bowtie2 (v.
2. 3. 5. 1), mitochondrial DNA was removed with Samtools (v. 1. 9) and Picard (v. 2. 24. 2) was used to remove PCR duplicates. Finally, Samtools allowed to remove non-unique alignments. Bigwig files were generated from BedGraph files using a bdg2bw script (https://gist. github. com/taoliu/2469050). These BedGraph files were created from the second series of bam files. All peaks detected by Genrich from exAM, texAM and host AM BAL samples were pooled and merged together with bedtools (v. 2. 26. 0). Then a table of counts was generated by bedtools from these peaks. The signal was computed from the second series of bam files. Then peaks with >10 counts in at least one condition were kept for the following analysis. Opened and closed peaks were detected with the DESeq package (v. 1. 20). Peaks with adjusted P values <0.
05 were declared as differentially open between two conditions. Comparisons were texAMs versus host AMs, exAMs versus host AMs and exAMs versus texAMs. Then log 2 (count + 1) was represented in scatter-plots comparing conditions two by two. Peaks with log 2 FC computed by DESeq2 as >1. 5 (and <1. 5, respectively) were highlighted in red (respectively in blue). For motif enrichment analysis, 1,492 peaks were selected as gained in exAMs (as compared to host) and 937 were selected as lost. Then motif analysis was performed on these two groups of peaks using HOMER (v. 4. 10) (http://homer. salk. edu/homer). HOMER, inter alia, is a discovery motif analyzer and screens for enrichment of known motifs. HOMER perl script findMotifsGenome. pl was used with the mm10 mouse genome as a background (random genomic sequences samples according to GC content of input sequences).
Venn diagrams were plotted using VennDiagram package (v. 1. 6. 20). GO analysis and gene associations with peaks were performed with GREAT tool (v. 4. 0. 4) using default options to assign regions and using minimum region-based fold enrichment of 1. 2 and term annotation counts (5-1,000) were set as global controls for GO term analysis. For peak visualization reference ATAC-seq sample bigwig files were downloaded from the Immgen consortium 50 (GSE100738). Peaks were visualized on our samples independently or in parallel with the Immgen references using IGV (v. 2. 10. 2) and were auto-scaled by group. Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Acknowledgements We thank  L. Chasson  for histology (CIML),  M. Barad ,  S. Bigot ,  A. Zouine  (CIML),  A. Gompf  and  K. Bernhardt  ( CMCB Technology Platform, TU Dresden ) for flow cytometry,  M. Sohn  for sequencing and the  K. Rajewsky laboratory  for hosting and general support (MDC). We also thank  L. Razafindramana ,  CIML and CRTD Animal  facilities for mouse husbandry and the light microscopy facility of the  CMCB technology platform, TU Dresden. This study was supported by institutional grants from  TU Dresden , the  Institut National de la Santé et de la Recherche Médicale ,  Centre National de la Recherche Scientifique  and  Aix-Marseille University  and grants to M. H. S. from the  Agence Nationale pour la Recherche  ( ANR-17-CE15-0007-01  and  ANR-18-CE12-0019-03 ),  Fondation ARC pour la Recherche sur le Cancer  ( PGA1 RF20170205515 ), an  INSERM-Helmholtz cooperation  and the  European Research Council  under the European Union's  Horizon 2020 Research and Innovation program  (grant agreement no.
5 × IQR whiskers). Individual samples are colored according to the human islet dataset they belong to (see color legend in Fig. 1a ). Nominal QTL p-values are provided. d Enrichment of sQTL and eQTL variants in different functional genomic annotations using GREGOR. The x-axis represents GREGOR fold change of observed vs. expected number of SNPs at each functional annotation. The dotted line represents 1. 5-fold change. e Percentage of eGenes and Junctions with eQTLs or sQTLs at FDR ≤ 1%, respectively, shared in different number of GTEx tissues. We used significant eQTLs and sQTLs from 47 distinct GTEx V8 release tissues a significant islet-specific component of sQTLs (Fig. 2e ). Taken together, our results reveal two separable layers of genetic regulation of the human islet transcriptome.
We applied colocalization as implemented in gwas-pw  [24] , which draws upon the original coloc algorithm but does not rely on user-defined priors. This identified candidate effector transcripts with high colocalization support (posterior probability of shared association between both phenotypes ≥ 0. 8) at 9 independent T2D GWAS signals using sQTLs, and 25 using eQTLs (Additional files 7 and 8: Table  S6 -S7, Additional file 1: Fig. S6a ). At loci Fig. 3  Role of human islet splicing variation in T2D susceptibility. a Quantile-quantile (QQ) plot showing observed T2D association p-values in human islet sQTLs (orange dots) and eQTLs (blue dots) against p-values under the null hypothesis. The grey-shaded region represents 1000 p-value distributions (in the -log 10 scale) of random sets of control sQTL variants (see the "Methods" section).
In this approach, which so far has not been applied to human islet RNA-seq data, genetic effects on splicing or expression are used to impute transcript variation in cases versus controls from GWAS datasets. This allowed us to identify splicing or expression changes associated with T2D and related traits. More specifically, we used the FUSION algorithm  [25]  and GWAS summary statistics  [1]  to identify T2D associations with islet junction usage or gene expression (sTWAS and eTWAS, respectively). Because TWAS findings do not necessarily distinguish between shared genetic effects and linkage  [26, 27] , we focused on TWAS signals showing colocalization with the GWAS phenotype (PP4 ≥ 0. 6), thus minimizing confounding effects from linkage. This identified 27 genes (42 splicing events) showing significant sTWAS with T2D risk, and 29 genes with eTWAS after multiple testing correction (Bonferroni p= 8.
Islet QTLs provide insights into T2D pathways The availability of an expanded list of candidate effector genes, as opposed to lists of genes that are located in the vicinity of associated SNPs, allowed us to explore the hypothesis that a subset of genetic signals could influence T2D susceptibility by acting on specific cellular pathways in pancreatic islet cells. We thus compiled 106 putative effector genes for T2D or glycemic traits (FG, FI) reported here as well as previously reported co-localizing islet eQTLs  [12]. This list excluded those exclusively detected in non-endocrine cells after analysis of single cell RNA-seq datasets (Additional file 13: Table  S12 ), as well as non-coding transcripts and genes without known function. Functional gene annotations revealed notable enriched pathways, including regulators of fatty acid biosynthesis (SCD5, FADS1, GCDH, BDH2, PAM; GO:0006636, GO:1901570, ENRICHR adjusted p = 0.
04), and genes upregulated by hypoxia or mTORC1 activation (MsigDB Hallmark Hypoxia and mTORC1 signaling; q < 0. 05 and < 0. 001, respectively). We further used the list of 106 genes to identify networks using STRING  [59]  v11. 5, using co-expression, experimental, and functional annotation databases. We allowed for minimal network inflation (≤5 interactions) and omitted text-mining to preclude bias arising from publications that name genes at significant GWAS loci. The resulting network exhibited 1. 8-fold greater proteinprotein interactions than random sets (PPI enrichment p < 10 -4 ) and contained two distinct subnetworks (Fig. 5a ). One subnetwork contained components of the eIF3 translational initiation complex (FDR = 0. 018). Several monogenic diabetes genes target translational initiation, in particular eIF2 complex (EIF2B1, EIF2S3, EIF2AK3), due to their impact on endoplasmic reticulum stress in β-cells  [60, 61]  (Fig.
Discussion This study adds splicing variation in pancreatic islets to the spectrum of molecular mechanisms that underlie T2D predisposition. Earlier studies had examined expression QTLs in human pancreatic islets  [10-13, 18, 19]. The current study offers the first systematic analysis of how common genetic variants influence RNA splicing in human islets. Parallel profiling of splicing and expression QTLs in the same dataset demonstrated that these represent two distinct mechanisms through which genetic variation can influence islet biology and disease. We found islet sQTLs that impact genes that have major roles in islet cell function as well as in the pathophysiology or treatment of various forms of diabetes. Furthermore, we observed a selective inflation of T2D association p-values among sQTLs, some of which showed fine co-localization with T2D variants.
Human pancreatic islet datasets We compiled 447 human islets samples, which comprised 89 samples from GEO accession number GSE50244  [19] , 118 from EGA accession number EGAD00001001601  [18] , 112 samples from T2DSystems EGA accession number EGAS00001005535  [13] , and 128 samples from CRG cohort (EGAS00001006440)  [70]. Out of 447 samples, 399 samples passed genotype and RNA-seq QC filters described below. For the CRG cohort, human pancreatic islets from organ donors without a history of glucose intolerance were purified using established isolation procedures  [71] [72] [73] [74] , shipped in culture medium and re-cultured at 37 °C in a humidified chamber with 5% CO 2 in glucose-free RPMI 1640 supplemented with 10% fetal calf serum, 100 U ml -1 penicillin, 100 U ml -1 streptomycin, and 11 mM glucose for 3 days before analysis.
Genotype processing The 128 CRG cohort samples were genotyped with Illumina Infinium OmniExpress 12 v1 and HumanOmni 2. 5-8v1 on a total of 624K SNPs. The 89 samples from Fadista et al. (GSE50244)  [19]  were genotyped using Illumina HumanOmniExpress 12v1 C on a total of 609K SNPs. The 118 samples from Van de Bunt et al. (EGAD00001001601)  [18]  were genotyped using Illumina Omni2. 5+Exome array on a total of 2. 6M SNPs. Finally, the 112 samples from T2DSystems (EGAS00001005535)  [13]  were genotyped using Illumina's Human Omni 2. 5 exome array on a total of 2. 6M SNPs. A three-step quality control of genotype data, involving two stages of SNP removal and one intermediate stage of sample exclusion, was conducted in each cohort. Genotyped SNPs were filtered if (i) minor allele frequency (MAF) < 0.
HRC-1000G-check-bim. pl script with the -n option (to turn off the removal of variants showing MAF differences between the reference panel and the study genotypes) was used. We submitted resulting VCF files to the Michigan Imputation Server (https:// imput ation server. sph. umich. edu/ index. html): EAGLE2  [75]  was used for phasing, minimac3  [76]  for genotype imputation with the HRC  [77]  r1. 1 and the 1000G Phase 3 release reference panels  [78] , independently. For each dataset of imputed genotypes, we excluded variants with (i) MAF < 1%, (ii) imputation-quality R 2 < 0. 7, and/or (iii) HWE p-value ≤ 1×10 -6. We extracted indels from the 1000G Phase3 imputed results, filtered them using the aforementioned criteria, and merged with the filtered HRC imputed dataset.
Principal component analysis of genotypes To identify individuals of divergent ancestry and to characterize population structure of 399 human islet samples, we first selected a subset of genotyped SNPs that were common in all 4 cohorts, that also passed all our QC filters (see the "Genotype processing" section), and with MAF ≥ 1% and missingness < 5% across all the samples. We also excluded SNPs in high LD (pairwise r 2 ≤ 0. 1 within 1 Mb window), C/G and A/T SNPs to avoid strand mismatches, and those located in previously reported regions with long-range LD. We aggregated the 1000 Genomes Phase3 reference dataset using the set of overlapping variants. flashPCA  [82]  tool was used to calculate genetic principal components (PCs) (Additional file 1: Fig. S1c-d ). First four genetic PCs are used as covariates in subsequent sQTL and eQTL mapping analysis.
Gene expression quantification In-house developed human pancreatic islet transcriptome annotations built with a combination of CAGE as well as short and long sequence reads were used to quantify total gene expression  (Atla, G. , Beucher, A. et al. , unpublished). FeatureCounts  [83]  was used to obtain gene level qualifications using default parameters except using appropriate strandedness flag for each dataset. Genes with less than 5 raw reads mapped in less than 10% of the samples within each cohort were removed. CPM normalization was performed using edgeR  [84]  cpm function and then log2 transformed. Combat  [85]  was used to remove known sequencing batch effects. Principal components (PCs) were calculated on Combat corrected gene expression values using prcomp function in R (Additional file 1: Fig.
Splicing quantification To quantify splicing activity, we used the annotation free method, leafcutter  [20]. Briefly bam2junc. sh script from leafcutter was used to quantify junction spanning reads based on spliced alignments from bam files. We removed junctions that are not supported by at least 5 spliced reads in at least 10% of samples. We then used leaf-cutter_cluster. py script with options -m 30 -l 500000 to cluster the junctions anchoring on shared splice sites. This identifies local splicing events which are a cluster of alternatively used junctions. prepare_phenotype_table. py script was used to get relative junction usage (ratios) across samples. The relative junction usage is also referred to as percent spliced in (junction PSI) in the manuscript. Combat  [85]  was used to remove known sequencing batch effects and principal components were generated using prcomp function in R (Additional file 1: Fig.
cis-eQTL mapping cis-eQTL mapping was performed using QTLtools  [86]  for 399 samples with available QCed genotype and RNA-seq data using a cis-window of 500 kb up-and downstream of the transcription start site (TSS). Fifteen PCs derived from gene expression and 4 genetic PCs were used as covariates in the linear model. In order to identify best associated cis eQTL SNP-eGene pairs, QTLtools was run using the permutation pass mode using default parameters and --permute 1000 --window 500000 --seed 123456. Beta approximated permutation p-values were adjusted for multiple testing correction using Storey q-values implemented in the qvalue R package. We set the significance threshold at FDR q-value ≤ 0. 01. This resulted in 3433 genes (eGenes) with significant eQTLs (Additional file 2: Table  S1 ).
Annotation of sQTL junctions As leafcutter identifies junctions de-novo, we used transcriptome annotation GTF files as a base to annotate the junctions with respective genes. We used gtf2leafcutter. pl script from leafcutter's leafviz module (https:// github. com/ david aknow les/ leafc utter/ tree/ master/ leafv iz) to obtain the intron coordinates of all genes. The sQTL junctions were then mapped to the intron coordinates and annotated with respective gene names.
Comparison of eQTLs and sQTLs Our maps of genetic effects on splicing were compared to our eQTLs but also to previously reported exonQTLs and eQTLs in the largest study to date in human pancreatic islets by the InsPIRE consortium  [12]. To estimate the gain of novel information uniquely provided by our sQTL analysis, we first determined the overlap between sGenes and eGenes in our study, and eQTLs and exon QTLs mapped by InsPIRE. For those genes that contain both sQTLs and either eQTLs or e/exon QTLs in InsPIRE, we calculated the LD (r 2 measure) between the lead sQTL and all other SNPs within 1Mb using PLINK  [88]  (v1. 9 --ld-window-kb 1000 --ld-window 99999 --ld-window-r2 0). As a reference dataset for LD calculations, we used our high-quality genotypes in 399 samples. Then, for each gene, we plot the LD r 2 distribution between the lead sQTL and the lead QTL from the dataset being compared, respectively.
Genomic enrichment analysis of eQTLs and sQTLs We used GREGOR  [89]  to perform enrichment analysis of lead sQTLs and eQTLs in different genomic annotations using the following options R2THRESHOLD=0. 7, LDWINDOWSIZE=50000 (for sQTLs), LDWINDOWSIZE=100000 (for eQTLs), MIN_NEIGHBOR_NUM=200, and POPULATION=EUR. Human islet regulatory annotations were described previously  [11]. For genic annotations, we used Gencode  [22]  v34 GTF file. We used gtf2leafcutter. pl script to obtain all exon, intron, 5′ and 3′ coordinates from GTF file. 5′ and 3′ coordinates were extended to +3bp into exon and +2bp into intron. All annotations contain mutually exclusive genomic space. In parallel, we used GARFIELD  [90]  for enrichment analysis among genomic annotations. We used QTLtools beta approximated permutation p-values of lead eQTL and sQTLs from all genes and junctions tested in cis-eQTL and sQTL analysis.
Comparison with GTEx We obtained summary statistics data for eQTLs and sQTLs from 49 GTEx tissues  [17]. We first obtained eGenes and junctions with significant eQTLs and sQTLs, respectively, at FDR 0. 01, consistent with our significance threshold. For the resulting significant variant-phenotype associations for each of the 49 tissues, variant position and junction coordinates (for sQTLs) were lifted down from hg38 to hg19 using liftOver  [91]. For each GTEx tissue, we looked at the variant-phenotype (eGenes or Junctions) overlap with our islet eQTL and sQTLs, using nominal QTL variants. For example, for GTEx x eGene in j tissue, if any of the GTEx significant variants mapped any of our nominal eQTL variants for that x eGene, we considered that islet eQTL signal to be shared with that given j tissue.
To estimate T2D risk inflation across islet-selective sQTLs and eQTLs, we grouped sQTLs and eQTLs that were shared in ≤ 6 GTEx tissues (see the "Comparison with GTEx" section) as isletselective QTL effects. To assess the strength of T2D-risk inflation, we leveraged 1000 permutations of control-sets of sQTL-like variants as described above. on CLIP experiments on 82 unique RBPs (ENCODE and other CLIP datasets), which is implemented as Seqweaver  [46]. We performed in silico mutagenesis on both eQTL and sQTL credible set variants using both DeepSEA and Seqweaver models and obtained Disease Impact Scores (DIS) values from each respective model  [21]. Prior to in silico mutagenesis, the strand information was added to sQTL credible sets based on the orientation of the target gene. Credible sets were further stratified based on their location in the genome and DIS values from both models for each category of variants were shown as boxplots.
Acknowledgements We thank the  Genomics Facilities  from the  Center of Genomic Regulation and NIHR Imperial Biomedical Research Centre. Members of the  T2DSystems Consortium. Miriam Cnop: ULB Center for Diabetes Research, Medical Faculty, Université Libre De Bruxelles, Brussels, Belgium. Lena Eliasson: Lund University Diabetes Centre ;  Department of Clinical Sciences, Malmö, Lund University ; and  Clinical Research Centre, Skåne University Hospital, Malmö, Sweden. Jonathan Lou S. Esguerra: Lund University Diabetes Centre ;  Department of Clinical Sciences, Malmö, Lund University ; and  Clinical Research Centre, Skåne University Hospital, Malmö ,  Sweden. Décio L. Eizirik :  ULB Center for Diabetes Research, Medical Faculty ,  Université Libre De Bruxelles, Brussels, Belgium. Anna L. Gloyn: Oxford Centre for Diabetes, Endocrinology and Metabolism, Radcliffe Department of Medicine, University of Oxford, Oxford, UK.
Availability of data and materials All processed data files for the analyses described here, including all nominal sQTL and eQTL results, and fine-mapped eQTL and sQTL data, are available at Zenodo DOI: 10. 5281/zenodo. 6546807 [21]. Raw RNA-seq and genotyped SNP array data are available from the European Genome-phenome Archive (EGAS00001006440)  [70].