The first category of summary statistics is Asian ancestry GWAS. The datasets are from the Biobank of Japan (BBJ)¹ (Ishigaki et al., 2020). We focus on the CLD-related phenotype that contain allele information and variant ID and that contain effect size and its standard error. With the two criteria, we obtained two GWAS summary statistics: cirrhosis (n = 212,453, prevalence = 1.03%) and HCC (n = 197,611, prevalence = 0.94%). Here, cirrhosis and HCC in BBJ were adjusted for age, sex, and top five genotype PCs (Ishigaki et al., 2020). The details of the two GWAS data are provided in Supplementary Table 1. Based on Asian ancestry from the 1000 Genome Project (1000 GP), we filtered out variants with minor allele frequency (MAF) < 0.01 and Hardy–Weinberg equilibrium (HWE) < 10^–6 (Auton et al., 2015). After these quality control (QC) steps, we finally obtained 7,246,475 and 7,246,543 SNPs from the two datasets.

The second category of GWAS summary statistics is from European ancestry. The dataset is from GeneATLAS website² (Canela-Xandri et al., 2018). We focus on the CLD-related phenotype that contain allele information and variant ID and that contain effect size and standard error. With the two criteria, we obtain one GWAS summary statistics: cirrhosis (n = 452,264, prevalence = 1.99%). This cirrhosis GWAS data was adjusted for sex, array batch, UK Biobank Assessment Center, age, age2 (Sudlow et al., 2015), and the top 20 genotype PCs as computed by UK Biobank. The details of these data are also provided in Supplementary Table 1. Based on European ancestry from the 1000 Genome Project, we filtered out variants with MAF < 0.01 and HWE < 10^–6 (Auton et al., 2015). After these QC steps, we finally obtained 7,636,847 SNPs from this dataset.

We treated the phase 3 of the 1000 Genome Project as the reference panel (Auton et al., 2015). Here, we collected 503 European individuals and 504 East Asian individuals with 81,271,745 SNPs. We used PLINK to calculate Pearson’s r² values of pairwise SNPs for RolyPoly with the default 1 MB window size (Chang et al., 2015). In LDSC-cts, we set the window size to 1 centiMorgan to estimate LD scores (Finucane et al., 2018).

Considering the cirrhosis and HCC data acquired from GWAS, we searched the GEO database for related scRNA-seq data and obtained one data for liver cirrhosis and two for HCC, whose raw counts data are available (Barrett et al., 2012; Ramachandran et al., 2019; Losic et al., 2020). In addition, to verify the specificity of the outcomes, we also downloaded an idiopathic Parkinson’s disease (IPD) data. The details are provided in Supplementary Table 2. Following the original study, we performed QC and clustering for each scRNA-seq data. Note that scRNA-seq data usually have the potential to have its clusters continuously subdivided, but we just controlled the cell type number of each data within 10–20 depending on the features and quality of each data. The specific processing details of each data are as follows: After demultiplexing, aligning, and estimating cell-containing partitions and associated UMIs, a cirrhosis dataset (GSE136103) consisting of CD45 + and CD45-, blood and liver, healthy and cirrhosis, and human and mice samples were downloaded (Ramachandran et al., 2019). Here, we only chose nine human cirrhotic samples, including five CD45 + and four CD45- samples, for downstream analysis.

For scRNA-seq data analysis, we first removed potential doublets, and then excluded the cells that expressed fewer than 300 genes or mitochondrial gene content >30% of the total UMI count (Ramachandran et al., 2019). We also excluded genes expressed in fewer than three cells. We followed the analysis flow in Seurat (Stuart et al., 2019): (1) used SCTransform, a new strategy to remove the influence of technical characteristics while preserving biological heterogeneity via regularized negative binomial regression, to normalize and scale scRNA-seq data (Hafemeister and Satija, 2019); (2) used default setting of IntegrateData to remove the batch effect (Butler et al., 2018); (3) performed unsupervised clustering and differential gene expression analyses on the integrated data; (4) used principal component analysis (PCA) for linear dimension reduction, and then used shared nearest neighbor (SNN) graph-based clustering, in which the graph was constructed using the top 30 principal components; and (5) used UMAP to visualize by the same number of principal components (PCs) as the associated clustering, with perplexity ranging from 30 to 300 according to the number of cells in the dataset or lineage. The details of data processing are shown in Figure 1.

In cell type definition, we referred to marker genes that are widely recognized and those from the original research. We used BuildClusterTree to assess cluster similarity by constructing the phylogenetic tree (Stuart et al., 2019). Totally, we identified 20 clusters on 23,184 cells (Supplementary Table 2 and Figure 2). Marker genes used for cell type definition are shown in Supplementary Table 3.

