MicroPredict: predicting species-level taxonomic abundance of whole-shotgun metagenomic data using only 16S amplicon sequencing data.

BACKGROUND: The importance of the human microbiome in the analysis of various diseases is emerging. The two main methods used to profile the human microbiome are 16S rRNA gene sequencing (16S sequencing) and whole-genome shotgun sequencing (WGS). Owing to the full coverage of the genome in sequencing, WGS has multiple advantages over 16S sequencing, including higher taxonomic profiling resolution at the species-level and functional profiling analysis. However, 16S sequencing remains widely used because of its relatively low cost. Although WGS is the standard method for obtaining accurate species-level data, we found that 16S sequencing data contained rich information to predict high-resolution species-level abundances with reasonable accuracy. OBJECTIVE: In this study, we proposed MicroPredict, a method for accurately predicting WGS-comparable species-level abundance data using 16S taxonomic profile data. METHODS: We employed a mixed model using two key strategies: (1) modeling both sample- and species-specific information for predicting WGS abundances, and (2) accounting for the possible correlations among different species. RESULTS: We found that MicroPredict outperformed the other machine learning methods. CONCLUSION: We expect that our approach will help researchers accurately approximate the species-level abundances of microbiome profiles in datasets for which only cost-effective 16S sequencing has been applied.

Characterization of Th17 tissue-resident memory cells in non-inflamed intestinal tissue of Crohn's disease patients.

Crohn's disease (CD) is a chronic inflammatory disorder affecting the bowel wall. Tissue-resident memory T (Trm) cells are implicated in CD, yet their characteristics remain unclear. We aimed to investigate the transcriptional profiles and functional characteristics of Trm cells in the small bowel of CD and their interactions with immune cells. Seven patients with CD and four with ulcerative colitis as controls were included. Single-cell RNA sequencing and paired T cell receptor sequencing assessed T cell subsets and transcriptional signatures in lamina propria (LP) and submucosa/muscularis propria-enriched fractions (SM/MP) from small bowel tissue samples. We detected 58,123 T cells grouped into 16 populations, including the CD4+ Trm cells with a Th17 signature and CD8+ Trm clusters. In CD, CD4+ Trm cells with a Th17 signature, termed Th17 Trm, showed significantly increased proportions within both the LP and SM/MP areas. The Th17 Trm cluster demonstrated heightened expression of tissue-residency marker genes (ITGAE, ITGA1, and CXCR6) along with elevated levels of IL17A, IL22, CCR6, and CCL20. The clonal expansion of Th17 Trm cells in CD was accompanied by enhanced transmural dynamic potential, as indicated by significantly higher migration scores. CD-prominent Th17 Trm cells displayed an increased interferon gamma (IFNγ)-related signature possibly linked with STAT1 activation, inducing chemokines (i.e., CXCL10, CXCL8, and CXCL9) in myeloid cells. Our findings underscored the elevated Th17 Trm cells throughout the small bowel in CD, contributing to disease pathogenesis through IFNγ induction and subsequent chemokine production in myeloid cells.

PASTRY: achieving balanced power for detecting risk and protective minor alleles in meta-analysis of association studies with overlapping subjects.

BACKGROUND: Meta-analysis is a statistical method that combines the results of multiple studies to increase statistical power. When multiple studies participating in a meta-analysis utilize the same public dataset as controls, the summary statistics from these studies become correlated. To solve this challenge, Lin and Sullivan proposed a method to provide an optimal test statistic adjusted for the correlation. This method quickly became the standard practice. However, we identified an unexpected power asymmetry phenomenon in this standard framework. This can lead to unbalanced power for detecting protective minor alleles and risk minor alleles. RESULTS: We found that the power asymmetry of the current framework is mainly due to the errors in approximating the correlation term. We then developed a meta-analysis method based on an accurate correlation estimator, called PASTRY (A method to avoid Power ASymmeTRY). PASTRY outperformed the standard method on both simulated and real datasets in terms of the power symmetry. CONCLUSIONS: Our findings suggest that PASTRY can help to alleviate the power asymmetry problem. PASTRY is available at https://github.com/hanlab-SNU/PASTRY .

Transcriptomic Profiling and Cellular Composition of Creeping Fat in Crohn's disease.

