Genomic Convergence Analysis of Schizophrenia: mRNA Sequencing Reveals Altered Synaptic Vesicular Transport in Post-Mortem Cerebellum

doi:10.1371/journal.pone.0003625

PLoS ONE. 2008; 3(11): e3625.

Published online 2008 November 5. doi: 10.1371/journal.pone.0003625

PMCID: PMC2576459

Copyright Mudge et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Genomic Convergence Analysis of Schizophrenia: mRNA Sequencing Reveals Altered Synaptic Vesicular Transport in Post-Mortem Cerebellum

Joann Mudge,¹ Neil A. Miller,¹ Irina Khrebtukova,² Ingrid E. Lindquist,¹ Gregory D. May,¹ Jim J. Huntley,¹ Shujun Luo,² Lu Zhang,² Jennifer C. van Velkinburgh,¹ Andrew D. Farmer,¹ Sharon Lewis,¹ William D. Beavis,¹ Faye D. Schilkey,¹ Selene M. Virk,¹ C. Forrest Black,¹ M. Kathy Myers,¹ Lar C. Mader,¹ Ray J. Langley,¹ John P. Utsey,¹ Ryan W. Kim,¹ Rosalinda C. Roberts,⁵ Sat Kirpal Khalsa,⁴ Meredith Garcia,⁴ Victoria Ambriz-Griffith,⁴ Richard Harlan,⁴ Wendy Czika,⁶ Stanton Martin,⁶ Russell D. Wolfinger,⁶ Nora I. Perrone-Bizzozero,³ Gary P. Schroth,² and Stephen F. Kingsmore¹^*

¹National Center for Genome Resources, Santa Fe, New Mexico, United States of America

²Illumina Inc., Hayward, California, United States of America

³Department of Neurosciences, University of New Mexico, Albuquerque, New Mexico, United States of America

⁴Northern New Mexico College, Española, New Mexico, United States of America

⁵Department of Psychiatry, University of Alabama at Birmingham, Birmingham, Alabama, United States of America

⁶SAS Institute, Cary, North Carolina, United States of America

Rodolfo Aramayo, Editor

Texas A&M University, United States of America

* E-mail: sfk/at/ncgr.org

Conceived and designed the experiments: GPS SFK. Performed the experiments: JM NAM IK IEL GM JJH LS LZ JCvV ADF RJL RWK. Analyzed the data: JM NAM IK IEL GM JJH JCvV ADF SL WB FDS SMV CFB MKM LCM RJL SKK MG VAG REH WC SM RDW NPB GPS SFK. Contributed reagents/materials/analysis tools: NAM GM ADF SL WB FDS SMV CFB MKM LCM JPU RWK RCR WC SM NPB GPS. Wrote the paper: JM SKK MG REH RDW NPB GPS SFK.

Received July 21, 2008; Accepted October 10, 2008.

This article has been cited by other articles in PMC.

Abstract

Schizophrenia (SCZ) is a common, disabling mental illness with high heritability but complex, poorly understood genetic etiology. As the first phase of a genomic convergence analysis of SCZ, we generated 16.7 billion nucleotides of short read, shotgun sequences of cDNA from post-mortem cerebellar cortices of 14 patients and six, matched controls. A rigorous analysis pipeline was developed for analysis of digital gene expression studies. Sequences aligned to approximately 33,200 transcripts in each sample, with average coverage of 450 reads per gene. Following adjustments for confounding clinical, sample and experimental sources of variation, 215 genes differed significantly in expression between cases and controls. Golgi apparatus, vesicular transport, membrane association, Zinc binding and regulation of transcription were over-represented among differentially expressed genes. Twenty three genes with altered expression and involvement in presynaptic vesicular transport, Golgi function and GABAergic neurotransmission define a unifying molecular hypothesis for dysfunction in cerebellar cortex in SCZ.

Introduction

