dataesr/scientific-paragraphs-categorization huggingface dataset:
Data analysis Amplification products were scored using the GeneMapper software version 3. 7 (Applied Biosystems) and only unambiguous amplification products were considered for data analysis. Genetic diversity parameters were calculated for each locus and each sampling location using the GenAlEx software version 6. 502  [88]. Expected heterozygosity (He) was calculated using the formula: He ¼ 1 - X p i 2 where p i is the frequency of each allele. The fixation index (F) was calculated as: F ¼ 1 - H 0 He where Ho is observed heterozygosity and He is expected heterozygosity  [89]. F-statistics analysis assessing genetic differentiation (Fst), genetic identity, number of migrants (Nm)  [90]  and non-hierarchical analysis of molecular variance (AMOVA) for estimating genetic differentiation within and among locations were performed with the same software.
Allelic richness was calculated using the SPA-GeDi software version 1. 5 (http://ebe. ulb. ac. be/ebe/SPA-GeDi. html [91] ;). Consecutively to K determination (see below), successive hierarchical AMOVA analyses were carried out with K = 2 and K = 3. The Mantel permutation test was used for assessing the correlation between genetic and geographic distances between sampling sites  [92, 93]. Two Principal Coordinates Analyses (PCoA) enabling the visualization of genetic variation distribution across individuals and sampling sites, respectively, were performed using GenAlEx. The STRUCTURE software version 2. 3. 4  [94]  was used for the determination of the most probable number of clusters for population structure (K value). Using the admixture model, eight simulations were performed for each inferred K value, with a running length composed of 300,000 burn-in periods and 50,000 Markov chain Monte Carlo (MCMC) replicates.
The output from this analysis was then used as input in the Structure HARV ESTER online program version 0. 6. 94 (http://taylor0. biology. ucla. edu/structureHarvester/) to determine the optimal value of K using the ΔK method of Evanno et al. [95]  and allowing for different estimates of K in accordance with Janes et al  [50]. Based on the resulting values of K, a clustering analysis of the studied sampling sites was performed and graphical output was generated using CLUMPAK's main pipeline (http://clumpak. tau. ac. il  [96] ; ). In order to further assess genetic clustering, a UPGMA tree based on Fst values using 1000 bootstrap replications was constructed using the POPTREE2 software  [97]. MJK, KA, SS, CM, LZ and TB processed and analyzed data. All authors contributed to writing and revising the manuscript.
Acknowledgements The Authors wish to thank  Emira Cherif ,  Karina Castillo  and  Frédérique Aberlenc  ( IRD, UMR DIADE, Montpellier, France ) for sharing their unpublished microsatellite markers data for transferability tests, and  Dr. Jaume Pellicer  ( Comparative Plant & Fungal Biology, Jodrell Laboratory, Royal Botanical Gardens, Kew, UK ), for allowing us to use the palm flow cytometry data.
Funding The work described in this article was funded through travel grants to MJK and KA under the framework of the  MooSciTIC  project granted to EJ by  Agropolis Fondation  ( ID 1501-011 , " Investissements d'Avenir" Program -Labex Agro :  ANR-10-LABX-0001-01 ). Additional funding was provided by the  Sud Expert Plantes -Développement Durable  ( SEP2D ) programme to KA (GenPhyB project, ID  AAP3-64 ). MJK is the recipient of a PhD fellowship from the  French Embassy in Benin. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Operational Taxonomic Units (OTUs) The sequencing data was processed into OTUs using Mothur (Version v. 1. 35. 1) following the MiSeq standard operational procedure  (Schloss et al. , 2009) , which are described in detail by  Raju et al. (2018). In summary, the sequences were filtered on quality and the read pairs were assembled. The assembled reads were aligned to the SILVA 16S rRNA gene database (Version V119) and clustered into OTUs with a >98% similarity cut-off to enable species' resolution. In total 6536 OTUs were identified, whereas approximately 1000 of the most abundant OTUs were included in the statistical analyses, depending on the specific sample set. The data were normalized with a subsampling threshold of 500. Samples below this threshold were excluded from the analyses. We rarefied data to normalize the sequenced depth variability.
After rarefaction, we calculated the Shannon index to assess richness  (Wagner et al. , 2018)  and inverse Simpson's index to assess richness and evenness of the microbiota for alpha diversity  (Simpson, 1949). The Bray-Curtis dissimilarity index (calculated from counts) was used for assessing compositional differences of the saliva microbiota between eating habits and breakfast and dinner patterns. The Bray-Curtis index of 0 means that two samples have the same composition, while the Bray-Curtis index of 1 means that two samples do not share any species together. Data for the abundance analyses were pre-filtered to include the OTUs with a minimum of 20 sequences.
, 2008). We adjusted for gender  (Raju et al. , 2019) , age  (Yatsunenko et al. , 2012; Lif Holgerson et al. , 2015) , and BMI  (Piombino et al. , 2014) , which have found to associate with oral and saliva microbiota. Furthermore, we adjusted for sequencing depth since it affects the data quality even after rarefaction. Additionally, we adjusted for eating habits when assessing breakfast and dinner patters. Interactions between gender and each of the dietary variables (eating habits and breakfast and dinner patterns) were tested. The comparisons of composition between eating habits and breakfast and dinner patterns were conducted by a permutational analysis of variance (PERMANOVA) for the matrix of Bray-Curtis index and illustrated by Principal Coordinates Analysis (PCoA) plots. For the analyses, we used adonis and betadisper functions in Phyloseq and Vegan R-packages.
The DESeq2 R-package (Version 3. 2. 4) tested the differential abundance of OTUs using a negative binomial model. In addition, the multiple comparison method was used to compare between two of the three eating habit groups: FV avoiders, healthy and unhealthy. We also tested the differential abundances of OTUs between irregular breakfast vs. regular breakfast patterns and irregular dinner vs. regular dinner patterns. We used false discovery rate (FDR) adjustment for p-value. All analyses were performed in R and SPSS version 24. 0 (IBM Corp. , Armonk, NY, United States). We illustrated the differential abundances OTUs with volcano plots for: (a) FV avoiders vs. healthy, (b) unhealthy vs. healthy, (c) FV avoiders vs. unhealthy, (d) irregular breakfast vs. regular breakfast patterns, and (e) irregular dinner vs.
Patient cohort We sequenced 1,318 individuals recruited as part of the PID domain of the United Kingdom NIHR BioResource -Rare Diseases program (NBR-RD; Extended Data Fig. 2 ; Supplementary Methods). The cohort comprised of both sporadic and familial PID patients (N=974) and family members. Of the patients, 886 were index cases who fell into one of the diagnostic categories of the European Society for Immunodeficiencies (ESID) registry diagnostic criteria (Fig. 1a ; Extended Data Table  1 ). This cohort represents a third of CVID and half of CID patients registered in the UK  9. Clinical phenotypes were dominated by adult-onset sporadic AD-PID: all had recurrent infections, 28% had autoimmunity, and 8% had malignancy (Fig. 1a-b , Extended Data Table  2 ), mirroring the UK national PID registry  6.
3a ) that has no coverage in the ExAC database (http://exac. broadinstitute. org). Two unaffected first-degree relatives were heterozygous for the frameshift variant, and two for the promoter deletion (Fig. 3b ), confirming compound heterozygosity in the patient. Western blotting demonstrated complete absence of ARPC1B and raised ARPC1A in platelets  22  (Fig. 3c ). ARPC1B mRNA was almost absent from mononuclear cells in the patient and was reduced in a clinically unaffected sister carrying the frameshift mutation (Supplementary Note 1). An allele specific expression assay demonstrated that the promoter deletion essentially abolished mRNA expression (Supplementary Note 1). ARPC1B is part of the Arp2/3 complex necessary for normal actin assembly in immune cells  23  , and monocyte-derived macrophages from the patient had an absence of podosomes, phenocopying deficiency of the Arp2/3 regulator WASp (Fig.
GWAS data allows identification of candidate monogenic PID genes and disease-modifying variants To investigate whether loci identified by GWAS of AD-PID and other IMD might be used to prioritize novel candidate monogenic PID genes, we used the data-driven pcHiC omnibus gene score (COGS) approach  21  (Methods, Supplementary Table  4 ). We selected six protein-coding genes with above average prioritisation scores in one or more diseases (Fig. 4b ), and identified a single protein truncating variant in each of ETS1, SOCS1 and PTPN2 genes, all occurring exclusively in PID patients. SOCS1 and PTPN2 variants were analysed further. SOCS1 limits phosphorylation of targets including STAT1, and is a key regulator of IFN-γ signalling  27. The patient with a heterozygous de-novo protein-truncating SOCS1 variant (p.
In summary, we show that cohort-based WGS in PID is a powerful approach to provide diagnosis of known genetic defects, and discover new coding and non-coding variants associated with disease (comparison of WGS with other methodologies; Supplementary Note 5). Improved analysis methodology and better integration of parallel datasets, such as GWAS and cell surface or metabolic immunophenotyping, will allow further exploration of the non-coding space, enhancing diagnostic yield. Such an approach promises to transform our understanding of genotype-phenotype relationships in PID and related immune-mediated conditions, and could redefine the clinical boundaries of immunodeficiency, add to our understanding of human immunology, and ultimately improve patient outcomes. Figure Legends.
PID cohort The PID patients and their family members were recruited by specialists in clinical immunology across 26 hospitals in the UK, and one each from the Netherlands, France and Germany. The recruitment criteria were intentionally broad, and included the following: clinical diagnosis of common variable immunodeficiency disorder (CVID) according to internationally established criteria (Extended Data Table  1 ); extreme autoimmunity; or recurrent and/or unusual severe infections suggestive of defective innate or cell-mediated immunity. Patients with known secondary immunodeficiencies caused by cancer or HIV infection were excluded. Although screening for more common and obvious genetic causes of PID prior to enrolment into this WGS study was encouraged, it was not a requirement. Consequently, a minority of patients (16%) had some prior genetic testing, from single gene Sanger sequencing or MLPA to a gene panel screen.
Paediatric and familial cases were less frequent in our cohort, in part reflecting that genetic testing is more frequently performed in more severe cases: 31% of paediatric onset cases had prior genetic testing compared to 10% of adult index cases (Extended Data Fig. 2 ). To expedite recruitment a minimal clinical dataset was required for enrolment, though more detail was often provided. There was a large variety in patients' phenotypes, from simple "chest infections" to complex syndromic features, and the collected phenotypic data of the sequenced individuals ranged from assigned disease category only to detailed clinical synopsis and immunophenotyping data. The clinical subsets used to subdivide PID patients were based on ESID definitions, as shown in Extended Data Table  1. The final PID cohort that we sequenced comprised of 886 index cases, 88 affected relatives, and 344 family members unaffected at the time of recruitment.
WGS data processing Details of DNA sample processing, whole genome sequencing, data processing pipeline, quality checks, alignment and variant calling, ancestry and relatedness estimation, variant normalisation and annotation, large deletion calling and filtering, and allele frequency calculations, are described in  38. Briefly, DNA or whole blood EDTA samples were processed and quality checked according to standard laboratory practices and shipped on dry ice to the sequencing provider (Illumina Inc, Great Chesterford, UK). Illumina Inc performed further QC array genotyping, before fragmenting the samples to 450bp fragments and processing with the Illumina TruSeq DNA PCR-Free Sample Preparation kit (Illumina Inc. , San Diego, CA, USA). Over the three-year duration of the sequencing phase of the project, different instruments and read lengths were used: for each sample, either 100bp reads on three HiSeq2500 lanes; or 125bp reads on two HiSeq2500 lanes; or 150bp reads on a single HiSeq X lane.
Each delivered genome had a minimum 15X coverage over at least 95% of the reference autosomes. Illumina performed the alignment to GRCh37 genome build and SNV/InDel calling using their Isaac software, while large deletions were called with their Manta and Canvas algorithms. The WGS data files were received at the University of Cambridge High Performance Computing Service (HPC) for further QC and processing by our Pipeline team. For each sample, we estimated the sex karyotype and computed pair-wise kinship coefficients (full methods described in  47  ), which allowed us to identify sample swaps and unintended duplicates, assign ethnicities, generate networks of closely related individuals (sometimes undeclared relatives from across different disease domains) and a maximal unrelated sample set (for the purposes of allele frequency estimation and control dataset in case-control analyses).
Variants in the gVCF files were normalised and loaded into an HBase database, where Overall Pass Rate (OPR) was computed within each of the three read length batches, and the lowest of these OPR values (minOPR) assigned to each variant. The rare variant analyses presented here are based on SNVs/InDels with minOPR>0. 98. Variants were annotated with Sequence Ontology terms according to their predicted consequences, their frequencies in other genomic databases (gnomAD, UK10K, 1000 Genomes), if they have been associated with a disease according to the HGMD Pro database, and internal metrics (AN, AC, AF, OPR). Large deletions (those >50bp in length, defined by Illumina) were merged and analysed collectively, as described in  38. Briefly, sample-level calls by the two algorithms, Manta (which uses read and mate-pair alignment information) and Canvas (which relies on read depth and is optimised for calls >1kb in length), were combined according to a set of rules  38  to generate a high quality set for each sample (and a large number across the project was visually inspected to ensure reasonably high specificity).
Diagnostic reporting We screened all genes in the International Union of Immunological Societies (IUIS) 2015 classification for previously reported or likely pathogenic variants. SNVs and small InDels were filtered based on the following criteria: OPR>0. 95; having a protein-truncating consequence, gnomAD AF<0. 001 and internal AF<0. 01; or present in the HGMD Pro database as DM variant. Large deletions called by both Canvas and Manta algorithms, passing standard Illumina quality filters, overlapping at least one exon, and classified as rare by the SVH method were included in the analysis. In order to aid variant interpretation and consistency in reporting, phenotypes were translated into Human Phenotype Ontology (HPO) terms as much as possible. Multi-Disciplinary Team (MDT) then reviewed each variant for evidence of pathogenicity and contribution to the phenotype, and classified them according to the American College of Medical Genetics (ACMG) guidelines  11.
008 one-sided Wilcoxon rank-sum test). BeviMed's sensitivity in prioritizing genes as causal, even if variants exist in only a few cases, is demonstrated by the observation that of the 8 IUIS-defined causal PID genes in the top 50 (all with a BeviMed PPA>0. 2), 3 are driven by 2 or 3 cases, while 5 have between 4 and 16. As allele frequency datasets for non-Europeans are much smaller than for Europeans, potential false positives may be induced by the unintentional inclusion of rare variants observed only in non-European populations  41. Furthermore, whilst the BeviMed analysis was restricted to the set of cases and controls carefully filtered to minimise relatedness, it remains possible that some associations could be false positives due to residual population stratification. We addressed this by flagging variants whose prioritisation was dependent upon cases with non-European ancestry.
When trained on our dataset, SKAT and BeviMed both identified NKFB1 as the gene with the strongest association signal, but BeviMed placed 8 IUIS 2017 PID genes in the top 50 results whilst SKAT placed 5, and ARPC1B was ranked 38th by BeviMed and 289th by SKAT (out of a total of 31,350 tested genes), consistent with the superiority of BeviMed over SKAT and related methods demonstrated in Greene et al. 1.
Phasing of DOCK8 variants In order to confirm the phase of two variants detected in the DOCK8 gene of a single individual, chr9:g. 306626-358548del and chr9:463519G>A, long read sequencing was performed using the Oxford Nanopore Technologies PromethION platform. The DNA sample was prepared using the 1D ligation library prep kit (SQK-LSK109), and genomic libraries were sequenced using a R. 9. 4. 1 PromethION flowcell. Raw signal data in FAST5 format was base called using Guppy (v2. 3. 5) to generate sequences in FASTQ format, which were then aligned against the GRCh37/hg19 human reference genome using minimap2 (v2. 2). Average coverage was 14x and median read length was 4,558 ± 4,007. A high quality set of heterozygous genotypes for the sample was created by using only variants from the short read Illumina WGS data with a phred score of <20 (probability of correct genotype > 0.
99). Haplotyping was then performed with Whatshap (v0. 14. 1) by using the long Nanopore reads to bridge across the informative genotypes from the short read data (https://whatshap. readthedocs. io/en/latest/index. html). We obtained a single high confidence haplotype block spanning the large deletion and the rare missense variant and showing that they were in trans (Extended data Fig. 5 ).
AD-PID GWAS GWAS was performed both on the whole PID cohort (N cases = 886) and on a subset comprising AD-PID cases (N cases = 733); the results of the AD-PID analysis were less noisy, and had increased power to detect statistical associations despite a reduced sample size (Extended Data Fig. 6 ). We used 9,225 unrelated samples from non-PID NBR-RD cohorts as controls. Variants selected from a merged VCF file were filtered to include bi-allelic SNPs with overall MAF>=0. 05 and minOPR=1 (100% pass rate across all WGS data for over 13,000 NBR participants). We ran PLINK logistic association test under an additive model. We adjusted for read length to guard against technical differences in genotype calls across the samples sequenced using 100bp, 125bp and 150bp reads, as Illumina chemistries changed throughout the duration of the project.
We also used sex and first 10 principal components from the ethnicity analysis as covariates, to mitigate against any population stratification effects. After filtering out SNPs with HWE p<10 -6 , we were left with the total of 4,993,945 analysed SNPs. There was minimal genomic inflation of the test statistic (lambda = 1. 022), suggesting population substructure and sample relatedness had been appropriately accounted for. Linear mixed model (LMM) analysis, as implemented in the BOLT-LMM package  43  , is an alternative method of association testing correcting for population stratification. It was used to confirm the observed associations (Extended Data Table  3 ). After genomic control correction  44  the only genome-wide significant (p<5x10 -8 ) signal was at the MHC locus, with several suggestive (p<1x10 -5 ) signals (Extended Data Fig.
6 ). We repeated the analysis with more relaxed SNP filtering criteria using 0. 005 < MAF < 0. 05 and minOPR>0. 95 (Extended Data Fig. 6 ). The only additional signal identified were the three TNFRSF13B variants shown in Supplementary Note 3. We obtained summary statistics data from the Li et al. CVID Immunochip case-control study  8  and, after further genomic control correction (lambda = 1. 039), performed a fixed effects meta-analysis on 95,417 variants shared with our AD-PID GWAS. Genome-wide significant (p<5x10 -8 ) signals were seen at the MHC and 16p13. 13 loci, with several suggestive (p<1x10 -5 ) signals (Extended Data Table  3 ). After metaanalysis, we conditioned on the lead SNP in each of the genome-wide and suggestive loci by including it as an additional covariate in the logistic regression model in PLINK, to determine if the signal was driven by single or multiple hits at those loci.
MHC locus analyses We imputed classical HLA alleles using the method implemented in the SNP2HLA v1. 0. 3 package  45  , which uses Beagle v3. 0. 4 for imputation and the HapMap CEU reference panel. We imputed allele dosages and best-guess genotypes of 2-digit and 4-digit classical HLA alleles, as well as amino acids of the MHC locus genes HLA-A, HLA-B, HLA-C, HLA-DRB1, HLA-DQA1 and HLA-DQB1. We tested the association of both allele dosages and genotypes using the logistic regression implemented in PLINK, and obtained similar results. We then used the best-guess genotypes to perform the conditional analysis (see above), since conditioning is not implemented in PLINK in a model with allele dosages. We repeated the conditional analyses as described above. The results of the sequential conditioning on the two lead classical alleles and amino acids within the Class I and Class II regions are shown in Extended Data Fig.
In order to compute the PPV for a given threshold on the posterior probability of association (PPA), we computed PPAs for 10,000 datasets without permutation of the case/control labels and 10,000 further datasets with a permutation of the case/control labels. We then sampled 1,000 PPAs from the permuted set and 10 PPAs from the non-permuted set to compute the PPV obtained when the PP threshold was set to achieve 100% power. The mean over 2,000 repetitions of this procedure is shown on the y-axis. The xaxis shows the number of cases in a hypothetical cohort. As the number of cases increases from 1,000 to 5,000, the PPV increases above 87. 5% irrespective of the proportion of cases with the same genetic aetiology. This demonstrates the utility of expanding the size of the PID case collection for detecting even very rare aetiologies resulting in the same broad phenotype as cases with different aetiologies.
Extended Data Table 2 -Description of the NIHR BioResource -Primary Immunodeficiency cohort. High-level clinical description and relevant clinical features were provided by recruiting clinicians. Index cases are patients recruited as sporadic cases or probands in pedigrees, and determined to be genetically unrelated by pairwise comparisons of common SNP genotypes in the WGS data. Numbers in brackets refer to the percentage of index cases in each category. Total number of patients is the sum of index cases and any affected relatives sequenced in this study.