Prioritizing Autism Risk Genes using Personalized Graphical Models Estimated from Single Cell RNA-seq Data

Abstract

Hundreds of autism risk genes have been reported recently, mainly based on genetic studies where these risk genes have more de novo mutations in autism subjects than healthy controls. However, as a complex disease, autism is likely associated with more risk genes and many of them may not be identifiable through de novo mutations. We hypothesize that more autism risk genes can be identified through their connections with known autism risk genes in personalized gene-gene interaction graphs. We estimate such personalized graphs using single cell RNA sequencing (scRNA-seq) while appropriately modeling the cell dependence and possible zero-inflation in the scRNA-seq data. The sample size, which is the number of cells per individual, ranges from 891 to 1,241 in our case study using scRNA-seq data in autism subjects and controls. We consider 1,500 genes in our analysis. Since the number of genes is larger or comparable to the sample size, we perform penalized estimation. We score each gene's relevance by applying a simple graph kernel smoothing method to each personalized graph. The molecular functions of the top-scored genes are related to autism diseases. For example, a candidate gene RYR2 that encodes protein ryanodine receptor 2 is involved in neurotransmission, a process that is impaired in ASD patients. While our method provides a systemic and unbiased approach to prioritize autism risk genes, the relevance of these genes needs to be further validated in functional studies.

Keywords: Cell dependence; Hurdle model; Poison-LogNormal distribution; Zero-inflation.

Fig. 1

Left: Histograms of the expression proportions of the top 1500 genes which are expressed in most cells in the Velmeshev dataset (individual number 5387) and the Gierahn dataset. Middle: The actual zero proportions versus the expected zero proportions under the fitted PLN model. Right: Histograms of the P-values of Pearson correlation tests for all cell pairs.

Fig. 2

Performance of different graph estimation methods under the non-zero-inflated setting (15). (a) Banded graph: μ₁ = ⋯ =μ_p = 2, c = 2, 157 edges; (b) Hub graph: μ₁ = ⋯ = μ_p = 2, c = 2, 72 edges; (c) Random graph: μ₁ = ⋯ = μ_p = 3, c = 1/2, 79 edges. The x-axes and y-axes represent FPR and TPR respectively. The sparsity of estimated graphs by each method varies by its specific tuning parameter.

Fig. 3

Performance of different graph estimation methods under the zero-inflated setting (16). (a) Banded graph: μ₁ = ⋯ = μ_p = 3, c = 2, γ₀ =1, γ₁ = 0.5, 157 edges; (b) Hub graph: μ₁ = ⋯ = μ_p = 3, c = 2, γ₀ = 0, γ₁ = 0.5, 72 edges; (c) Random graph: μ₁ = ⋯ =μ_p = 3, c = 0.5, γ₀ = 0.5, γ₁ = 0.5, 79 edges. The x-axes and y-axes represent FPR and TPR respectively. The sparsity of estimated graphs by each method varies by its specific tuning parameter.

Fig. 4

Accuracy evaluation of graph estimation with real scRNA-seq datasets based on “catalysis precedes” relationship benchmark (left: with dataset from Velmeshev et al. (2019); right: with dataset from Gierahn et al. (2017))

Fig. 5

Use gene-gene interaction graph estimated from scRNA-seq data of 1,241 neuron cells of individual 5278 (Velmeshev et al., 2019) and autism risk genes by Simons Foundation Autism Research Initiative (SFARI) to prioritize novel autism risk genes. (A) SFARI risk score after neighborhood aggregation by 1-step neighbors (SFARIscore_1s) or 2-step neighbors (SFARIscore_2s). (B) A scatter plot of SFARI risk scores after 1-step versus 2-step neighborhood aggregation. (C-F) Four genes with degree larger than 5 and SFARIscore_1s larger or equal to 1.