Schizophrenia (SCZ [MIM #181500]) is a common, severe psychiatric disorder with a strong genetic component and complex inheritance [1]–[5]. Cytogenetic, linkage, positional cloning, candidate gene association and genome-wide studies have identified several credible SCZ risk genes. However, few have yet been replicated or translated into causal alleles, in vitro diagnostics or therapeutics [6]–[8]. In many studies of SCZ, genetic analysis has been impeded by phenotypic definition based upon multiple, subjectively ascertained, behavioral parameters that lack reference to biological mechanism [9]–[12]. In addition, more than 20 whole-genome linkage scans have demonstrated heterogeneity of linkage [13], suggesting the existence of genocopies (similar phenotypes that are determined by distinct risk loci) [12]. Evidence exists that the genetic architecture of SCZ may be further obscured by allelic heterogeneity (additive genetic variance segregating in the population at causative loci), epistasis (different combinations of loci producing a phenotype in different pedigrees), pleiotropy (loci that affect more than one phenotype) and phenocopies (heterogeneous environmental factors that mimic allelic effects) [14]. In addition, contributions of risk alleles to complex traits may not fit basal multiplicative and/or threshold models, and studies performed to date may not have had sufficient power, or appropriate theory, to assess non-linear (i.e., epistatic and genotype-by-environment) models. As a result, case-control association studies have identified numerous significantly associated susceptibility loci, but lack of replication among studies is widespread [3], [4], [8], [15]–[18]. Furthermore, the most validated loci were largely selected based on involvement in networks implicated in SCZ (such as dopaminergic and glutamatergic neurotransmission), introducing bias and limiting identification of novel risk factors.

Absence of a clear understanding of the molecular basis of SCZ imposes significant challenges to timely diagnosis and prognostic or therapeutic categorization [12], [19]. Supplementation of diagnostic criteria with biomarkers that are causally related to SCZ or endophenotypes may allow definition of homogeneous subgroups that are predictive of progression and therapeutic response in individual patients [20] and would serve as a starting point for development of therapeutics directed at causal variants. An alternative approach for molecular dissection of SCZ is identification of altered gene expression in affected tissues. Because gene expression reflects both genetic and environmental influences, it may be particularly useful for identifying risk factors for a complex disorder such as SCZ, which is believed to have a multifactorial etiology [21]. Two factors have hitherto limited the effectiveness of gene expression analysis in SCZ: Firstly, mRNA analyses in post-mortem brains in SCZ is challenging due to type I and type II errors resulting from variation in cause of death (affecting agonal gene expression), postmortem interval (affecting RNA quality), concurrent medication, substance abuse, age, sex, race and duration of illness [22]. Secondly, in common with genome-wide association studies, gene expression comparisons employing available cohort sizes sail between the Scylla of many false-positives due to multiple comparisons and the Charybdis of insufficient power to detect true-positives following statistical correction [23]. Recently, however, studies of mRNA expression in post mortem brains in SCZ that account for these variables have started to be reported [7], [21], [22].

An elegant, new approach to navigate Scylla and Charybdis and, thereby, accomplish molecular definition of SCZ may be genomic convergence analysis [24]. Predicated on an implication of the central dogma of molecular biology, genomic convergence analysis posits that clinically relevant nucleotide variation should result in detectable cis- and trans-effects in messenger RNA (mRNA) that amalgamate into functional changes in networks and pathways. Importantly, genomic convergence analysis provides a strategy to collectively interpret and employ the massive, disease-related data sets produced by unbiased (i.e. non-hypothesis driven) linkage and expression studies. Indeed, integration of gene expression and genetic linkage data has shown promise in several neurologic disorders [24]–[27] and has started to be applied to SCZ [7], [28]–[30].

mRNA sequencing with shotgun, massively parallel sequencing platforms has recently shown utility for measurement of transcript abundance, splice isoforms and allelic influence on gene expression [31]–[50]. mRNA abundance is determined by sequencing either 3′ end tags or random cDNA fragments (digital transcript expression, DTE), followed by read alignment to reference databases and calculation of aligned read frequencies. Potential advantages of DTE in comparison to array hybridization include: single molecule sensitivity (corresponding to approximately 1 mRNA molecule per 30 cells; Hayashizake, personal communication); absolute, rather than relative, measurement of transcript abundance; sequence verification for each measurement; comprehensive detection of both known and novel, unannotated transcripts and isoforms; applicability to any eukaryotic species; very little technical imprecision; absence of interference from abundant transcripts (e.g. globin); and extensibility to concomitant measurement of non-coding RNA and to detection of nucleotide and structural variation.

As the first stage of a genomic convergence analysis of SCZ, we describe shotgun mRNA sequencing of an affected tissue (post-mortem cerebellar cortex) of patients and controls, together with a rigorous and systematic approach to DTE. The approach presented is conservative due to application of statistical and bioinformatic methods that substantially reduce type I error rates by using statistical significance criteria rather than fold-change values and by incorporation of potential confounding variables, such as psychotropic medication and cause of death. It is also more comprehensive than prior microarray studies of SCZ. The stringency of this approach, in combination with genome-wide genotypes from the same individuals, is anticipated to enable a genomic convergence analysis that advances understanding of the molecular basis of SCZ.

Results

Sequencing-by-synthesis of cerebellar mRNA

16.7 billion nucleotides of shotgun, full length cDNA sequence data were generated using Illumina Genome Analyzer platforms with sequencing-by-synthesis (SBS) chemistry from 20 mRNA samples [39], [42], [43]. mRNA samples were isolated post-mortem from the lateral hemispheres of the cerebellar cortices of 14 patients with SCZ and 6 control individuals [51], [52] (Table 1). Unrelated subjects were chosen to facilitate sampling of genetic heterogeneity [53], [54]. Cases and controls were approximately matched for age (cases, 45.2±11.8 years; controls 41.3±9.2 years), sex (all male), race, post-mortem interval (cases, 12.2±5.0 hours; controls 17.7±3.3 hours), cause of death, autopsy brain pH (cases, 6.54±0.19; controls 6.46±0.10) and RNA integrity number (cases, 8.06±0.53; controls, 7.87±0.41) (Table 1). 12.5–38.7 million, high quality sequences of length 32–36 bp were generated per sample (Table 2). Sequences were aligned to the human genome and RefSeq transcript databases using the algorithm GMAP, which allowed <2 (≤6%) mismatches [55]. There was little intra-sample variability in the number of sequences aligned to each locus from run-to-run or instrument-to-instrument (Fig. 1A; all source coefficient of variation 3.4%). 43.5±6.7% of sequences aligned to a transcript and 69.4±9.6% to the genome, evidence that annotation of mRNA isoforms in Homo sapiens is incomplete [56] (Table 2). 91% of alignments were unique (Table 2). Reads aligning to more than one location contained repetitive, paralogous, polymorphic or low complexity sequences and primarily mapped to untranslated regions or highly polymorphic gene families, such as major histocompatibility genes [48]. Unmapped sequences did not align to mitochondrial or 1879 viral genomes, offering negative evidence of chronic viral etiology for SCZ in these patients.

Table 1

Clinical Features of Patients and mRNA Samples.

Table 2

cDNA Sequencing and Alignment Statistics.

Figure 1

Characteristics of cDNA sequencing-by-synthesis with Illumina/Solexa instrument.

In order to further investigate differences in proportions of sequences aligning to different reference databases, an additional 2.2 billion nucleotides of shotgun, full length, paired cDNA sequence data were generated from mRNA sample 3S using an Illumina Genome Analyzer II platform and aligned to three reference databases (the human genome, RefSeq transcript and UniGene transcript databases): 46.1% of reads aligned to all three databases, representing known exons. 31.4% of reads aligned to the genome alone, representing novel exons, splice isoforms [57] or transcripts [56]. 5.5% mapped to RefSeq transcript alone, and 6.5% to UniGene transcript alone, representing sequences that span known exon junctions [57]. 10.6% of reads did not align to any of these reference sets, representing reads that span boundaries of novel exons or splice isoforms, novel genomic sequence, or poor quality sequence.

Coverage and Composition of cDNA Sequences

The number of transcripts detected in cerebellar cortex mRNA differed little between samples (33,200±1,000; Table 2), corresponding to 85±3% of RefSeq transcript entries. While 12.5 million sequences per sample was sufficient to reach a plateau in the number of transcripts detected, deeper sequence generation resulted in linear increase in average depth of coverage (Fig. 1D). As anticipated with hexamer-primed cDNA synthesis, the distribution of sequence alignments along transcripts appeared random. Transcript coverage showed a moderate 3′ bias, particularly among long transcripts, attributable to only very slight levels of degradation in these samples together with two rounds of poly-A⁺ RNA selection (Fig. 1B). Decreased coverage was observed at 5′ and 3′ termini due to edge effects of random priming very close (<30 nts) to the ends of the mRNA. No compositional bias was detected. The average depth of coverage achieved was 8.0-fold (Table 2). The distribution of transcript abundance in cerebellar cortex was interesting: 76–81% of expressed transcripts had an abundance of at least 1 read per million. However, as abundance increased, the number of transcripts declined greater than exponentially (Fig. 1C, r²

0.97). Thus, only 0.3–0.5% of expressed transcripts had an abundance of ≥1000 reads per million.

Comparison of Digital Transcript Expression and Array Hybridization

Gene expression was evaluated in mRNA samples from the 20 cerebellar cortices using both oligonucleotide array hybridization, the current standard, and tag frequencies from the Illumina cDNA sequencing assay (aligned reads per million). The dynamic range of read frequencies was two orders of magnitude greater than array hybridization (log₁₀ dynamic range DTE–4.37; log₁₀ dynamic range array hybridization–2.33; Fig. 2). While 41% of arrayed oligonucleotide probes had hybridization intensities greater than the conventional signal: noise threshold (Fig. 2), 85±3% transcripts had aligned reads. Array hybridization signals are transformed and often normalized prior to evaluation of inter-sample differences. For read frequencies, log transformation, but not normalization, improved overlaid kernel density estimates, univariate distribution results, and Mahalanobis distances (Fig. 3 and and4).4). Therefore, log transformed values were used in all subsequent expression analyses. Aligned read frequencies, adjusted for transcript length, correlated weakly with oligonucleotide array hybridization signals (r²

0.35; Fig. 2). Read frequencies exhibited much higher correlation coefficients in pair wise sample comparisons (r²

0.93–0.99) than array hybridization (r²

0.83–0.88; Fig. 5 and and6).6). Furthermore, pair wise sample correlations of genome- and transcript-aligned read frequencies had very similar correlation coefficients (data not shown). SCZ patients could readily be distinguished from controls by unsupervised principal component analysis (by Pearson product-moment correlation, Fig. 7B) or Ward hierarchical clustering of Pearson product-moment correlations of read frequencies (data not shown). In contrast, patients and controls could not be separated with these methods on the basis of array hybridization signals (Fig. 7A and data not shown). Log transformed genome- and transcript-aligned read frequencies showed identical hierarchical clustering and Pearson product-moment correlations of samples (data not shown).

Figure 2

Comparison of gene expression levels measured by two methods.

Figure 3

Comparison of Mahalanobis Distances of Gene Expression by Array Hybridization and Sequence Read Frequencies.

Figure 4

Overlayed kernel density estimates of Gene Expression by Array Hybridization and Sequence Read Frequencies.

Figure 5

Pairwise sample correlations of log10-transformed, genome-aligned read frequencies, showing pairwise correlation coefficients.

Figure 6

Pairwise sample correlations of log10-transformed, array hybridization signals, showing pairwise correlation coefficients.

Figure 7

Unsupervised principal component analysis of log10 transformed array hybridization signals (A) and log10 transformed read frequencies (B).