BACKGROUND AND AIMS: Creeping fat [CF] is a poorly understood feature of Crohn's disease [CD], characterized by the wrapping of mesenteric adipose tissue [MAT] around the inflamed intestine. The aim of this study was to investigate the transcriptional profile and compositional features of CF. METHODS: We collected 59 MAT samples: 23 paired samples from patients with CD (CF [CD-CF] and MAT around the uninflamed intestine [CD-MAT]) and 13 MAT samples from non-CD patients [Con-MAT]. Differentially expressed gene [DEG], functional pathway, cell deconvolution, and gene co-expression network analyses were performed. RESULTS: By comparing three different MAT samples, we identified a total of 529 DEGs [|log2FoldChange| > 1.5; false discovery rate < 0.05]. Of these, 323 genes showed an incremental pattern from Con-MAT to CD-MAT, and to CD-CF, while 105 genes displayed a decremental pattern. Genes with an incremental pattern were related to immune cell responses, including B- and T-cell activation, while genes with a decremental pattern were involved in cell trafficking and migration. Cell deconvolution analysis revealed significant changes in cellular composition between the CD-CF and Con-MAT groups, with increased proportions of B-cells/plasma cells [p = 1.16 × 10-4], T-cells [p = 3.66 × 10-3], and mononuclear phagocytes [p = 3.53 × 10-2] in the CD-CF group. In contrast, only the B-cell/plasma cell component showed a significant increase [p = 1.62 × 10-2] in the CD-MAT group compared to Con-MAT. CONCLUSION: The distinct transcriptional profiles and altered cellular components of each MAT found in our study provide insight into the mechanisms behind CF and highlight its possible role in the pathogenesis of CD.

Building an optimal predictive model for imputing tissue-specific gene expression by combining genotype and whole-blood transcriptome data.

Accurate imputation of tissue-specific gene expression can be a powerful tool for understanding the biological mechanisms underlying human complex traits. Existing imputation methods can be grouped into two categories according to the types of predictors used. The first category uses genotype data, while the second category uses whole-blood expression data. Both data types can be easily collected from blood, avoiding invasive tissue biopsies. In this study, we attempted to build an optimal predictive model for imputing tissue-specific gene expression by combining the genotype and whole-blood expression data. We first evaluated the imputation performance of each standalone model (using genotype data [GEN model] and using whole-blood expression data [WBE model]) using their respective data types across 47 human tissues. The WBE model outperformed the GEN model in most tissues by a large gain. Then, we developed several combined models that leverage both types of predictors to further improve imputation performance. We tried various strategies, including utilizing a merged dataset of the two data types (MERGED models) and integrating the imputation outcomes of the two standalone models (inverse variance-weighted [IVW] models). We found that one of the MERGED models noticeably outperformed the standalone models. This model involved a fixed ratio between the two regularization penalty factors for the two predictor types so that the contribution of the whole-blood transcriptome is upweighted compared with the genotype. Our study suggests that one can improve the imputation of tissue-specific gene expression by combining the genotype and whole-blood expression, but the improvement can be largely dependent on the combination strategy chosen.

Tutorial: a statistical genetics guide to identifying HLA alleles driving complex disease.

