Single-cell multiome of the human retina and deep learning nominate causal variants in complex eye diseases

Abstract

Genome-wide association studies (GWASs) of eye disorders have identified hundreds of genetic variants associated with ocular disease. However, the vast majority of these variants are noncoding, making it challenging to interpret their function. Here we present a joint single-cell atlas of gene expression and chromatin accessibility of the adult human retina with more than 50,000 cells, which we used to analyze single-nucleotide polymorphisms (SNPs) implicated by GWASs of age-related macular degeneration, glaucoma, diabetic retinopathy, myopia, and type 2 macular telangiectasia. We integrate this atlas with a HiChIP enhancer connectome, expression quantitative trait loci (eQTL) data, and base-resolution deep learning models to predict noncoding SNPs with causal roles in eye disease, assess SNP impact on transcription factor binding, and define their known and novel target genes. Our efforts nominate pathogenic SNP-target gene interactions for multiple vision disorders and provide a potentially powerful resource for interpreting noncoding variation in the eye.

Graphical abstract

Figure 1

Transcriptional profiles from joint scRNA- and ATAC-seq identify major cell types of the human retina (A) Schematic of the human retina, depicting the cell types analyzed in this study. (B) Uniform manifold approximation and projection (UMAP) plot of the 51,645 human retinal cells detected by scRNA-seq after quality control filtering and removal of putative doublets. Eight postmortem retinas from four donors were profiled. A total of 22 clusters were resolved and assigned to 13 cell types. (C) Frequency of different cell types in the human retina as determined by scRNA-seq. Numbers above each bar denote absolute counts out of 51,645. (D) Dot plot visualizing the normalized RNA expression of selected marker genes by cell type. The color and size of each dot correspond to the average expression level and fraction of expressing cells, respectively.

Figure 2

Chromatin accessibility profiles from joint scRNA- and ATAC-seq of the human retina reveal cell-type-specific epigenetic landscapes (A) Number of chromatin accessibility peaks for each cell type as determined by scATAC-seq. Peaks were required to be present in a least two pseudo-bulk ATAC replicates (n = 2 for astrocyte and microglia, n = 5 for all other cell types). (B) Overlap of scATAC peaks with peaks from published human retina bulk ATAC-seq data. Overlapping was defined as peaks with any overlapping bases. (C) Heatmap of scATAC marker peaks enriched in each cell type. Each column represents a marker peak. (D) Sequencing tracks of chromatin accessibility near selected marker genes by cell type. Each track represents the aggregate scATAC signal of all cells from the given cell type normalized by the total number of reads in TSS regions. Genes in the sense direction (TSS on the left) are shown in red, and genes in the antisense direction (TSS on the right) are shown in blue. Coordinates for each region are as follows: PDE6A (chr5:149924792–149964793), GRIK1 (chr21:29905031−29955033), RLBP1 (chr15:89201750−89241751), GRM6 (chr5:178975297−179015298), NIF3L1 (chr2:200874325−200914327), ARR3 (chrX:70248304−70288305), GAD2 (chr10:26186306−26246307), ONECUT1 (chr15:52781076−52821078), SLC6A9 (chr1:44005465−44035467), NEFL (chr8:24937109−24977110), CALB2 (chr16:71323711−71368713), PAX2 (chr10:100715602−100755603), and HLA-DRA (chr6:32419841−32459842).

Figure 3

Motif analysis of accessible DNA regions in the human retina predicts cell-type-specific TFs (A) Heatmap of selected TF binding motifs enriched in each cell type. Darker colors indicate more significant enrichment. (B) Footprinting analysis of selected TFs across cell types. Footprints were corrected for Tn5 insertion bias by subtracting the Tn5 insertion signal from the footprinting signal.

Figure 4