Patient, sample, and experimental parameters were examined to quantify sources of variability in read-frequency-based DTE. Decomposition of principal components of variance showed that DTE attributed a greater proportion of the total variance to diagnosis (SCZ versus control, 45.3%) than array hybridization (14.1%) (Fig 8). The largest sources of variability of array hybridization results were mRNA quality metrics (Fig. 8A), whereas DTE results were only marginally influenced by mRNA quality. In contrast, large sources of variability in DTE results were year of sequence generation, sequencing instrument-to-instrument variation and cause of death (Fig. 8B). Variance component categories were identical in genome- and transcript-aligned read frequencies (data not shown).

Figure 8

Principal components of variance of log10 transformed array hybridization signals (A) and log10 transformed read frequencies (B).

Digital Transcript Expression Results

Differences between SCZ patients and controls in gene expression in cerebellar cortex were identified with analysis of variance with diagnosis as the discriminatory effect and the major non-diagnosis components of variance as fixed effects (brain pH, post-mortem interval, RNA integrity number and RNA isolation date for array hybridization, and cause of death, post-mortem interval, sequencing instrument and year sequenced for read frequencies). Following FDR correction, no differences were identified between SCZ patients and controls by array hybridization, due to the magnitude of non-diagnosis components of variance. In contrast, 88 genes exhibited FDR-corrected differences in cerebellar expression between SCZ patients and controls in genome-aligned read frequencies and 152 genes differed significantly in transcript-aligned reads (Fig. 9; Table 3). Between these two sets of genes, 25 were identified in common to both genome- and transcript-aligned reads. 24 of these 25 genes showed congruent direction of change in genome- and transcript-aligned reads. 95% of genes exhibiting significant change in expression with only one alignment had congruent but non-significant change in the second alignment or were absent in one of the reference sets. The correlation coefficient of the case-control ratio of the 215 differentially expressed genes between genome- and transcript-aligned reads was 0.73 and between transcript-aligned reads and array hybridization was 0.46 (data not shown). A majority of genes with disparity in the magnitude of gene expression change between alignments had greater read counts for genome- than transcript-alignments, which, in some cases, appeared to be due to unannotated exons and novel splice isoforms (see below).

Figure 9

Volcano plot of analysis of variance of log10-transformed, transcript- (A) and genome-aligned read frequencies (B), showing genes with FDR-corrected differences in LSMeans.

Table 3

Genes Exhibiting Significantly Altered Aligned Read Frequencies between SCZD Cases and Controls.

Genes with the largest increase in expression in SCZ samples were the ISL2 transcription factor, LIM/homeodomain, >800% in genome and transcript alignments; AMME complex, gene 1, 580% in transcript alignments; 3-hydroxy-3-methylglutaryl-Coenzyme A synthase, 372% in transcript alignments; kelch domain containing 1, 360% in transcript alignments; interphotoreceptor matrix proteoglycan 2, 321% in transcript alignments; and NK2 transcription factor related, locus 3, 320% in genome alignments. Genes with the largest decrease in expression in SCZ samples were the solute carrier family 25 (ornithine transporter), member 2, 59% in transcript alignments and LOC126170, similar to peptidylprolyl isomerase A isoform 1, 73% in transcript alignments.

An SCZ candidate gene showing differential expression by both measures was the gamma-aminobutyric acid (GABA)-mediated neuroinhibitory receptor, GABRA1 (increased expression in SCZ by 52% [transcript alignment] and 46% [genome alignment], Table 3). Eleven other genes involved in GABA neurotransmission showed non-significant increases in expression: GABRA2 (59%), GABRA4 (72%), GABRR1 (69%), GABRB1 (168%), GABRG2 (42%), GABRE (35%), GABRB2 (33%), GABBR2 (27%), GABARAPL2 (35%), GABRA6 (14%) and GABRB3 (8%). Four GABAergic genes showed non-significant decreases in expression: GABBR1 (22%), GABRG3 (54%), GABRR2 (26%) and SLC6A1 (GAT-1; 16%).

GO annotation of genes with significantly altered expression revealed over-representation of membrane associated genes (Χ² test, p<0.01), genes involved in Zinc binding or transport (p<0.02), regulation of transcription (p<e⁻⁶), Golgi apparatus (p<e⁻⁶) and vesicle-mediated transport (p<0.01). 13 of 20 genes involved in vesicular transport or Golgi apparatus were up-regulated.

Further annotation of these 20 genes revealed up-regulation of:

