Transcriptome outlier analysis implicates schizophrenia susceptibility genes and enriches putatively functional rare genetic variants

Abstract

We searched a gene expression dataset comprised of 634 schizophrenia (SZ) cases and 713 controls for expression outliers (i.e., extreme tails of the distribution of transcript expression values) with SZ cases overrepresented compared with controls. These outlier genes were enriched for brain expression and for genes known to be associated with neurodevelopmental disorders. SZ cases showed higher outlier burden (i.e., total outlier events per subject) than controls for genes within copy number variants (CNVs) associated with SZ or neurodevelopmental disorders. Outlier genes were enriched for CNVs and for rare putative regulatory variants, but this only explained a small proportion of the outlier subjects, highlighting the underlying presence of additional genetic and potentially, epigenetic mechanisms.

Figure 1.

Workflow-integrative analysis of transcriptome (RNAseq) outliers and target resequencing data.

Figure 2.

Frequency of genes containing outlier subjects, and the correlation of the number of such genes with gene expression values. A total of 8355 autosomal genes were expressed in all analyzed subjects. Histograms showed distributions of the number of such genes by the number of outlier subjects in the tails of each gene at 2SD (A) and 3SD (B) cut-offs. The XY plots show no correlation between the number of outlier subjects for a gene and its gene expression value (expressed as mean log₂RPKM) at either 2SD (C) or 3SD (D) cut-offs.

Figure 3.

Copy number variant detection sensitivity of expression outlier analysis. (A) Percentages of detected outliers among genes known to be fully deleted or duplicated [Affymetrix 6.0 SNP array intensity data from the MGS sample (11)], at 2SD or 3SD cut-offs. The data are based on a total of 95 genes known to be deleted (Supplementary Material, Table S1) and 218 genes known to be duplicated (Supplementary Material, Table S2). (B) Expression Z-scores of CNV genes detected as having outlier subjects (2SD cut-off) at lower tails for those with known 1–2 exons deletions or whole gene deletions, and at upper tails for those with known 1–2 exons duplications or known whole gene duplications. Shown are the average expression Z-score ± SD of different groups of genes with outlier subjects: 7 outlier events for genes known to have 1–2 deleted exons, 2 outlier events for genes known to have 1–2 duplicated exons, 80 lower-tail outlier events for genes known to be deleted and 176 upper tail outlier events for genes known to be duplicated. Deletion and duplication calls were made from Affymetrix 6.0 SNP array intensity data for the MGS sample (11). Whole gene deletions showed stronger effects on expression (i.e., larger Z-score) than 1-2 exon deletions (P = 0.04; mean Z-scores of 5.4 versus 3.8; 2-tailed unpaired Student's t-Test) and than whole-gene duplication (P = 1.4 × 10⁻¹⁶, 2-tailed unpaired Student's t-Test). Perhaps due to the small number of events, 1–2 exon duplications did not show significant differences from any other groups.

Figure 4.

Rare to low-frequency regulatory variants are enriched in sequenced outliers of 17 selected genes, and two frameshift variants identified in two ADA expression outliers cause NMD of the transcript. (A) Venn diagram shows the overlap of common variants (MAF ≥ 5%) identified in outlier DNA sequencing and in public whole-genome sequencing data sets (UK10K-TwinsUK and 1000 Genomes). (B) Variants in DNaseI HS and nonsense/frameshift indels are enriched in our outlier-seq data set versus UK10K-TwinsUK date set for rare to low frequency (MAF < 5%) variants, but not for those common variants (MAF ≥ 5%) (C). Fisher's exact test was used and two-tailed P values are presented. 1000 Genomes was not used for comparisons in (B) and (C) because it is a small data set (only 85 CEU subjects). (D) Gene structure of ADA and its exons where the two frameshift mutations (chr20_43249673_CTCTT_del and chr20_43254291_T_ins) reside. (E) The two frameshift mutations were detected by multiplex outlier DNA sequencing (middle panel) and confirmed by Sanger sequencing (bottom panel), but the frameshift allele was absent in RNAseq (top panel). T allele of rs61737144 adjacent to the insertion chr20_43254291_T_ins and in the same haplotype was also not detected at the transcript level. Both findings support that the frameshift variants create premature stop codons that caused NMD of the transcript. (F) Distribution of Z-scored expression of ADA in 312 SZ cases and 322 controls. Two outlier subjects with the frameshift variants are indicated by arrows at the at lower 2SD tail.