Cohorts
Recruitment
Blood samples were collected from 2,339 consenting individuals from two Australian biomedical research cohorts: the Tasmanian Ophthalmic Biobank (TOB) (ethics protocol 2020/ETH02479) and the BioHEART study (ethics protocol 2019/ETH08376) cohort.
The TOB cohort consisted of 983 individuals (562 female, 421 male) recruited between January and February 2018 by the Royal Hobart Hospital and Menzies Institute of Medical Research in Tasmania, Australia. The selection criteria for inclusion in the cohort were as follows: (1) age 18 years and over, (2) no signs of ocular disease and (3) self-reported British, Scottish or Irish ancestry.
The BioHEART63 cohort included 1,356 individuals (561 female and 795 male) recruited between 29 March 2016 and 2 May 2023, presenting to study sites for clinically indicated cardiovascular imaging or acutely with ST elevation myocardial infarction. The selection criteria for inclusion in the cohort were as follows: (1) age 18 years and over and (2) presenting with cardiac disorders (or for the investigation of suspected cardiac disorders) including: acute myocardial infarction, suspected coronary artery disease, patients undergoing cardiothoracic or vascular surgery, or plastic surgery involving resection of vascularized tissue.
After donor, variant and cell quality control, including selection to individuals of European ancestry only, 1,925 individuals remained: 950 from TOB and 975 from BioHEART (Supplementary Fig. 1).
Ethics oversight
This study complied with all relevant ethical regulations and was overseen by the Human Research Ethics Committees at the Royal Children’s Hospital (Melbourne, Australia), Northern Sydney Local Health District (Sydney, Australia), St. Vincent’s Hospital (Sydney, Australia) and University of Tasmania (Hobart, Australia). Informed consent was obtained from all participants.
Blood collection and DNA extraction
For the TOB16 cohort, blood samples were collected from each participant and stored in 1× 8-ml heparin cell preparation tubes and 1× 9-ml EDTA tubes. For the BioHEART cohort, aliquots were stored at −80 °C.
DNA extraction from BioHEART buffy coat samples was performed by the Australian Genome Research Facility using QIAamp96 DNA QIAcube HT kit. DNA extraction from TOB samples was performed at the Menzies Institute of Medical Research (University of Tasmania) using the QIAmp Blood Mini Kit. For a subset of TOB samples that failed initial quality control checks during sequencing, DNA was extracted from PBMCs a second time using the Qiagen DNeasy Blood and Tissue Kit and eluted in an AE buffer by the Garvan Molecular Genetics laboratory service. DNA concentration was measured for all samples using fluorometric quantitation with the Quantifluor from Promega before sequencing.
WGS data
WGS data were generated for all individuals from the TOB cohort (n = 983) and BioHEART (n = 1,356).
Sequencing library preparation and WGS
Sequencing library preparation was performed using the KAPA PCR-free v.2.1 kit (Roche) for the TOB samples, and the DNA PCR-free kit (Illumina) for the BioHEART samples. Short-read WGS was performed using Illumina 2 × 150 bp chemistry on the NovaSeq 6000 instrument to achieve mean 30× coverage. Library preparation and sequencing were performed by the Garvan Sequencing Platform.
Alignment and variant calling
Alignment and variant calling were performed in line with GATK-DRAGEN Best Practices. Reads were aligned from FASTQ (BioHEART) or CRAMs (TOB) against the GRCh38 reference genome (masked version, as available at https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle) using a fork of DRAGMAP v.1.3.1 (https://github.com/populationgenomics/DRAGMAP). When aligning from CRAMs, bazam (https://github.com/ssadedin/bazam) was used to first convert the reads back to FASTQ format. Calling of SNPs and indels was then performed using GATK4 HaplotypeCaller in DRAGEN mode. The resulting gVCF files were combined using Hail’s VDS Combiner and stored in the variant dataset (VDS) format.
Sample quality control
Hard sample filtering was applied to remove individuals with: high rates of contamination (FREEMIX ≥ 5%); an increased percentage of chimeras (PCT_CHIMERAS ≥ 6.3%); median insert size less than 250 bp and sex aneuploidy or discordance between reported and observed sex. The dataset was also inspected for outlier samples based on the number of singletons, heterozygous /homozygous variant ratio, transition/transversion ratio and insertion/deletion ratio, taking care not to exclude ancestry outliers on the basis of these metrics.
A set of high-quality sites, defined as an LD-pruned subset of autosomal, biallelic, allele frequency greater than 0.01, call rate greater than 99% and InbreedingCoeff greater than −0.25, were extracted from 1,000 Genomes and the Human Genome Diversity Project. Samples with call rate less than 0.85 across these sites were also excluded.
Following hard sample filtering, related individuals were identified using Hail’s implementation of PC-Relate and confirmed using somalier relate v.0.2.15 (https://github.com/brentp/somalier). For family sets consisting of second-degree or closer relatives, as defined by a kinship coefficient greater than 0.088, the sample with the highest quality scRNA-seq data was retained, as determined by the number of cells and cell-type composition. When equal quality single-cell data were available, the sample with the highest mean genome sequencing depth on chromosome 20 was retained.
Ancestry inference was then performed on the high-quality, unrelated set of samples, by combining this set with background samples from 1,000 Genomes and the Human Genome Diversity Project, then computing principal components using gnomAD’s run_pca_with_relateds method and the same set of high-quality sites defined for call rate estimation above. Ancestry labels were assigned on the basis of genetic similarity, as identified through principal component analysis, to those of known ancestry in the background datasets. For association testing, only samples assigned to the European cluster were retained (Supplementary Fig. 39).
Variant quality control
Allele-specific Variant Quality Score Recalibration or AS-VQSR (GATK v.4.2.6.1) was run using all quality control-pass, unrelated samples (regardless of ancestry), using the default training resources and priors. For SNPs this is hapmap (prior, 15), omni (prior, 12), 1000 G (prior, 10) and dbsnp (prior, 7). For indels this is mills (prior, 12), axiomPoly (prior, 10) and dbsnp (prior, 2). The following features were included for all SNPs and indels: AS_ReadPosRankSum, AS_MQRankSum, AS_QD, AS_FS, AS_SOR. AS_MQ was used for SNPs only. All sites failing on the site-level filter calculated by AS-VQSR were removed.
Within the European-only association callset, we then removed any variants that had no high-quality genotype calls (defined as genoype quality (GQ) ≥ 20, depth of coverage (DP) ≥ 10 and allele balance greater than 0.2 for heterozygotes) or with an excess of heterozygotes at that site compared with Hardy–Weinberg expectations (InbreedingCoeff < −0.3). A final filter to remove sites with missingness greater than 15% was then applied before association testing.
scRNA-seq data
scRNA-seq data were collected from PBMCs from samples from individuals from the TOB (OneK1K16; n = 1,036) and BioHEART63 (n = 1,012) cohorts. Supplementary Fig. 40 illustrates the experimental workflow and key steps.
Sample pooling and cell count estimation
TOB samples were previously pooled during preparation of the OneK1K project, with pools made up of 1,214 individuals16. To capture more cells per individual (aiming for up to 5,000 cells per individual), two pools (n = 24–28) at a time were further pooled (weighted by cell viability) and loaded into three sequencing lanes (aiming for 40,000 cells each, so 120,000 in total). In addition, ‘maxi pools’ of four combined pools (96–112 individuals) were sequenced in two separate lanes to compensate for possible low cell numbers.
By contrast, BioHEART samples were pooled (weighted by cell viability) for this project; each pool comprised 8 individuals and was sequenced in a separate sequencing lane. We aim to include 40,000 cells (roughly 5,000 cells per individual) per pool. In some cases, pools with low cell viability were sequenced again to rescue cell numbers and get closer to the target.
Capture and sequencing
Cells were captured using the single-cell 3′ gene expression high-throughput (3′ GEX HT) kit v.3.1 from 10X on a Chromium X machine. For the TOB cohort, cells were recaptured starting with the original biomaterial (OneK1K individuals) that we had biobanked and using the new 10X chemistry with the new high-throughput kit. The long-term cryopreservation of these samples led to a smaller proportion of monocytes in the TOB cohort compared with BioHEART (Supplementary Note 1). Library preparation and short-read sequencing were performed using Illumina’s NovaSeq 6000 machine.
Alignment
scRNA-seq reads were aligned to a genome reference (GRCh38, Gencode release 44, Ensembl 110) using STAR64, as implemented within 10X’s CellRanger software (v.7.2.0) that also estimates the number of cells captured (Supplementary Fig. 41). CellBender65 (v.0.3.0) estimated numbers of cells differed substantially from those estimated by CellRanger, leading us to discard these results (Supplementary Note 2).
Demultiplexing and doublet detection
To assign cells to their donor of origin, as well as to remove doublets (droplets containing more than one cell), we used Demuxafy66 (v.2.1.0), a tool to runs a number of demultiplexing and double detection software and outputs a consensus call for each individual cell. Using Demuxafy, we ran vireo67 to demultiplex cells to their donor of origin within each sequencing library, using genotypes extracted from either SNP array data (available as part of the OneK1K data; for TOB) or using common (MAF > 5%) biallelic SNPs from WGS (for BioHEART). A handful of individuals did not have available genotype data, but vireo allows us to include the expected number of individuals in each pool with or without genotypes, which makes the donor assignment and doublet detection more accurate. For doublet detection, we used Demuxafy to run scds68 and ScDblFinder69 and used majority voting between all three tools (vireo, scds and scDblFinder) to determine the final set of doublets for each sequencing library (Supplementary Fig. 42).
Cell typing
To assign individual cells to specific immune cell types we follow the recommendations of the single-cell eQTLGen consortium28. We adopt a combined approach that uses scPred70, hierarchical progressive learning, scHPL71 and Azimuth as the reference dataset72 as detailed in https://powellgenomicslab.github.io/WG2-pipeline-classification-docs/index.html. To run this step, we converted our single-cell Scanpy object into Seurat73 objects. We also used CellTypist74,75 (v.1.6.2, model used ‘Immune_All_High’) to confirm cell typing of individual cells, with good concordance (Supplementary Fig. 43). As the former approach provided more granularity (28 cell types versus 18 identified by CellTypist), we proceeded with the 28 cell types for the remainder of the article.
Quality control
Next, we performed quality control with Scanpy (v.1.8.2)76. For each library, we removed cells with more than 20% mitochondrial reads, fewer than 800 reads detected and fewer than 1,000 or more than 10,000 genes detected (Supplementary Fig. 44). We also removed individuals with abnormal cell-type composition and individuals with fewer than 100 cells.
Overall, after cell and sample quality control (based on scRNA-seq and WGS), a total of 5,438,679 cells from 1,925 individuals were left for consequent analyses (Supplementary Fig. 1).
Association testing
SAIGE-QTL
To test for the effects of sets of both common and rare-variant associations with single-cell gene expression, we used the recently proposed SAIGE-QTL method27. In brief, SAIGE-QTL implements a Poisson mixed model, thus effectively modelling single-cell counts and repeat cells for the same individual. Furthermore, several computational and statistical speed-ups are implemented to improve speed, scalability and storage usage. Finally, implementation of both single-variant tests (in which effects of the genotype at a single locus are assessed), as well as set-based tests (in which signals from several variants within a region are assessed together in their combined effect), is available, as briefly described below.
We mapped eQTLs separately in each of the 28 cell types. We confirmed that our tests were calibrated by shuffling the individual IDs and reperforming the test for our most abundant cell type, CD4TCM cells, the results of which show no evidence of false positives for both common and rare-variant tests (Fig. 2a).
Fitting the null
For every gene and cell type, we fit a single null model using SAIGE-QTL, the same for both subsequent tests (common and rare-variant tests alike). Briefly, we fit a Poisson model using the raw counts from the single-cell expression data as the outcome variable (phenotype), and the following as covariates: sex and age of the individuals, seven genotype principal components, five cell-type-specific expression principal components, cohort status (TOB or BioHEART), percentage of counts from mitochondrial genes and total counts. We provide all extra parameters and flags as part of the available code in GitHub. We only test genes that are expressed in at least 1% of cells in each cell type. For a subset of gene–cell-type combinations the model failed to converge (0.27%) and these were excluded from the analysis (Supplementary Note 3).
Common-variant eQTLs
We test for the effects of common genetic variants on single-cell expression using the single-variant test as implemented in SAIGE-QTL. We consider all SNPs and indels with MAF ≥ 1% and missingness less than 15% in our sample. We test for cis (variants within a ±100-kb window around the gene body) effects only. We apply a two-step multiple testing correction as recommended in the SAIGE-QTL paper27, using the Aggregated Cauchy Association Test77 at gene level and Storey’s q value method (qvalue package in R) across genes. At FDR < 5%, we identify a total of 154,932 eQTLs in at least 1 of the 28 cell types, from 17,674 unique genes (83% of all 21,404 genes tested).
The number of eQTLs identified varied substantially across cell types, ranging from 11,491 eQTLs in CD4TCM cells (85% of all genes tested in this cell type) to 449 eQTLs in innate lymphoid cells (3%) as shown in Fig. 2b. This was probably explained by differences in the statistical power across cell types, largely due to the number of cells available for each cell type (correlation between log10 number of cells and number of eGenes detected, 0.98; Extended Data Fig. 1).
The same analysis was run again for all cell types (common variants only) using a ±1-Mb cis window around the gene instead. These results are not used for any of the analyses described in the main paper (except Supplementary Fig. 45), but we made the summary statistics available, see ‘Data availability’ for details.
Rare-variant eQTLs
We test for the effects of rare variants on single-cell expression using the set-based test as implemented in SAIGE-QTL. We consider all SNPs and indels with MAF < 1% and missingness less than 15%. We test for the combined effects of all variants within a cis region around the gene (±100 kb around the gene body). As described in the SAIGE-QTL publication, we consider a burden test and the SKAT test and use the Cauchy method to combine P values at the gene level. Similar to the approach used for the common-variant results above, we use the q value method and report significant results at FDR < 5%.
Overall, we observe 34,376 rare-variant signals across all cell types, for 10,005 unique genes after FDR correction using Storey’s q value on the Cauchy combined P values from the two tests (the burden test and SKAT). The percentage of genes with a rare-variant combined signal out of all gene tested ranged from 0.09% (innate lymphoid cells) to 39% (CD4TCM cells; Fig. 2b). Once again, the tests were well calibrated when shuffling individual IDs (Fig. 2a) and the number of cells for each cell type was highly correlated with the number of eGenes identified (Extended Data Fig. 1). We discarded results from the Aggregated Cauchy Association Test-V test because we found them to be not calibrated in this setting (Supplementary Note 4).
An equivalent analysis was performed after restricting the association to only (cis rare) variants that were functionally annotated (details below).
Single-cell ATAC-seq processing
Open-chromatin status annotation of case examples
Single-cell ATAC-seq data were obtained from a subset of the TOB cohort that includes 30,829 nuclei from 8 individuals26. Sequencing reads were aligned to the Human Genome Reference GRCh38-1.2.0 (gencode v.44) using Cell Ranger ATAC (v.2.1.0) and the chromatin peaks were called using the built-in MACS2 (ref. 78). Cell-type labels were annotated based on a CITE-seq reference72 with a combined Multiome bridging dataset26, consistent with the 28 PBMC cell types annotated for scRNA-seq data. After quality control, 17,099 nuclei were retained to annotate the chromatin peaks in each cell type for the eQTL signals.
Open-chromatin status annotation of all rare variants for association set tests
To annotate all rare variants before association testing (in the run using only variants in likely-to-be-functional regions), we used the TenK10K scATAC-seq data to annotate variants found to be in open-chromatin regions for the relevant cell type. We first obtained the open-chromatin peak list for each of 26 cell types (2 removed due to super rarity; innate lymphoid cells; n = 33 and CD8Proliferating; n = 4) from ref. 26. For each cell type, we re-called peaks from the raw fragment files using the MACS3 tool (v.3.0.3) integrated into SnapATAC2 (v.2.8.0). The number of peaks identified (q < 0.05) per cell type ranged from 43,014 in CD4Proliferating to 269,371 in CD16Mono, with a median of 181,702. Then we set up a searching function to check whether the rare variant is located within any peaks in the corresponding cell type.
Functional annotation of variants
Publicly available resources
We used publicly accessible resources to annotate our variants, including ENCODE79,80,81 (v.3) ciscandidate regulatory elements (cCREs), University of California, Santa Cruz (UCSC) genome annotation database82, Activity-by-Contact (ABC)83,84 scores, FANTOM585 enhancer annotations, functional assay and conservation or constraint results from a deep literature search, chromHMM86 predictions and enhancers from the VISTA Enhancer Database87 with positive signal. Any annotations not from genome build GRCh38 were converted using the UCSC LiftOver tool. Any LiftOver segments mapping to non-canonical chromosomes and/or fragmented less than 100 bp in GRCh38 were dropped. All overlapping intervals within categories were merged using bedtools merge88. Detailed annotation information and references are summarized in Supplementary Table 3.
PromoterAI and SpliceAI to predict promoter and splicing activity of variants
PromoterAI is a convolutional neural network that uses roughly 20 kb of sequence context around a promoter variant to predict its effect on gene expression. It is pretrained on histone modifications, DNA accessibility, transcription factor occupancy and strand-specific CAGE value around TSS at base-pair resolution and fine-tuned with a curated set of rare non-coding promoter variants associated with unusually high or low expression in GTEx v.8. PromoterAI scores range from −1 to 1, with negative values indicating a predicted under-expression effect and positive values indicating a predicted overexpression effect. More details on PromoterAI can be found in the publication32 and common variant lead eQTLs are annotated with their PromoterAI score (when applicable) in Supplementary Table 4.
SpliceAI35 was run on all coding tested variants (common and rare) and the scores are available in Zenodo89. AlphaGenome38 and Pangolin90 were also used to analyse the possible mis-splicing outcomes of two variants.
Annotations for rare-variant set-test analyses
We used an aggregated ‘functional and open’ flag to label any variant sitting in a region considered to be in any of a promoter, enhancer or actively transcribed region, in a transcription factor binding site (Supplementary Table 3), predicted to be splicing-altering or found to reside in an open-chromatin region in the relevant cell type based on our own scATAC-seq (details above).
Fine mapping of common-variant eQTLs
We ran SuSiE29 (R package ‘susieR’, v.0.14.2) for fine mapping. Specifically, we used the ‘susie_rss’ function using L = 10, Z scores from the SAIGE-QTL summary statistics, effective N (number of individuals tested, which varied across cell types, Fig. 2b) and correlation matrices we calculated for all-cis variants for each gene (noting that SAIGE-QTL using a Poisson model violates SuSiE’s linear assumption, Supplementary Note 5). We considered each credible set as an ‘independent signal’, and we defined as a ‘primary’ signal the credible set containing the variant with the smallest P value, which is the same as the main eQTL for each eGene reported throughout the rest of the article. For genes with more than one independent signal, all other credible sets were defined as ‘non-primary’. The credible set size in Fig. 2c refers to the number of variants contained in each credible set. For the comparison between primary and non-primary signals (Supplementary Fig. 6) all variants with PIP > 0.9 in either the primary credible set or any of the non-primary credible sets for every gene were considered in terms of both their distance from the TSS and their functional annotations.
All cell types pseudobulk analysis
We aggregated expression values by taking the mean value at the level of gene and individual for each cell type. We then used each cell-type pseudobulk as a biological replicate and ran SAIGE-QTL (Poisson) using all the same parameters as for the main results presented in the article. In the comparisons to the single-cell cell-type-specific analyses (Extended Data Figs. 4 and 8b), only genes tested in the pseudobulk analysis as well were retained (n = 21,192).
Cell-type specificity analysis
To quantify cell-type specificity, for each of the 154,932 identified eQTLs, we first performed pairwise comparisons of association estimates in the primary cell type and in each of the other 27 cell types. We classified these pairwise comparisons into five scenarios based on eQTL association, eGene expression and direction of allelic effect on gene expression in the two cell types: (1) eGene is expressed only in the primary cell type, (2) eGene is expressed in both cell types, but eQTL identified only in the primary cell type, (3) independent eQTL variants (LD r2 < 0.5) for the same eGene are identified in each cell type, (4) same eQTL variants or different variants in LD r2 ≥ 0.5) for the same eGene are identified in each cell type with concordant direction of effect and (5) same eQTL variants or different variants in LD r2 ≥ 0.5) for the same eGene are identified in each cell type with discordant direction of effect. These scenarios are summarized in Supplementary Table 8.
Further, we used a pairwise meta-analysis approach for all eQTLs that are initially classified as in scenario 2 to take into account potential type 2 error due to insufficient statistical power for association analysis in rarer cell types. This was done for each significant eQTL (FDR < 5%), by comparing its effect size in the primary cell type to its effect sizes in each of the other cell types. Meta-analysis was performed using a random-effects model with the metafor package in R91 (v.4.0.0). We compared the resulting meta-analysis P value to the primary Pn value, denoting an empirically chosen threshold of 105-fold reduction (Supplementary Fig. 38) as evidence of subdiscovery eQTL signal in cell type 2, in which case the eQTL in question would be reclassified as scenario 4 or 5 depending on concordance of allelic effect direction.
Next, we calculated an aggregate score for each of the 17,674 common-variant eGenes by taking the maximum value of the scenarios above across all pairs of cell types for each eGene (results presented in Fig. 4a, bar plot on the right). Supplementary Table 1 describes interpretation of the scenario or aggregated values in each of the comparisons.
Given that correlation of effect estimates are proportional to LD r2 between variants, we noted that variants in lower LD are more likely to show a discordant allelic effect direction and therefore may result in overestimation of aggregated count for scenario 5. We repeated this approach using the gene-level rare-variant association results. For each gene in a given cell type identified at FDR < 5% based on the Cauchy-aggregation P value (from burden test and the SKAT test), we compared the estimates for the same gene against those derived from each of the other 27 cell types. Similar to the common-variant results above, each pairwise comparison was initially classified into one of the following scenarios: (1) a gene was expressed in the primary cell type only, (2) a gene was expressed in both cell types, but the rare-variant association signal was identified only in cell type 1 and (3) the rare-variant association signal was identified in both cell types. To account for potential misclassification of scenario 2 genes in which a rare-variant association signal was absent due to lack of statistical power in the association test in the secondary cell type, we meta-analysed effect estimates from burden test results of both cell types using a random-effects model. Genes with a meta-analysis P value 105-fold smaller than the original meta-analysis were reclassified as probably shared eGenes. Unlike the common-variant results, we did not perform a comparison of direction of effect sizes given that the discovery was done using unsigned P values from Cauchy aggregation of burden and SKAT tests.
Statistical colocalization
We use the COLOC92 package to perform statistical colocalization between our eQTLs and a number of traits and diseases using GWAS summary stats from publicly available sources.
Specifically, we downloaded GWAS summary statistics for the following 16 disease traits: Alzheimer’s disease93, asthma94, breast cancer95, COVID-1996, Crohn’s disease97, IgA nephropathy98, IBD97, lung cancer99, malignant lymphoma36, multiple sclerosis100, Parkinson’s disease101, prostate cancer102, rheumatoid arthritis103, SLE104, type 1 diabetes105 and ulcerative colitis97.
In addition, we downloaded UKBB GWAS summary statistics for 44 blood cell count and serum traits40 (SNPs and indels only) (from https://gymreklab.com/science/2023/09/08/Margoliash-et-al-paper.html).
Colocalization was performed using the coloc.abf function from COLOC (v.5.2.3)92. Colocalization was performed using SNP GWAS summary statistics and common-variant eQTLs in our dataset with eQTL type set to ‘quant’ and providing N and MAF rather than sdY as single-cell expression was not normalized and modelled as a Poisson distribution by SAIGE-QTL.
Cell-state eQTL analysis
Cell-state definition
We computed cell function scores as cell states using scDeepID, a multi-task transformer, in which a biologically informed masked-encoder was used to generate cell function embeddings along the supplied database, followed by combining a
To define the cell states in the main NK cells of the TenK10K cohort, we first trained our scDeepID model on a publicly available dataset that contains different NK subtypes52, then performed label transfer by the trained model on our data. Specifically, we used Scanpy to automatically select roughly 600 highly variable genes (min_mean=0.0125, max_mean=3, min_disp=0.5) in the reference dataset, then only kept the shared genes between the reference and TenK10K cohort to train the model (that is, 560 genes). In particular, we further defined NKintermediate, NKadaptive and NKmature cells in the NK cell type. Then we trained a separate scDeepID model with the three fine-grained cell types as well as NKCD56bright and NKProliferating cells (with 2,315 highly variable genes, selected by Scanpy with min_mean=0.0125, max_mean=3, min_disp=0.5). We used the human gene ontology biological process database as a guide to train the scDeepID models above. We then extracted the attention-based latent space and normalized each cell function score with a zero mean and unit variance.
Cell-state eQTL mapping
We used the cell states defined above to test for genotype × cell-state interactions. We focused on all eQTLs mapped in any of the NK cells subpopulations (NK, NKProliferating and NKCD56bright cells n = 16,759 unique variant–gene pairs) and used lme4 (ref. 107) to test the significance of the interaction term cell state × genotype for each eQTL and each of the five cell states, while correcting for age, sex, cohort (BioHEART or TOB), total counts as fixed effect covariates and including an individual term as a random effect.
The alternative model tested was:
$${\bf{sc.\; exprs}}={\rm{covs}}+{\bf{cell\; state}}+{\bf{genotype}}+{\bf{cell\; state}}\times {\bf{genotype}}+{\bf{individual}}+{\varepsilon }$$
where sc. exprs is the single-cell expression profile of the target eGene, ‘cell state’ is one of the cell-state functions estimated using scDeepID, genotype is the genotype vector at the variant tested, encoded as the number of minor alleles for each individual (0, 1, 2; expanded such that each cells from a given individual has the genotype of its individual of origin), ‘individual ’ is a random-effect term to account for the fact that each cell is not an independent observation and \(\varepsilon \) is a Gaussian error.
csaQTL analysis
csaQTL mapping
We performed csaQTL associations using GeNA56 v.1.0.0. We generated single cell objects for each main cell type group and applied the following processing to each object. Following a standard Scanpy pipeline, we normalized read counts per 10,000, log-transformed normalized counts and selected roughly 3,000 highly variable genes per main cell type (min_mean=0.0125, max_mean=3, min_disp=0.5). The scanpy.pp.regress_out function was used to obtain residualized counts, removing technical variation associated with total counts (library size) and percentage of mitochondrial reads. Expression for each gene was scaled to unit variance across cells. We then created MultiAnnData objects containing the processed counts and donor covariate information. We applied principal component analysis with Scanpy using default parameters. Harmony108 (harmonypy v.0.0.10) was run with default parameters, using sequencing library pools as batch assignment to correct principal components for technical batch effects between sequencing libraries. The harmonized principal components were used to construct a nearest neighbour graph and uniform manifold approximation and projection (UMAP) using Scanpy default parameters.
Next, we generated a neighbourhood abundance matrix (NAM) and its principal components with CNA109 (v.0.1.6) using default parameters. We then applied GeNA, controlling for age, sex, cohort and the top seven genotype principal components. We tested associations for all common variants (MAF ≥ 0.01 and percentage missing genotypes less than 15%). To calibrate the model, we also ran GeNA on permuted genotypes generated by random label swapping (one permutation performed per cell type). On the basis of this permutation testing, we observed poor calibration for variants below 0.05 MAF and removed them from the analysis. We considered variants passing the genome-wide significance threshold (P < 5 × 10−8) as statistically significant. Lead SNPs were defined using the clumping method as described in the GeNA publication56. In brief, at each locus the most significant SNP was considered the lead SNP, and all SNPs with LD > 0.8 to the lead SNP or within a 1-Mb window were considered LD-SNPs. We considered the number of unique lead SNPs for csaQTL loci as the total independent csaQTLs discovered. To visualize cell counts adjusted for covariates, we calculated the fraction of cells within each main subtype belonging to each minor subtype. Using these cell fractions, we fit a linear model with the lm function in R, including age, sex, cohort and the first seven genotype principal components as covariates, and visualized the residuals from the fitted model.
Characterization of csaQTL-associated cell states
To characterize the cell-state phenotypes detected by GeNA, we used methods previously described in ref. 56 briefly, we calculated the Pearson correlation between GeNA neighbourhood-level phenotypes and expression of highly variable genes. Highly variable genes were then ranked by Pearson’s r for each csaQTL and used as input to GSEA. The fgsea (v.1.32.2) R package110 was used to run GSEA for all Molecular Signatures Database hallmark gene sets. Gene sets with Benjamini–Hochberg-adjusted P < 0.05 were considered statistically significant. To visualize the expression of enriched gene sets, module scores were calculated with the Scanpy score_genes function and overlayed on to UMAP plots.
Colocalization analysis between csaQTLs and blood cell count GWAS
To obtain GeNA summary statistics compatible with colocalization analysis, we calculated fixed-phenotype summary statistics for each significant lead csaQTL using methods described in ref. 56. Colocalization was run with the coloc.abf function from the COLOC R package92, using default parameters. We provided COLOC with these fixed-phenotype summary statistics, as well as GWAS from two blood cell count and serum trait GWAS studies. We used the UKBB GWAS summary statistics for 44 blood cell count and serum traits described above, supplemented with an extra 26 blood trait GWAS performed in ref. 58. Colocalizations with PPH4 > 0.8 were considered significant.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