Nine genes encoding proteins involved in transport from the trans-Golgi network to the synaptic vesicle (golgi autoantigen, subfamily a, member 1, GOLGA1, 28% in transcript alignments, 14% in genome alignments; solute carrier family 35 (UDP-N-acetylglucosamine transporter), member A3, SLC35A3, 69% in transcript alignments, 18% in genome alignments; component of oligomeric golgi complex 6, COG6, 30% in transcript alignments, 37% in genome alignments; thyroid hormone receptor interactor 11, TRIP11, 48% in transcript alignments, 23% in genome alignments; adaptor-related protein complex 1, gamma 1 subunit, AP1G1, 6% in transcript alignments,−13% in genome alignments; ADP-ribosylation factor guanine nucleotide-exchange factor 2, ARFGEF2, 2% in transcript alignments,7% in genome alignments; vesicle docking protein p115, USO1, 36% in transcript alignments, 25% in genome alignments; Rho-associated, coiled-coil containing protein kinase 1, ROCK1, 85% in transcript alignments, 60% in genome alignments; RAB9B, 35% in transcript alignments, 23% in genome alignments and vacuolar protein sorting 35 homolog,
VPS35 (52% in transcript alignments, 31% in genome alignments), which is involved in retrograde transport to the Golgi apparatus,
Early endosome antigen 1, 162 kD, EEA1 (88% in transcript alignments, −1% in genome alignments), which is involved in homotypic fusion of early endosomes,
Synaptotagmin I, SYT1 (8% in transcript alignments, 15% in genome alignments), which is involved with synaptic vesicle exocytosis, and
AP2 associated kinase 1, AAK1 (64% in transcript alignments, 34% in genome alignments), which is involved in receptor-mediated endocytosis (Fig. 10, Table 3).
Figure 10
Cartoon illustrating functions and/or synaptic locations of 23 proteins corresponding to genes with altered expression in SCZ.

Down-regulated transcripts corresponded to:

Two genes involved in transport from the trans-Golgi network to the synaptic vesicle (syntaxin 10, STX10, −17% in transcript alignments, −19% in genome alignments and ADP-ribosylation factor related protein 1, ARFRP1, −10% in transcript alignments, −12% in genome alignments),
Two genes involved in transport from the endoplasmic reticulum to the Golgi (golgi phosphoprotein 2, GOLM1, −14% in transcript alignments, −17% in genome alignments and synaptic glycoprotein 2, GPSN2, −8% in transcript alignments, −5% in genome alignments),
Two genes involved with synaptic vesicle exocytosis (synaptic vesicle glycoprotein 2A, SV2A, −23% in transcript alignments, −26% in genome alignments and neurochondrin, NCDN, −12% in transcript alignments, −17% in genome alignments), and
Sorting nexin 17, SNX17 (−13% in transcript alignments, −21% in genome alignments) which is involved in retrograde transport back to the Golgi apparatus (Fig. 10, Table 3).

In addition, three post-synaptic membrane genes showed altered expression: GABRA1 (see above), ZACN (ligand-gated ion channel, zinc activated 1; down-regulated 21% in transcript alignments and 5% in genome alignments) and CACNG2 (calcium channel, voltage-dependent, gamma subunit 2; up-regulated 169% in transcript alignments and 10% in genome alignments) (Fig 10, Table 3).

Identification of Novel Splice Isoforms

10,022 genes exhibited ≥20% more reads aligned to genomic loci than to the corresponding, annotated transcripts, evidence for novel, unannotated exons. Putative use of alternative 5′ and 3′ terminal exons, novel internal exons and read-through of introns were all detected [57]. A SCZ susceptibility locus that exhibited intron-read-through in cerebellar mRNA samples was proline oxidase (PRODH, EC 1.5.99.8, SCZ4): 421 cerebellar cortex 41 reads aligned to the sole annotated PRODH transcript, while 838 aligned to the genomic locus. 33 of the latter mapped to intron 14 (nucleotides 17280876–17280950 on Chromosome 22), resulting in insertion of 25 new, in-frame amino acids adjacent to the N-terminus (Fig. 11). In addition to reads mapping within this intron, reads showing splicing out of the intron were also observed (Fig. 11). A further 112 reads aligned to 660 nucleotides of intron 13, generating an alternative N-terminus and premature stop codon (Fig. 11). Presence of these intronic PRODH sequences was observed in all cerebellum mRNA samples and is supported by independently isolated Sanger EST sequences (Genbank accession numbers CN362766, AA300535, AA322439, CD671570, CN362770).

Figure 11

Novel alternative splicing of the PRODH locus.

As noted above, the presence of unannotated exons and novel splice isoforms appeared to explain disparity in the magnitude of gene expression change between cases and controls when using genome- and transcript-alignments. Dedicator of cytokinesis 3 (DOCK3), for example, was significantly up-regulated in SCZ using genome alignments (225 reads/million in cases versus 200 reads per million in controls, −log10(p-value)

3.82) but non-significantly decreased using transcript alignments (138 reads/million in cases versus 147 reads per million in controls). Inspection of genomic read alignments revealed the existence of several, putative, unannotated exons. For example 26 reads in SCZ case 5S and 17 reads in SCZ case 5 aligned uniquely to DOCK3 intronic sequences corresponding to nucleotides 5,130,7158–51307328 on Chromosome 3 (data not shown).

Discussion

As the first stage of a genomic convergence analysis of SCZ, we sought gene expression changes between post-mortem cerebellar cortices of 14 patients and 6 controls using shotgun, clonal mRNA sequencing, together with a rigorous and systematic analysis approach. Cerebellum was chosen for study since it has been shown to be affected in SCZ: Patients exhibit impaired sensory-motor integration, reduced cerebellar volume and grey matter, increased cerebellar blood flow and glucose uptake, both at rest and following cognitive challenge, greatest approximate conditional likelihood score in cerebellar vermis, and gene expression changes in the lateral hemispheres of the cerebellar cortex consistent with decreased GABAergic and increased glutamatergic neurotransmission [32], [51], [52], [58]–[66]. In addition, the simple organization and lack of perturbation by antipsychotic medication of cerebellar cortical circuits relative to other brain regions affected by SCZ was appealing for the current study. Specifically, granule cells, which are important in processing sensory information, are by far the most abundant cerebellar glutamatergic neuron, comprising the majority of neurons in the brain [67].

Cases and controls were all male and were matched reasonably well for age, race, cause of death and RNA quality metrics. Since only cases had received psychotropic medication, this variable could not be matched. There is currently considerable excitement about stratification of patients with SCZ on the basis of clinical or neuroimaging endophenotypes [20]. Unfortunately, that was not possible in this series.

Approximately 22 million random, cDNA sequence tags were generated from each of 20 samples, providing the deepest mRNA coverage reported to date and providing an unparalleled assessment of the transcriptional complexity of the human cerebellar cortex. Sequences did not align to 1879 viral genomes, offering negative evidence of chronic viral etiology for SCZ in these patients. 85% of annotated transcripts were expressed (based on molecular counting, Table 2), of which 63% were detected at a level of at least 1-fold coverage (equivalent to 56 reads per million or 1.1 copies per million). This pervasive transcription accords with other recent studies [32], [44], [56], [68]. In contrast to 3′ end tag sequencing, random-primed mRNA sequencing provided good coverage of coding domains (Fig. 1B), allowing assessment of exon utilization in cerebellar cortex. Alternative splicing is thought to generate more than 5 transcripts per human locus, most of which have not yet been annotated [56]. 85% of loci expressed at a level of at least 1-fold coverage showed evidence for utilization of unannotated exons, based on excess genomic read alignments. Proline oxidase, for example, is a very well studied gene that appeared to utilize hitherto unannotated exons in human cerebellum (Fig. 11). For this reason, sequences were aligned both to a reference transcript database (for estimation of expression of well-annotated transcripts) and to the human genome (for estimation of total expression of all exons, whether annotated or not). With additional improvements in library preparation and analysis software it should be possible to measure digital gene expression on an exon-by-exon basis in a strand-specific manner [31]–[50], [69]. The number of transcripts with aligned reads declined greater than exponentially as a function of level of expression (Fig 1B). As a result, 29 million aligned reads provided an impressive linear dynamic range (4.37×log10, Fig. 2). Despite lower read lengths and quality scores in this study (Table 2) than currently obtained with the Illumina GA II instrument, only 9% of reads mapped non-uniquely. Finally, technical imprecision of DTE measurements was paltry (CV 3.4%). In summary, our technical assessment indicated deep mRNA sequencing to be a powerful tool for digital gene expression and empiric annotation of transcribed elements in genomes, in accord with a rapidly growing literature [31]–[50], [69].

Gratifyingly, digital transcript expression values, expressed as aligned read frequencies, were amenable to very similar data transformation, quality control, pattern discovery, row-by-row modeling and annotation analyses as array hybridization datasets. This enabled direct comparison of digital transcript expression and array hybridization data metrics (Figs. 4–9).9). Like array hybridization signals, aligned read frequencies benefited from log transformation as evidenced by improved Mahalanobis distances, overlaid kernel density estimates, univariate distribution results, and unsupervised principal component analysis. Aligned read frequencies had greater correlation coefficients in pair wise sample comparisons than array hybridization, and much greater ability to distinguish SCZ cases from controls, as evaluated by unsupervised principal component analysis (PCA), Ward hierarchical clustering of Pearson product-moment correlations of read frequencies, and the magnitude of the component of variance attributable to diagnosis. DTE was clearly better than array hybridization for PCA-based segregation of cases and controls. Read frequencies based upon alignment to a reference transcript database and to the human genome were almost indistinguishable by these metrics. Partitioning of variance components of clinical, sample and experimental metadata revealed the largest components of array hybridization variability to be metrics of mRNA quality, in accord with recent studies of post-mortem brain tissue [70]. DTE was much less influenced by RNA quality. Non-diagnosis DTE variability included year of sequence generation, sequencing instrument variation and cause of death (Fig. 8). This is not surprising since substantial improvements occurred in sequencing instrument specifications, reagents and base-calling software during the project. Given continued, rapid evolution of generation II sequencing technologies, it is important to keep instrument specification, reagents and base-calling software relatively stable during mRNA sequencing projects to minimize non-hypothesis-related variability. Importantly, variance decomposition allowed incorporation of large clinical, sample and experimental effects as fixed effects in analyses of variance, minimizing detection of gene expression changes associated with confounding variables. Incorporation of fixed effects in array hybridization, but not DTE, analyses of variance resulted in no expression changes meeting significance cutoffs. Coefficients of correlation between aligned read frequencies, adjusted for transcript length, and array hybridization signals were 0.35 for all values and 0.46 for fold-change in genes with significant differences in expression. These values are somewhat lower than similar comparisons between array hybridization platforms and are probably reflective of the exacting RNA source [71]. It should be noted that aligned read frequencies correlate very well with quantitative PCR results [57], [72]. In conclusion, mRNA sequencing appears to be superior to array hybridization for gene expression analysis and is amenable to standard statistical procedures for quality assessment and gene expression analysis.