The first HCC dataset (GSE149614) contains 21 primary tumor, portal vein tumor thrombus (PVTT), metastatic lymph node, and non-tumor liver samples from 10 HCC patients. We downloaded the raw count data, which have been processed and aligned by Cell Ranger, and chose only 10 primary tumor samples for downstream analysis (Zheng G.X.Y. et al., 2017). After processing and clustering, we totally identified 14 cell types on 30,983 cells in this dataset (Supplementary Table 2).

Another HCC dataset (GSE112271) contains three and four tumor samples coming from different regions of two different individuals, and we included all seven samples for downstream analysis. After data processing, we totally identified 13 clusters on 32,200 cells in this dataset (Losic et al., 2020; Supplementary Table 2).

We downloaded the processed and aligned IPD dataset (GSE157783), which contains samples from six control and five idiopathic Parkinson’s disease cases. We chose only five disease samples for downstream analysis and totally identified 12 clusters on 19,002 cells following our procedure (Supplementary Table 2).

We used RolyPoly and LDSC-cts to define the specific cell types associated with cirrhosis and HCC (Calderon et al., 2017; Finucane et al., 2018). Based on polygenic model, RolyPoly treats the variance of each gene as the linear combination of each cell type and estimates the coefficients by method-of-moment. Then, RolyPoly uses block bootstrap to estimate the variance for the cell type effects, then construct t-statistics to test them (Efron and Tibshirani, 1986). By utilizing GWAS summary statistics for all SNPs near protein-coding genes, the model performed joint analysis with gene expression of a variety of cell types simultaneously, to define prioritized trait-relevant cell types (Calderon et al., 2017). We extracted the log-normalized matrix from each processed data and averaged the expression across each identified cell-type classes. We also scaled the expression data, and then took the absolute expression values, so as to form the input of RolyPoly (Calderon et al., 2017). We referred to the Ensembl database (GRCh37) and defined a 10-kb window center around the transcription start site (TSS) of a gene as its transcribed region, to construct a block annotation as recommended that could link the location of GWAS variants with related genes. Of note, we only retained genes on autosomes (Calderon et al., 2017). We used the default parameters and set 1,000 times bootstrap to obtain robust standard errors.

Based on partition heritability, LDSC-cts needs the top upregulated genes list of each cell type rather than the expression data (Finucane et al., 2018). Here, we used Wilcoxon rank sum test embedded in Seurat to find the DEGs for each cell type with all remaining clusters as control. Following Finucane et al. (2018), we extracted the top 10% upregulated genes ranked by P value from each cell type. DEGs were identified as genes expressed in at least 0.1% total cells and with log-transformed fold change above 0 in the target cluster under comparison, so as to ensure a sufficient number of genes could be obtained from each cluster. DEGs lists of each scRNA-seq data used for LDSC-cts analysis are summarized in Supplementary Tables 4,8. We referred to the Ensembl database (GRCh37) and defined the region from the TSS to the transcription end sites (TES) of a gene as its transcribed region (Yates et al., 2019). We also added 100-kb windows on either side of the transcribed region of each gene. Finally, we applied LDSC-cts by jointly modeling the annotation that corresponded to each cell type, a common annotation that included all of the genes, and the 52 annotations in the default “baseline model,” to identify CLD-specific cell types (Finucane et al., 2018).

We also made a sensitivity analysis. Specifically, we changed the resolution used in clustering to obtain a coarser cell type list for analysis. In particular, since LDSC-cts is sensitive to the gene list used for analysis, we simultaneously changed the number of genes included in LDSC-cts to the top 5% upregulated ones.

Bonferroni correction was used for multiple tests (P < 0.1/n, where n = 4 or three is the number of cell type groups, including epithelial cell, non-parenchymal cell, lymphatic immune cell, myeloid immune cell for liver tissue, or gliocyte, neuron, and vascular cell for the brain tissue, Supplementary Table 9) (Hao et al., 2020).

We used scDblFinder package (version 1.4.0), Seurat package (version 1.4.0), biomaRt package (version 2.45.6), and RolyPoly package (version 0.1.0) in R software (version 3.6.3) (R Core Team, 2020). We used PLINK (version 2.0) (Chang et al., 2015) to analyze GWAS data. We also used LDSC-cts (version 1.0.1) in python software (version 2.7.18) (Van Rossum and De Boer, 1991).