The human leukocyte antigen (HLA) locus is associated with more complex diseases than any other locus in the human genome. In many diseases, HLA explains more heritability than all other known loci combined. In silico HLA imputation methods enable rapid and accurate estimation of HLA alleles in the millions of individuals that are already genotyped on microarrays. HLA imputation has been used to define causal variation in autoimmune diseases, such as type I diabetes, and in human immunodeficiency virus infection control. However, there are few guidelines on performing HLA imputation, association testing, and fine mapping. Here, we present a comprehensive tutorial to impute HLA alleles from genotype data. We provide detailed guidance on performing standard quality control measures for input genotyping data and describe options to impute HLA alleles and amino acids either locally or using the web-based Michigan Imputation Server, which hosts a multi-ancestry HLA imputation reference panel. We also offer best practice recommendations to conduct association tests to define the alleles, amino acids, and haplotypes that affect human traits. Along with the pipeline, we provide a step-by-step online guide with scripts and available software ( https://github.com/immunogenomics/HLA_analyses_tutorial ). This tutorial will be broadly applicable to large-scale genotyping data and will contribute to defining the role of HLA in human diseases across global populations.

Case-case genome-wide association analysis identifying genetic loci with divergent effects on Crohn's disease and ulcerative colitis.

Crohn's disease (CD) and ulcerative colitis (UC), two major subtypes of inflammatory bowel disease, show substantial differences in their clinical course and treatment response. To identify the genetic factors underlying the distinct characteristics of these two diseases, we performed a genome-wide association study (GWAS) between CD (n = 2359) and UC (n = 2175) in a Korean population, followed by replication in an independent sample of 772 CD and 619 UC cases. Two novel loci were identified with divergent effects on CD and UC: rs9842650 in CD200 and rs885026 in NCOR2. In addition, the seven established susceptibility loci [major histocompatibility complex (MHC), TNFSF15, OTUD3, USP12, IL23R, FCHSD2 and RIPK2] reached genome-wide significance. Of the nine loci, six (MHC, TNFSF15, OTUD3, USP12, IL23R and CD200) were replicated in the case-case GWAS of European populations. The proportion of variance explained in CD-UC status by polygenic risk score analysis was up to 22.6%. The area under the receiver-operating characteristic curve value was 0.74, suggesting acceptable discrimination between CD and UC. This CD-UC GWAS provides new insights into genetic differences between the two diseases with similar symptoms and might be useful in improving their diagnosis and treatment.

Exploration of shared features of B cell receptor and T cell receptor repertoires reveals distinct clonotype clusters.

Although B cells and T cells are integral players of the adaptive immune system and act in co-dependent ways to orchestrate immune responses, existing methods to study the immune repertoire have largely focused on separate analyses of B cell receptor (BCR) and T cell receptor (TCR) repertoires. Based on our hypothesis that the shared history of immune exposures and the shared cellular machinery for recombination result in similarities between BCR and TCR repertoires in an individual, we examine any commonalities and interrelationships between BCR and TCR repertoires. We find that the BCR and TCR repertoires have covarying clonal architecture and diversity, and that the pattern of correlations appears to be altered in immune-mediated diseases. Furthermore, hierarchical clustering of public B and T cell clonotypes in both health and disease based on correlation of clonal proportion revealed distinct clusters of B and T cell clonotypes that exhibit increased sequence similarity, share motifs, and have distinct amino acid characteristics. Our findings point to common principles governing memory formation, recombination, and clonal expansion to antigens in B and T cells within an individual. A significant proportion of public BCR and TCR repertoire can be clustered into nonoverlapping and correlated clusters, suggesting a novel way of grouping B and T cell clonotypes.

FastRNA: An efficient solution for PCA of single-cell RNA-sequencing data based on a batch-accounting count model.

Almost always, the analysis of single-cell RNA-sequencing (scRNA-seq) data begins with the generation of the low dimensional embedding of the data by principal-component analysis (PCA). Because scRNA-seq data are count data, log transformation is routinely applied to correct skewness prior to PCA, which is often argued to have added bias to data. Alternatively, studies have proposed methods that directly assume a count model and use approximately normally distributed count residuals for PCA. Despite their theoretical advantage of directly modeling count data, these methods are extremely slow for large datasets. In fact, when the data size grows, even the standard log normalization becomes inefficient. Here, we present FastRNA, a highly efficient solution for PCA of scRNA-seq data based on a count model accounting for both batches and cell size factors. Although we assume the same general count model as previous methods, our method uses two orders of magnitude less time and memory than the other count-based methods and an order of magnitude less time and memory than the standard log normalization. This achievement results from our unique algebraic optimization that completely avoids the formation of the large dense residual matrix in memory. In addition, our method enjoys a benefit that the batch effects are eliminated from data prior to PCA. Generating a batch-accounted PC of an atlas-scale dataset with 2 million cells takes less than a minute and 1 GB memory with our method.

A theory-based practical solution to correct for sex-differential participation bias.

Most genomic cohorts are retrospective where the exposures and outcomes are predetermined prior to sample collection. Therefore, a spurious association between an exposure and an outcome can arise if both variables affect study participation. Such concerns were raised in previous studies questioning the representativeness of the UK Biobank. Recently, a genome-wide association study (GWAS) on biological sex found many autosomal hits and non-negligible autosomal heritability which the authors attribute to selection bias. In this study, we propose a simple and a practical method that can overcome sex-driven selection bias based on theoretical analysis and simulations.

Identification of shared loci associated with both Crohn's disease and leprosy in East Asians.

Genome-wide association studies (GWAS) of Crohn's disease (CD) in European and leprosy in Chinese population have shown that CD and leprosy share genetic risk loci. As these shared loci were identified through cross-comparisons across different ethnic populations, we hypothesized that meta-analysis of GWAS on CD and leprosy in East Asian populations would increase power to identify additional shared loci. We performed a cross-disease meta-analysis of GWAS data from CD (1621 cases and 4419 controls) and leprosy (2901 cases 3801 controls) followed by replication in additional datasets comprising 738 CD cases and 488 controls and 842 leprosy cases and 925 controls. We identified one novel locus at 7p22.3, rs77992257 in intron 2 of ADAP1, shared between CD and leprosy with genome-wide significance (P = 3.80 × 10-11) and confirmed 10 previously established loci in both diseases: IL23R, IL18RAP, IL12B, RIPK2, TNFSF15, ZNF365-EGR2, CCDC88B, LACC1, IL27, NOD2. Phenotype variance explained by the polygenic risk scores derived from Chinese leprosy data explained up to 5.28% of variance of Korean CD, supporting similar genetic structures between the two diseases. Although CD and leprosy shared a substantial number of genetic susceptibility loci in East Asians, the majority of shared susceptibility loci showed allelic effects in the opposite direction. Investigation of the genetic correlation using cross-trait linkage disequilibrium score regression also showed a negative genetic correlation between CD and leprosy (rg [SE] = -0.40[0.13], P = 2.6 × 10-3). These observations implicate the possibility that CD might be caused by hyper-sensitive reactions toward pathogenic stimuli.

Exploration of errors in variance caused by using the first-order approximation in Mendelian randomization.

Mendelian randomization (MR) uses genetic variation as a natural experiment to investigate the causal effects of modifiable risk factors (exposures) on outcomes. Two-sample Mendelian randomization (2SMR) is widely used to measure causal effects between exposures and outcomes via genome-wide association studies. 2SMR can increase statistical power by utilizing summary statistics from large consortia such as the UK Biobank. However, the first-order term approximation of standard error is commonly used when applying 2SMR. This approximation can underestimate the variance of causal effects in MR, which can lead to an increased false-positive rate. An alternative is to use the second-order approximation of the standard error, which can considerably correct for the deviation of the first-order approximation. In this study, we simulated MR to show the degree to which the first-order approximation underestimates the variance. We show that depending on the specific situation, the first-order approximation can underestimate the variance almost by half when compared to the true variance, whereas the second-order approximation is robust and accurate.

MarcoPolo: a method to discover differentially expressed genes in single-cell RNA-seq data without depending on prior clustering.

The standard analysis pipeline for single-cell RNA-seq data consists of sequential steps initiated by clustering the cells. An innate limitation of this pipeline is that an imperfect clustering result can irreversibly affect the succeeding steps. For example, there can be cell types not well distinguished by clustering because they largely share the global structure, such as the anterior primitive streak and mid primitive streak cells. If one searches differentially expressed genes (DEGs) solely based on clustering, marker genes for distinguishing these types will be missed. Moreover, clustering depends on many parameters and can often be subjective to manual decisions. To overcome these limitations, we propose MarcoPolo, a method that identifies informative DEGs independently of prior clustering. MarcoPolo sorts out genes by evaluating if the distributions are bimodal, if similar expression patterns are observed in other genes, and if the expressing cells are proximal in a low-dimensional space. Using real datasets with FACS-purified cell labels, we demonstrate that MarcoPolo recovers marker genes better than competing methods. Notably, MarcoPolo finds key genes that can distinguish cell types that are not distinguishable by the standard clustering. MarcoPolo is built in a convenient software package that provides analysis results in an HTML file.

Dynamic increase of M2 macrophages is associated with disease progression of colorectal cancers following cetuximab-based treatment.

We aimed to investigate the dynamic changes of gene expression profiles and immune microenvironment linked to resistance to cetuximab-based treatments in patients with metastatic colorectal cancer (mCRC). A total of 106 patients with RAS-wild type mCRC who were treated with cetuximab-based treatments were included as the study population. RNA-sequencing and multiplexed immunohistochemistry were performed using paired or unpaired pre-treatment and post-treatment tumor tissues. Differentially expressed gene analysis of paired pre-treatment and post-treatment tumor tissues that develop acquired resistance (AR) identified the AR signature. Gene ontology analysis of the AR signature indicated enrichment of immune-related pathway genes. Among the immune subsets whose abundance was estimated by CIBERSORT, M2 macrophages showed the most prominent positive correlation with the expression of the AR signature. Among the post-treatment samples, progressive disease (PD) tumors showed a significantly higher abundance of M2 macrophages compared to non-PD tumors. These findings were validated by multiplexed immunohistochemistry analysis: the density of CD68+CD206+ M2 macrophages significantly increased at the time of PD following cetuximab-based treatment, whereas it did not consistently change in the tumor pairs of non-PD. In conclusion, a dynamic increase of M2 macrophages is associated with disease progression during cetuximab-based treatment of mCRCs. Targeting M2 macrophages is a promising immunotherapeutic strategy in this clinical context.

Phenome-wide association study of the major histocompatibility complex region in the Korean population identifies novel association signals.

Human leukocyte antigen (HLA) gene variants in the major histocompatibility complex (MHC) region are associated with numerous complex human diseases and quantitative traits. Previous phenome-wide association studies (PheWAS) for this region demonstrated that HLA association patterns to the phenome have both population-specific and population-shared components. We performed MHC PheWAS in the Korean population by analyzing associations between phenotypes and genetic variants in the MHC region using the Korea Biobank Array project data samples from the Korean Genome and Epidemiology Study cohorts. Using this single-population dataset, we curated and analyzed 82 phenotypes for 125 673 Korean individuals after imputing HLA using CookHLA, a recently developed imputation framework. More than one-third of these phenotypes showed significant associations, confirming 56 known associations and discovering 13 novel association signals that were not reported previously. In addition, we analyzed heritability explained by the variants in the MHC region and genetic correlations among phenotypes based on the MHC variants.

Expression of neurofibromin 1 in colorectal cancer and cetuximab resistance.

Neurofibromin 1 (NF1) is a tumor suppressor that has been previously reported to regulate RAS­MAPK signaling. The present study investigated the possible relationship between NF1 expression and anti­EGFR antibody (cetuximab) sensitivity in colorectal cancer cell lines. In addition, primary or metastatic colorectal cancer samples from patients treated with cetuximab were assessed for the association of cetuximab sensitivity. The quantities of the NF1 transcript, NF1­related pathway enrichment and NF1 mutation profile were measured and investigated using RNA sequencing and targeted DNA sequencing. Based on growth inhibition and colony formation assay results, cell lines were designated to be cetuximab­sensitive (NCI­H508 and Caco2) or cetuximab­resistant (KM12C and SM480). Western blotting revealed NF1 was highly expressed in cetuximab­sensitive cell lines whilst there was little expression in their cetuximab­resistant counterparts. Knocking down NF1 expression using small interfering RNA in the cetuximab­sensitive cell lines enhanced the phosphorylation of MEK and ERK according to western blotting. NF1 knockdown also reduced apoptosis, as observed by the decreased number of apoptotic bodies by DAPI nuclear staining and reduced cleavage of caspase and poly­(ADP ribose) polymerase. NF1 overexpression by transfection with GTPase­activating protein­related domain subunit rendered the cetuximab­resistant cell lines, KM12C and SW480, more susceptible to cetuximab­induced apoptosis. RNA sequencing of 111 RAS and BRAFV600 wild­type tumor samples collected from cetuximab­treated patients with metastatic colorectal cancer revealed that the pre­treatment NF1 expression levels were not associated with the cetuximab response. However, tumor samples obtained after cetuximab treatment displayed slightly lower NF1 transcript levels compared with those in the pre­treatment samples, suggesting that exposure to the anti­EGFR antibody may be associated with reduced NF1 expression levels. Next­generation sequencing revealed that the frequency of inactivating mutations in NF1 were rare (1.8%) in patients with colorectal cancer and were not associated with the protein expression levels of NF1 except for in a small number of cases (0.5%), where the biallelic inactivation of NF1 was observed. To conclude, the present study showed that modification of NF1 expression can affect sensitivity to cetuximab in colorectal cancer cell lines, though a limitation exists in terms of its potential application as a biomarker for RAS and BRAFV600 wild­type tumors.

A high-resolution HLA reference panel capturing global population diversity enables multi-ancestry fine-mapping in HIV host response.

Fine-mapping to plausible causal variation may be more effective in multi-ancestry cohorts, particularly in the MHC, which has population-specific structure. To enable such studies, we constructed a large (n = 21,546) HLA reference panel spanning five global populations based on whole-genome sequences. Despite population-specific long-range haplotypes, we demonstrated accurate imputation at G-group resolution (94.2%, 93.7%, 97.8% and 93.7% in admixed African (AA), East Asian (EAS), European (EUR) and Latino (LAT) populations). Applying HLA imputation to genome-wide association study data for HIV-1 viral load in three populations (EUR, AA and LAT), we obviated effects of previously reported associations from population-specific HIV studies and discovered a novel association at position 156 in HLA-B. We pinpointed the MHC association to three amino acid positions (97, 67 and 156) marking three consecutive pockets (C, B and D) within the HLA-B peptide-binding groove, explaining 12.9% of trait variance.

Rare variants regulate expression of nearby individual genes in multiple tissues.

The rapid decrease in sequencing cost has enabled genetic studies to discover rare variants associated with complex diseases and traits. Once this association is identified, the next step is to understand the genetic mechanism of rare variants on how the variants influence diseases. Similar to the hypothesis of common variants, rare variants may affect diseases by regulating gene expression, and recently, several studies have identified the effects of rare variants on gene expression using heritability and expression outlier analyses. However, identifying individual genes whose expression is regulated by rare variants has been challenging due to the relatively small sample size of expression quantitative trait loci studies and statistical approaches not optimized to detect the effects of rare variants. In this study, we analyze whole-genome sequencing and RNA-seq data of 681 European individuals collected for the Genotype-Tissue Expression (GTEx) project (v8) to identify individual genes in 49 human tissues whose expression is regulated by rare variants. To improve statistical power, we develop an approach based on a likelihood ratio test that combines effects of multiple rare variants in a nonlinear manner and has higher power than previous approaches. Using GTEx data, we identify many genes regulated by rare variants, and some of them are only regulated by rare variants and not by common variants. We also find that genes regulated by rare variants are enriched for expression outliers and disease-causing genes. These results suggest the regulatory effects of rare variants, which would be important in interpreting associations of rare variants with complex traits.

A systematic analysis of protein-altering exonic variants in chronic obstructive pulmonary disease.

Genome-wide association studies (GWASs) have identified regions associated with chronic obstructive pulmonary disease (COPD). GWASs of other diseases have shown an approximately 10-fold overrepresentation of nonsynonymous variants, despite limited exonic coverage on genotyping arrays. We hypothesized that a large-scale analysis of coding variants could discover novel genetic associations with COPD, including rare variants with large effect sizes. We performed a meta-analysis of exome arrays from 218,399 controls and 33,851 moderate-to-severe COPD cases. All exome-wide significant associations were present in regions previously identified by GWAS. We did not identify any novel rare coding variants with large effect sizes. Within GWAS regions on chromosomes 5q, 6p, and 15q, four coding variants were conditionally significant (P < 0.00015) when adjusting for lead GWAS single-nucleotide polymorphisms A common gasdermin B (GSDMB) splice variant (rs11078928) previously associated with a decreased risk for asthma was nominally associated with a decreased risk for COPD [minor allele frequency (MAF) = 0.46, P = 1.8e-4]. Two stop variants in coiled-coil α-helical rod protein 1 (CCHCR1), a gene involved in regulating cell proliferation, were associated with COPD (both P < 0.0001). The SERPINA1 Z allele was associated with a random-effects odds ratio of 1.43 for COPD (95% confidence interval = 1.17-1.74), though with marked heterogeneity across studies. Overall, COPD-associated exonic variants were identified in genes involved in DNA methylation, cell-matrix interactions, cell proliferation, and cell death. In conclusion, we performed the largest exome array meta-analysis of COPD to date and identified potential functional coding variants. Future studies are needed to identify rarer variants and further define the role of coding variants in COPD pathogenesis.