Expression of 215 genes differed significantly between cases and controls in post-mortem cerebellar cortex by analysis of variance (Table 3). Of these, 88 were identified in genome-aligned reads and 152 in transcript-aligned reads. Only 25 genes were common to both alignments, but 96% of the 215 had congruent direction of change and correlation between genome- and transcript-aligned fold-change was 0.73. Major contributors to disparity between alignments appeared to be detection of unannotated exons in genome alignments and forced mis-alignments in incomplete transcript datasets. The latter may largely be avoided by removal of non-unique alignments from calculations. As RefSeq transcript becomes more complete, this disparity should decrease. In the interim, it is wise to align both to a well-annotated transcript dataset and to the genome.

17 of the 215 genes with altered cerebellar expression have previously shown changes in SCZ in the dorsolateral prefrontal cortex or superior temporal gyrus (APBA2, BTG1, CACNG2, CAP2, DKFZP434A0131, GABRA1, GOLGA1, HSPBP1, KIAA0256, KPNA1, RPL3, SLC35A3, SLC39A7, SRRM2, TCF4, ZNF148 and ZNF195). Despite substantial differences in gene expression in cerebellar cortex and cerebral cortex [73], 13 genes were congruent with the direction of change reported previously. Three genes involved in GABAergic neurotransmission (GABRA6, GABRB3 and SLC6A1) showed changes in cerebellar cortex in SCZ that agreed with a previous report [52]. In the context of cerebellar cortical function in SCZ, alteration in expression of genes may be causal, consequential, compensatory, or the result of confounding factors. Because the substantial heritability of schizophrenia appears to be polygenic [1]–[5], expression differences might reflect SCZ-associated nucleotide variants that alter expression in cis (e.g. eSNPs). Indeed, seven of these genes have previously shown association with SCZ in genetic association studies (CACNG2, GABRA1, GPSN2, HIRA, PSAP, RANBP5 and TCF4) [74]. Alternatively, expression differences may represent compensatory changes in pathways and networks. GO annotation of genes with altered expression revealed over-representation of membrane association, Zinc binding or transport, regulation of transcription, Golgi apparatus and vesicle-mediated transport.