For the HCC GWAS data from BBJ, we totally retained 7,246,543 variants with HWE < 10^–6 and MAF > 0.01, as well as their annotation. For the scRNA-seq data (GSE149614), we identified 14 cell types on 30,983 cells (Supplementary Table 2 and Supplementary Figures 1,2). We further excluded cluster with less than 100 cells (63 mast cells) to avoid the interference of their unstable signal on the results. We also excluded the circulating cluster (2,510 cells), since it usually contains various immune cells from the circulation and may represent a mixed signal. Finally, we retained a total of 28,410 cells from 12 cell types. After integrative analysis, we identified B cell (β = 2.956 × 10^–4, se = 1.442 × 10^–4, P = 0.0228) as cell type associated with HCC in RolyPoly (Figure 3), whereas natural killer cell (NK), monocyte, CD4 + T cell, plasma, macrophage, hepatocyte, regulatory T cell (Treg), endotheliocyte, mesenchymal cell, CD8 + T cell, and dendritic cell (DC) showed no significance (P > 0.05). In LDSC-cts analysis, we also obtained B cell (β = 2.475 × 10^–9, se = 1.116 × 10^–9, P = 0.0133) as the significant cell type.

We used another HCC scRNA-seq data from GEO for verification. Totally, we recognized 12 cell types on 30,931 cells from the GSE112271 data with one circulating (1,192 cells) and one small cluster (77 liver sinusoids endothelial cells) excluded (Supplementary Table 2 and Supplementary Figures 3,4). We identified monocyte-derived macrophage (MDM, β = 1.665 × 10^–4, se = 6.098 × 10^–5, P = 0.0031), T cell (β = 1.732 × 10^–4, se = 7.170 × 10^–5, P = 0.0076), and natural killer cell (NK, β = 1.458 × 10^–4, se = 6.976 × 10^–5, P = 0.0191) as cell types significantly associated with HCC in RolyPoly (Figure 3), whereas the obtained NK (β = 2.331 × 10^–9, se = 1.118 × 10^–9, P = 0.0186) and B cell (β = 2.255 × 10^–9, se = 1.134 × 10^–9, P = 0.0234) as the significant cell types in LDSC-cts analysis.

We also integrated the two HCC scRNA-seq data and obtained a combined data consisting of 60,120 cells and 13 cell types for further analysis (Supplementary Figures 5,6). The RolyPoly analysis showed that B cell (β = 2.451 × 10^–4, se = 9.240 × 10^–5, P = 0.0040) was significantly associated with HCC (Figure 3), whereas the LDSC-cts identified no significant cell type.

We used scRNA-seq data from other disease to verify the specificity of our findings. To be specific, we downloaded one IPD (GSE157783) scRNA-seq data, and identified 12 cell types on 19,002 cells (Supplementary Table 2 and Supplementary Figures 7,8). After excluding clusters with too few cells (47 fibroblasts and 26 T cells), we identified no cell type significantly associated with HCC in either RolyPoly or LDSC-cts analysis (Figure 4).

We also made a sensitivity analysis by changing the resolution used in clustering and got nine, eight, and nine cell types for GSE149614, GSE112271, and their integrated data, respectively. Sensitivity analysis showed that B cell was still significantly associated with HCC in RolyPoly analysis on GSE149614 and the integrated data, as well as in LDSC-cts analysis on the integrated data. It also showed nominal significance (P < 0.1) in LDSC-cts analysis on GSE112271, and was the top cell type (P = 0.119) in the analysis on GSE149614 (Supplementary Figure 9).

For the cirrhosis GWAS data from BBJ of East Asian population, we totally retained 7,246,475 variants with their annotation. For the scRNA-seq data (GSE136103), we identified 20 cell types on 23,184 cells (Supplementary Table 2; Figure 2; Supplementary Figure 10), but further excluded circulating cluster (309 cells) and clusters with less than 100 cells (56 CLEC9A + dendritic cells and 31 mast cells). Finally, we retained a gene expression data of 17 cell types. RolyPoly showed that CD4 + T cell (β = 2.278 × 10^–4, se = 1.149 × 10^–4, P = 0.0259) was significantly associated with cirrhosis, whereas LDSC-cts identified no significant cell type (Figure 5).

We also used a cirrhosis GWAS summary data of European population from GeneATLAS website to verify the stability of our outcomes, in which a total of 7,636,847 variants was retained after QC. We identified natural killer T cell (NKT, β = 6.535 × 10^–10, se = 2.423 × 10^–10–1.110 × 10^–9, P = 0.0038) and hepatocyte (β = 2.891 × 10^–10, se = 1.364 × 10^–10, P = 0.0149) as cell types significantly associated with cirrhosis in RolyPoly, while we obtained no significant cell type in the LDSC-cts analysis (Figure 5).