Single-cell multiomics pinpoints the cellular targets of noncoding variants in eye diseases (A) Overview of SNP selection for interrogating ocular disease GWASs. Index SNPs obtained from GWASs of each disease were subjected to LD expansion, and the resulting noncoding SNPs intersected with scATAC peaks. (B) Percentage of LD expanded noncoding SNPs from each disease that overlapped with chromatin accessibility peaks for each cell type. (C) Number of scATAC peaks co-accessible with each promoter peak. Co-accessible was defined as scATAC peaks whose accessibility showed a correlation score greater than 0.3. (D) Number of predicted target genes for each scATAC peak. Predicted target genes were defined as genes whose RNA expression showed a correlation score >0.3 relative to the accessibility of the tested scATAC peak. (E) Sequencing tracks of chromatin accessibility near rs4821699 (chr22:37719685) and rs17421627 (chr5:88551768). Genes in the sense and antisense directions are shown in red and blue, respectively. The location of each SNP is depicted by a vertical gray line. Gray arcs indicate predicted target genes for the scATAC peak containing the SNP of interest.

Figure 5

Integration of the single-cell multiome with HiChIP and eQTL data prioritizes functional noncoding polymorphisms in the human retina (A) Overlap of H3K27ac HiChIP loop anchors (n = 2 biological replicates) with scATAC peaks. (B) Percentage of SNPs in scATAC peaks for each disease with available retina eQTL data. (C) Sequencing tracks of chromatin accessibility near rs9966620 (chr18:24100771), rs2730260 (chr7:159054238), and rs66475830 (chr6:116087639). Genes in the sense and antisense directions are shown in red and blue, respectively. The location of each SNP is depicted by a vertical gray line. Gray arcs indicate predicted target genes for the scATAC peak containing the SNP of interest. The black arc overlapping with rs9966620 indicates a H3K27ac HiChIP loop with the region encompassed by the opposite anchor, highlighted in purple. (D) Significance of SNP-gene associations for rs2730260 or rs66475830 and their nearby genes, as determined by retina eQTL analysis. Adjusted p values for each gene were calculated by multiplying the nominal p value listed in the EyeGEx database by the number of SNP-gene pairs tested for that SNP.

Figure 6

Integration of the single-cell multiome with base-resolution deep learning nominates functional mechanisms for disease-associated SNPs (A) Schematic of the CNN-based deep learning pipeline. (B) Percentage of noncoding index SNPs (n = 1,284), LD expanded SNPs (n = 7,034), LD expanded SNPs in scATAC peaks (n = 1,152), randomly selected GC-matched SNPs (n = 9,984), and randomly selected SNPs in scATAC peaks (n = 1,160) that were categorized as high-effect. (C) Top: predicted per-base accessibility for rs1532278 (chr8:27608798) and rs1874459 (chr16:65041801) in Müller glia and rod bipolar cells, respectively, as determined by deep learning models. A 100-bp window depicts the importance of each base to predicted accessibility at the SNP, and a 1,000-bp window depicts predicted per-base counts for the reference (blue) and alternate (orange) alleles. SNP bases are highlighted in purple. For rs1874459, similar changes in accessibility were predicted for OFF-cone bipolar, ON-cone bipolar, gly-amacrine, and AII-amacrine cells. Bottom: sequencing tracks of chromatin accessibility near rs1532278 and rs1874459. Genes in the sense and antisense directions are shown in red and blue, respectively. The location of each SNP is depicted by a vertical gray line. Gray arcs indicate predicted target genes for the scATAC peak containing the SNP of interest. (D and F) Significance of SNP-gene associations for rs1532278 (D) or rs1874459 (F) and their nearby genes, as determined by retina eQTL analysis. Adjusted p values for each gene were calculated by multiplying the nominal p value listed in the EyeGEx database by the number of SNP-gene pairs tested for that SNP. (E) Dot plot visualizing the normalized RNA expression of 40 different homeodomain TFs in Müller glia. The selected TFs correspond to the 40 homeodomain factors whose binding motifs were most significantly enriched in Müller glia, as determined by motif analysis (Data S4). (G) Dot plot visualizing the normalized RNA expression of neuroD and neurogenin family members by cell type.