Most striking were 23 genes involved in presynaptic vesicular transport / Golgi apparatus or post-synaptic neurotransmission, 15 of which were up-regulated and 8 down-regulated. Up-regulated genes included nine involved in transport from the trans-Golgi network to the synaptic vesicle (golgi autoantigen A1, UDP-N-acetylglucosamine transporter (SLC35A3), component of oligomeric golgi complex 6, thyroid hormone receptor interactor 11, adaptor-related protein complex 1G1, ADP-ribosylation factor guanine nucleotide-exchange factor 2, vesicle docking protein p115, Rho-associated, coiled-coil containing protein kinase 1 (ROCK1) and RAB9B), one involved in synaptic vesicle exocytosis (synaptotagmin I), one involved in clathrin-mediated endocytosis (AP2 associated kinase 1 (AAK1)), one involved in homotypic fusion of early endosomes (early endosome antigen 1 (EEA1), and one involved in retrograde transport to the Golgi apparatus (vacuolar protein sorting 35 (VPS35)). Down-regulated genes included two involved in transport from the trans-Golgi network to the synaptic vesicle (syntaxin 10 and ADP-ribosylation factor related protein 1), two involved with synaptic vesicle exocytosis (synaptic vesicle glycoprotein 2A and neurochondrin), two involved in transport from the endoplasmic reticulum to the Golgi (golgi phosphoprotein 2 and synaptic glycoprotein 2) and one involved in retrograde transport to the Golgi apparatus (sorting nexin 17) (Fig. 10, Table 3). Up-regulated post-synaptic membrane genes included GABRA1 (see above) and CACNG2 (calcium channel, voltage-dependent, gamma subunit 2). Several other genes involved in GABAergic neurotransmission showed non-significant changes in expression in cerebellar cortex in SCZ, in agreement with previous reports [52]. One down-regulated post-synaptic membrane gene was ligand-gated ion channel, zinc activated 1. Most of the up-regulated genes, but none of the down-regulated genes in this set, had change in expression of greater than 30%. The largest alterations were CACNG2 (+169%), EEA1 (+88%), ROCK1 (+85%), SLC35A3 (+69%), AAK1 (64%) and GABRA1 and VPS35 (both +52%). Several previous studies have shown decreased expression of genes involved in presynaptic vesicular transport in the dorsolateral prefrontal cortex and hippocampus in SCZ [75]–[78]. While the molecular determinants of synaptic vesicular transport have been intensively studied, it is not yet possible to integrate these transcriptional changes into a single biological outcome. Isolated elevation of AAK1, for example, should decrease neurotransmitter endocytosis [79], while elevated EEA1 should increase early endosomal fusion [80]. However, activity of most of these proteins is dependent either upon post-translational modification or presence of interacting proteins. Furthermore, there is likely to be heterogeneous expression of these proteins among cerebellar cortical neurons. Nevertheless, it appears clear that presynaptic vesicular transport, Golgi function and GABAergic neurotransmission are perturbed in cerebellar cortex in SCZ, in agreement with previous analyses of gene expression in post-mortem brain in SCZ [28], [52], [75] and genetic associations of nucleotide variants in synaptic vesicular transport genes with SCZ (namely CACNG2, COG2, STX1A, SNAP29, MUTED, NRG1, PLDN, NRG3, HOMER3, RTN4R, NRG2, RTN4, CPLX2, BLOC1S1, SYNJ1, HOMER1, CLINT1, SYP, BLOC1S2, SYN3, SYT11, SYNGR1, SNAPAP, STX7, BLOC1S3, DISC1, SNAP25, DLG4, CPLX1, CNO, SYN2, DTNBP1 and DAOA) [74]. Replication of alterations in presynaptic vesicular transport and Golgi function, particularly as affecting GABAergic neurotransmission, in additional cases or other brain regions–such as prefrontal cortex–would substantiate novel molecular mechanisms underpinning SCZ and establish novel targets for therapeutic intervention.

The current study reports the first phase of a genomic convergence analysis of SCZ, which is intended to integrate linkage and expression studies. In a subsequent manuscript, we will present results of expressed nucleotide variants that differ in allele frequency between cases and controls, and integrate nucleotide variant and expression analyses. In addition, a recent study has shown substantial variation in alternative splice isoform expression and alternative polyadenylation in cerebellar cortex between normal individuals [57], and it will be interesting to ascertain whether specific examples of alternative isoform expression show association with SCZ. Confirmation of causal components of the molecular basis of SCZ is anticipated to have significant impact on clinical practice, particularly with respect to timely diagnosis and prognostic or therapeutic categorization [12], [19].

Materials and Methods

Sample preparation

was as previously described [51]. Anonymized cerebellar samples were provided from the Maryland Brain Collection with permission from the Maryland Brain Collection Steering Committee. Permission for the study of brain tissues was provided by families post-mortem in accordance with the guidelines of the Uniform Anatomical Gift Act. Fourteen samples were from patients with a diagnosis of SCZ according to DSM-IV criteria and six were controls. Cortical areas corresponding to crus I/VIIa of cerebellar hemispheres were dissected at −20°C and frozen at −80°C, as described [52]. The average pH of the samples was 6.6±0.17. Cerebellum specimens were obtained according to NIH guidelines for confidentially and privacy. Genomic DNA and total RNA were isolated from samples using standard techniques (Qiagen, Valencia, CA).

Sequencing-by-synthesis

Cerebellum samples were sequenced using Illumina Genome Analyzers with modifications for mRNA samples [39], [42], [43], [52], [57]: Following quality assessment using a Bioanalyzer 2100 (Agilent Inc., Santa Clara, CA; Table 1), poly A+ RNA was isolated from 5–10 µg total RNA by two rounds of oligo-dT selection (Invitrogen Inc., Santa Clara, CA). mRNA was annealed to high concentrations of random hexamers and reverse transcribed. Following second strand synthesis, end repair, and A-tailing, adapters complementary to sequencing primers were ligated to cDNA fragment ends. Resultant cDNA libraries were size fractionated on agarose gels, 200 bp fragments excised, and amplified by 15 cycles of polymerase chain reaction. Following quality assessment using a Bioanalyzer 2100, single-stranded cDNA-adapter fragments were randomly annealed to the surface of a flow cell in a cluster station (Illumina Inc., San Diego, CA) via primers complementary to the adapters and incubated under conditions fostering annealing of the ends of cDNA-adapter fragments to adjacent complementary primers. Primers with annealed cDNA-adapter fragments were extended with DNA polymerase and unlabelled dNTPs in a solid-phase “bridge amplification”. Resultant double-stranded products were denatured, yielding 2 single stranded fragments. The latter steps were repeated for 35 cycles, generating ~40 million clusters of clonal amplicons. Subsequently, 32–36 cycles of sequencing-by-synthesis chemistry were performed in Illumina Genome Analyzer instruments with 4 dNTPs featuring cleavable dyes and reversible terminators. Following each base extension, 4 images (one for each nucleotide) are taken upon laser excitation. Incorporation of the next base occurred after removal of the blocked 3′ terminus and fluorescent tag of the previously incorporated nucleotide. Sequences were retained if of high quality (defined as those passing the default quality filtering parameters used in the Illumina GA Pipeline GERALD stage, i.e. clusters with intensities greater than 0.6-times the average of the highest and the sum of the two highest intensities for the first 12 cycles). Furthermore, sequencing runs with poor average nucleotide quality score graphs (associated with a minority of reads passing the default quality filtering parameters) were discarded.

Read Alignment-based gene expression profiling

High-quality reads were aligned to the human genome, Build 36.2, RefSeq transcript database, Release 22 [81], Unigene Hs, build 215 [82] using the algorithm GMAP [55] and the software system Alpheus® [48], with adjustments for short SBS reads (oligomer overlap interval

3 nt, identities ≥34/36 or 94%). A read was denoted ‘aligned’ to a locus if its sequence alignment to the genomic reference sequence (NCBI build 36.2) fell within the boundaries of the locus coordinates on the chromosome. Locus boundaries on the genome were defined by NCBI annotations, as reported through the Nucleotide database. Annotation data was downloaded on 2/2/2007 using the Batch Entrez query tools (http://www.ncbi.nlm.nih.gov/sites/batchentrez?dbNucleotide). Only the highest scoring alignment(s) was retained. Reads with a single best alignment or with equally good alignments to alternative transcripts of the same locus were considered uniquely aligned. Aligned read frequencies (per million aligned reads) were calculated for each sample and locus using Alpheus® [48].

Array-based gene expression profiling was performed using Affymetrix Human Genome U133 Plus 2.0 oligonucleotide arrays and standard procedures (Affymetrix, Santa Clara, CA). Array-based gene expression profiling of cerebellar cortex of unaffected individuals was performed using Affymetrix Human U95A oligonucleotide arrays as previously reported [83]. Primary image analysis was performed with Affymetrix GENECHIP 3.2. Images were scaled to an average hybridization intensity of 200, and the threshold for expression was an average hybridization intensity of ≥50.

Statistical Analysis

Array hybridization signals and read frequencies were log10 transformed prior to evaluation of inter-sample differences. Overlaid kernel density estimates, univariate distribution results, Mahalanobis distances, correlation coefficients of pair wise sample comparisons, unsupervised principal component analysis (by Pearson product-moment correlation) and Ward hierarchical clustering of Pearson product-moment correlations of read frequencies were performed with JMP Genomics, Version 3.2 (SAS Institute, Cary, NC). Decomposition of principal components of variance of array signals and read frequencies, FDR corrected analysis of variance (with inclusion of the non-diagnosis components of variance discussed above as fixed effects) and chi squared comparisons of GO annotations were also performed with JMP Genomics, Version 3.2. Patient, sample, and experimental parameters examined to quantify sources of variability in measurements were age, alignment percentage, brain pH, cause of death, cluster station, diagnosis, library creator, number of reads, post-mortem interval, psychotropic medication, race, read length, read quality score, RNA isolation date, RNA integrity number, sample storage duration, sequencing instrument and year sequenced.

Acknowledgments

A Deo lumen, ab amicis auxilium. We thank Steven Day, Dan Weems, Melodie J. Rice, Terri L. Gomez, Linda Julien, Kamal Gajendran for assistance with Alpheus software development. We thank Becky Archuleta, Dawn Doyle, Pauline Brito, and Leisha Pak for assistance with variant annotation. We thank Eric Vermaas, Victor Quijano, Juying Yan, Leonda Clendenen and Jimmy Woodward for assistance in generating some of the sequencing data. We thank the NIH Neuroscience DNA Microarray Consortium and Translational Genomics Research Institute (Phoenix, AZ) for help with gene expression profiling of human cerebellar tissue. Irina Khrebtukova, Shujun Luo, Lu Zhang and Gary P. Schroth are employees of Illumina Inc. Study data is available at NCBI Gene Expression Omnibus with accession number GSE12297 and at the Short Read Archive with accession number SRA002355.1 and is freely available for query at http://schizophrenia.ncgr.org.

Footnotes

Competing Interests: Irina Khrebtukova, Shujun Luo, Lu Zhang and Gary P. Schroth are employees of Illumina Inc. Wendy Czika, Stanton Martin and Russell D. Wolfinger are employees of SAS Institute, Cary, NC.

Funding: This work was supported by the National Science Foundation, State of New Mexico, National Institutes of Mental Health, National Center for Research Resources. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Faraone SV, Seidman LJ, Kremen WS, Toomey R, Pepple JR, et al. Neuropsychologic functioning among the nonpsychotic relatives of schizophrenic patients: the effect of genetic loading. Biol Psychiatry. 2000;48:120–126. [PubMed]

2. McGuffin P, Owen MJ, O'Donovan MC, Thapar A, Gottesman London: Geskell Press; 1994.

3. Lichtermann D, Karbe E, Maier W. The genetic epidemiology of schizophrenia and of schizophrenia spectrum disorders. Eur Arch Psychiatry Clin Neurosci. 2000;250:304–310. [PubMed]

4. Ross CA, Margolis RL, Reading SA, Pletnikov M, Coyle JT. Neurobiology of schizophrenia. Neuron. 2006;52:139–153. [PubMed]

5. Saha S, Chant D, Welham J, McGrath J. A systematic review of the prevalence of schizophrenia. PLoS Med. 2005;2:e141. [PMC free article] [PubMed]

6. Abusaad I, MacKay D, Zhao J, Stanford P, Collier DA, et al. Stereological estimation of the total number of neurons in the murine hippocampus using the optical disector. J Comp Neurol. 1999;408:560–566. [PubMed]

7. Harrison PJ, Weinberger DR. Schizophrenia genes, gene expression, and neuropathology: on the matter of their convergence. Mol Psychiatry. 2005;10:40–68; image 45. [PubMed]

8. Riley B, Kendler KS. Molecular genetic studies of schizophrenia. Eur J Hum Genet. 2006;14:669–680. [PubMed]

9. Diagnostic and statistical manual of mental disorders. Washington, DC: American Psychiatric Association; 1994.

10. Cardon LR, Palmer LJ. Population stratification and spurious allelic association. Lancet. 2003;361:598–604. [PubMed]

11. Sillanpaa MJ, Auranen K. Replication in genetic studies of complex traits. Ann Hum Genet. 2004;68:646–657. [PubMed]

12. Tsuang MT, Stone WS, Faraone SV. Toward reformulating the diagnosis of schizophrenia. Am J Psychiatry. 2000;157:1041–1050. [PubMed]

13. Lewis CM, Levinson DF, Wise LH, DeLisi LE, Straub RE, et al. Genome scan meta-analysis of schizophrenia and bipolar disorder, part II: Schizophrenia. Am J Hum Genet. 2003;73:34–48. [PubMed]

14. Rutter M, Moffitt TE, Caspi A. Gene-environment interplay and psychopathology: multiple varieties but real effects. J Child Psychol Psychiatry. 2006;47:226–261. [PubMed]

15. Carter CJ. Schizophrenia susceptibility genes converge on interlinked pathways related to glutamatergic transmission and long-term potentiation, oxidative stress and oligodendrocyte viability. Schizophr Res. 2006;86:1–14. [PubMed]

16. Crespi B, Summers K, Dorus S. Adaptive evolution of genes underlying schizophrenia. Proc Biol Sci. 2007;274:2801–2810. [PMC free article] [PubMed]

17. Kirov G, O'Donovan MC, Owen MJ. Finding schizophrenia genes. J Clin Invest. 2005;115:1440–1448. [PMC free article] [PubMed]

18. Owen MJ, Craddock N, O'Donovan MC. Schizophrenia: genes at last? Trends Genet. 2005;21:518–525. [PubMed]

19. Kendler KS, Neale MC, Walsh D. Evaluating the spectrum concept of schizophrenia in the Roscommon Family Study. Am J Psychiatry. 1995;152:749–754. [PubMed]

20. Gottesman, Gould TD. The endophenotype concept in psychiatry: etymology and strategic intentions. Am J Psychiatry. 2003;160:636–645. [PubMed]

21. Glatt SJ, Everall IP, Kremen WS, Corbeil J, Sasik R, et al. Comparative gene expression analysis of blood and brain provides concurrent validation of SELENBP1 up-regulation in schizophrenia. Proc Natl Acad Sci U S A. 2005;102:15533–15538. [PubMed]

22. Olsen L, Hansen T, Jakobsen KD, Djurovic S, Melle I, et al. The estrogen hypothesis of schizophrenia implicates glucose metabolism: association study in three independent samples. BMC Med Genet. 2008;9:39. [PMC free article] [PubMed]

23. Cheng C, Pounds S. False discovery rate paradigms for statistical analyses of microarray gene expression data. Bioinformation. 2007;1:436–446. [PMC free article] [PubMed]

24. Hauser MA, Li YJ, Takeuchi S, Walters R, Noureddine M, et al. Genomic convergence: identifying candidate genes for Parkinson's disease by combining serial analysis of gene expression and genetic linkage. Hum Mol Genet. 2003;12:671–677. [PubMed]

25. Li YJ, Oliveira SA, Xu P, Martin ER, Stenger JE, et al. Glutathione S-transferase omega-1 modifies age-at-onset of Alzheimer disease and Parkinson disease. Hum Mol Genet. 2003;12:3259–3267. [PubMed]

26. Noureddine MA, Li YJ, van der Walt JM, Walters R, Jewett RM, et al. Genomic convergence to identify candidate genes for Parkinson disease: SAGE analysis of the substantia nigra. Mov Disord. 2005;20:1299–1309. [PubMed]

27. Oliveira SA, Li YJ, Noureddine MA, Zuchner S, Qin X, et al. Identification of risk and age-at-onset genes on chromosome 1p in Parkinson disease. Am J Hum Genet. 2005;77:252–264. [PubMed]

28. Le-Niculescu H, Balaraman Y, Patel S, Tan J, Sidhu K, et al. Towards understanding the schizophrenia code: an expanded convergent functional genomics approach. Am J Med Genet B Neuropsychiatr Genet. 2007;144B:129–158. [PubMed]

29. Lipska BK, Mitkus SN, Mathew SV, Fatula R, Hyde TM, et al. Functional genomics in postmortem human brain: abnormalities in a DISC1 molecular pathway in schizophrenia. Dialogues Clin Neurosci. 2006;8:353–357. [PMC free article] [PubMed]

30. Middleton FA, Pato CN, Gentile KL, McGann L, Brown AM, et al. Gene expression analysis of peripheral blood leukocytes from discordant sib-pairs with schizophrenia and bipolar disorder reveals points of convergence between genetic and functional genomic approaches. Am J Med Genet B Neuropsychiatr Genet. 2005;136B:12–25. [PubMed]

31. Bainbridge MN, Warren RL, Hirst M, Romanuik T, Zeng T, et al. Analysis of the prostate cancer cell line LNCaP transcriptome using a sequencing-by-synthesis approach. BMC Genomics. 2006;7:246. [PMC free article] [PubMed]

32. Kim JB, Porreca GJ, Song L, Greenway SC, Gorham JM, et al. Polony multiplex analysis of gene expression (PMAGE) in mouse hypertrophic cardiomyopathy. Science. 2007;316:1481–1484. [PubMed]

33. Butz JA, Roberts KG, Edwards JS. Detecting changes in the relative expression of KRAS2 splice variants using polymerase colonies. Biotechnol Prog. 2004;20:1836–1839. [PubMed]

34. Butz JA, Yan H, Mikkilineni V, Edwards JS. Detection of allelic variations of human gene expression by polymerase colonies. BMC Genet. 2004;5:3. [PMC free article] [PubMed]

35. Cheung F, Haas BJ, Goldberg SM, May GD, Xiao Y, et al. Sequencing Medicago truncatula expressed sequenced tags using 454 Life Sciences technology. BMC Genomics. 2006;7:272. [PMC free article] [PubMed]

36. Chiu KP, Ariyaratne P, Xu H, Tan A, Ng P, et al. Pathway aberrations of murine melanoma cells observed in Paired-End diTag transcriptomes. BMC Cancer. 2007;7:109. [PMC free article] [PubMed]

37. Emrich SJ, Barbazuk WB, Li L, Schnable PS. Gene discovery and annotation using LCM-454 transcriptome sequencing. Genome Res. 2007;17:69–73. [PubMed]

38. Jones-Rhoades MW, Borevitz JO, Preuss D. Genome-wide expression profiling of the Arabidopsis female gametophyte identifies families of small, secreted proteins. PLoS Genet. 2007;3:1848–1861. [PubMed]

39. Lister R, O'Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, et al. Highly integrated single-base resolution maps of the epigenome in Arabidopsis. Cell. 2008;133:523–536. [PMC free article] [PubMed]

40. Marioni J, Mason C, Mane S, Stephens M, Gilad Y. RNA-seq: An assessment of technical reproducibility and comparison with gene expression arrays. Genome Res 2008 [PubMed]

41. Mikkilineni V, Mitra RD, Merritt J, DiTonno JR, Church GM, et al. Digital quantitative measurements of gene expression. Biotechnol Bioeng. 2004;86:117–124. [PubMed]

42. Morin RD, O'Connor MD, Griffith M, Kuchenbauer F, Delaney A, et al. Application of massively parallel sequencing to microRNA profiling and discovery in human embryonic stem cells. Genome Res. 2008;18:610–621. [PubMed]

43. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5:621–628. [PubMed]

44. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, et al. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science. 2008;320:1344–1349. [PMC free article] [PubMed]

45. Ng P, Tan JJ, Ooi HS, Lee YL, Chiu KP, et al. Multiplex sequencing of paired-end ditags (MS-PET): a strategy for the ultra-high-throughput analysis of transcriptomes and genomes. Nucleic Acids Res. 2006;34:e84. [PMC free article] [PubMed]

46. Nielsen KL, Hogh AL, Emmersen J. DeepSAGE–digital transcriptomics with high sensitivity, simple experimental protocol and multiplexing of samples. Nucleic Acids Res. 2006;34:e133. [PMC free article] [PubMed]

47. Ohtsu K, Smith MB, Emrich SJ, Borsuk LA, Zhou R, et al. Global gene expression analysis of the shoot apical meristem of maize (Zea mays L.). Plant J. 2007;52:391–404. [PMC free article] [PubMed]

48. Sugarbaker DJ, Richards WG, Gordon GJ, Dong L, De Rienzo A, et al. Transcriptome sequencing of malignant pleural mesothelioma tumors. Proc Natl Acad Sci U S A. 2008;105:3521–3526. [PubMed]

49. Toth AL, Varala K, Newman TC, Miguez FE, Hutchison SK, et al. Wasp gene expression supports an evolutionary link between maternal behavior and eusociality. Science. 2007;318:441–444. [PubMed]

50. Weber AP, Weber KL, Carr K, Wilkerson C, Ohlrogge JB. Sampling the Arabidopsis transcriptome with massively parallel pyrosequencing. Plant Physiol. 2007;144:32–42. [PubMed]

51. Paz RD, Andreasen NC, Daoud SZ, Conley R, Roberts R, et al. Increased expression of activity-dependent genes in cerebellar glutamatergic neurons of patients with schizophrenia. Am J Psychiatry. 2006;163:1829–1831. [PubMed]

52. Bullock WM, Cardon K, Bustillo J, Roberts R, Perrone-Bizzozero NI. Altered expression of genes involved in GABAergic transmission and neuromodulation of granule cell activity in the cerebellum of patients with schizophrenia. Am J Psychiatry In press 2008 [PubMed]

53. Freimer N, Sabatti C. The use of pedigree, sib-pair and association studies of common diseases for genetic mapping and epidemiology. Nat Genet. 2004;36:1045–1051. [PubMed]

54. McClellan JM, Susser E, King MC. Schizophrenia: a common disease caused by multiple rare alleles. Br J Psychiatry. 2007;190:194–199. [PubMed]

55. Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005;21:1859–1875. [PubMed]

56. Birney E, Stamatoyannopoulos JA, Dutta A, Guigo R, Gingeras TR, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature. 2007;447:799–816. [PMC free article] [PubMed]

57. Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, et al. Alternative Isoform Regulation in Human Tissue Transcriptomes. 2008. Nature Oct 22 Advance Online Publication.

58. Andreasen NC, Calage CA, O'Leary DS. Theory of mind and schizophrenia: a positron emission tomography study of medication-free patients. Schizophr Bull. 2008;34:708–719. [PMC free article] [PubMed]

59. Andreasen NC, Pierson R. The role of the cerebellum in schizophrenia. Biol Psychiatry. 2008;64:81–88. [PMC free article] [PubMed]

60. Bottmer C, Bachmann S, Pantel J, Essig M, Amann M, et al. Reduced cerebellar volume and neurological soft signs in first-episode schizophrenia. Psychiatry Res. 2005;140:239–250. [PubMed]

61. Edwards CR, Newman S, Bismark A, Skosnik PD, O'Donnell BF, et al. Cerebellum volume and eyeblink conditioning in schizophrenia. Psychiatry Res. 2008;162:185–194. [PMC free article] [PubMed]

62. Ho BC, Mola C, Andreasen NC. Cerebellar dysfunction in neuroleptic naive schizophrenia patients: clinical, cognitive, and neuroanatomic correlates of cerebellar neurologic signs. Biol Psychiatry. 2004;55:1146–1153. [PubMed]

63. Malaspina D, Harkavy-Friedman J, Corcoran C, Mujica-Parodi L, Printz D, et al. Resting neural activity distinguishes subgroups of schizophrenia patients. Biol Psychiatry. 2004;56:931–937. [PMC free article] [PubMed]

64. Okugawa G, Nobuhara K, Takase K, Kinoshita T. Cerebellar posterior superior vermis and cognitive cluster scores in drug-naive patients with first-episode schizophrenia. Neuropsychobiology. 2007;56:216–219. [PubMed]

65. Potkin SG, Alva G, Fleming K, Anand R, Keator D, et al. A PET study of the pathophysiology of negative symptoms in schizophrenia. Positron emission tomography. Am J Psychiatry. 2002;159:227–237. [PubMed]

66. Puri BK, Counsell SJ, Saeed N, Bustos MG, Treasaden IH, et al. Regional grey matter volumetric changes in forensic schizophrenia patients: an MRI study comparing the brain structure of patients who have seriously and violently offended with that of patients who have not. BMC Psychiatry. 2008;8(Suppl 1):S6. [PMC free article] [PubMed]

67. Chadderton P, Margrie TW, Hausser M. Integration of quanta in cerebellar granule cells during sensory processing. Nature. 2004;428:856–860. [PubMed]

68. Kapranov P, Cheng J, Dike S, Nix DA, Duttagupta R, et al. RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science. 2007;316:1484–1488. [PubMed]

69. de Hoon M, Hayashizaki Y. Deep cap analysis gene expression (CAGE): genome-wide identification of promoters, quantification of their expression, and network inference. Biotechniques. 2008;44:630, 632. [PubMed]

70. Popova T, Mennerich D, Weith A, Quast K. Effect of RNA quality on transcript intensity levels in microarray analysis of human post-mortem brain tissues. BMC Genomics. 2008;9:91. [PMC free article] [PubMed]

71. Chen JJ, Hsueh HM, Delongchamp RR, Lin CJ, Tsai CA. Reproducibility of microarray data: a further analysis of microarray quality control (MAQC) data. BMC Bioinformatics. 2007;8:412. [PMC free article] [PubMed]

72. Canales RD, Luo Y, Willey JC, Austermiller B, Barbacioru CC, et al. Evaluation of DNA microarray results with quantitative gene expression platforms. Nat Biotechnol. 2006;24:1115–1122. [PubMed]

73. Evans SJ, Choudary PV, Vawter MP, Li J, Meador-Woodruff JH, et al. DNA microarray analysis of functionally discrete human brain regions reveals divergent transcriptional profiles. Neurobiol Dis. 2003;14:240–250. [PMC free article] [PubMed]

74. Allen NC, Bagade S, McQueen MB, Ioannidis JP, Kavvoura FK, et al. Systematic meta-analyses and field synopsis of genetic association studies in schizophrenia: the SzGene database. Nat Genet. 2008;40:827–834. [PubMed]

75. Mirnics K, Middleton FA, Marquez A, Lewis DA, Levitt P. Molecular characterization of schizophrenia viewed by microarray analysis of gene expression in prefrontal cortex. Neuron. 2000;28:53–67. [PubMed]

76. Thompson PM, Egbufoama S, Vawter MP. SNAP-25 reduction in the hippocampus of patients with schizophrenia. Prog Neuropsychopharmacol Biol Psychiatry. 2003;27:411–417. [PubMed]

77. Vawter MP, Crook JM, Hyde TM, Kleinman JE, Weinberger DR, et al. Microarray analysis of gene expression in the prefrontal cortex in schizophrenia: a preliminary study. Schizophr Res. 2002;58:11–20. [PubMed]

78. Vawter MP, Thatcher L, Usen N, Hyde TM, Kleinman JE, et al. Reduction of synapsin in the hippocampus of patients with bipolar disorder and schizophrenia. Mol Psychiatry. 2002;7:571–578. [PubMed]

79. Henderson DM, Conner SD. A novel AAK1 splice variant functions at multiple steps of the endocytic pathway. Mol Biol Cell. 2007;18:2698–2706. [PMC free article] [PubMed]

80. Haas AK, Fuchs E, Kopajtich R, Barr FA. A GTPase-activating protein controls Rab5 function in endocytic trafficking. Nat Cell Biol. 2005;7:887–893. [PubMed]

81. Pruitt KD, Tatusova T, Maglott DR. NCBI Reference Sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res. 2005;33:D501–504. [PMC free article] [PubMed]

82. Pontius JU, Wagner L, Schuler GD. UniGene: a unified view of the transcriptome. The NCBI Handbook. Bethesda, MD: National Center for Biotechnology Information; 2003.

83. Su AI, Cooke MP, Ching KA, Hakak Y, Walker JR, et al. Large-scale analysis of the human and mouse transcriptomes. Proc Natl Acad Sci U S A. 2002;99:4465–4470. [PubMed]

Articles from PLoS ONE are provided here courtesy of
Public Library of Science