Comparison of strategies for identification of regulatory quantitative trait loci of transcript expression traits

Population and phenotypes

We used the Genetic Analysis Workshop 15 (GAW 15) Problem 1 microarray gene expression profiles for the analyses. Data were available for 14 three-generation CEPH Utah families. Expression levels of genes were obtained from lymphoblastoid cells using the Affymetrix Human Focus Arrays [2]. We were provided with data on 3554 transcripts that showed greater variation between individuals than within the same individual.

Family members were genotyped for 2882 autosomal and X-linked single-nucleotide polymorphisms (SNP) generated by the SNP Consortium http://snp.cshl.org/. Genetic map positions were obtained using the SNP Mapping web application developed by the University College Dublin Conway Institute of Biomolecular and Biomedical Research http://actin.ucd.ie/software.html, which uses the Rutgers Combined Linkage-Physical Map of the Human Genome and data taken from the NCBI dbSNP Build 123 (in Kosambi centimorgans). This information was used to calculate multipoint identity by descent matrices (MIBDs) with Merlin and Minx [3], after removal of Mendelian inconsistencies and double recombinants with Preswalk (based on Simwalk mistyping probabilities) [4]. MIBDs were used for linkage analyses.

Transcript distributions were normalized using an inverse normalization transformation of z-scores of individual transcripts regressed on the mean individual transcript level. We further adjusted for the effects of age, age², sex, age × sex and age² × sex interaction using predictive linear regression models in SAS 9.1 (Cary, NC). We generated these residuals as part of our processing of the transcripts for linkage analyses.

Heritability estimation and linkage analysis

Heritability was estimated using maximum likelihood variance decomposition methods in SOLAR [5, 6]. Genome scans were performed using multipoint variance-components models that test for linkage between traits and genetic variants by partitioning the variance of the expression level into its additive genetic and environmental variance components [7]. For transcripts with co-localized QTLs, we performed bivariate linkage analysis to identify shared genetic effects. The bivariate polygenic model estimates correlations caused by residual additive genetic effects (ρ_G) and correlations caused by random environmental effects (ρ_E) [8]. To test for additive genetic correlation among pairs of traits, the log likelihood of a model in which ρ_G is constrained to 0 (null hypothesis, no correlation) or ρ_G = 1 (null hypothesis, complete shared genetic effect) is compared to that of a model in which ρ_G is estimated for the traits. Significant differences among the models (ρ_G ≠ 0) suggest that some of the same genes influence both traits. We also performed linkage analysis using the factors obtained from the PCA in a sample of transcripts with co-localized QTLs.

Principal-component analysis

PCA was used to reduce the number of expression profiles into statistically meaningful groups while retaining the original total variance using all the expression profiles [9, 10]. We selected two different chromosomal regions of a length of 10 to 12 MB in which the QTLs of at least 10 transcripts were co-localized. Only transcripts of genes that were not located in these selected chromosomal regions were included in the analyses (trans-regulatory genes). Because of the small number of individuals in the study and concerns of overfitting the model, a maximum sample of 50 transcript values were considered at one time [11]. The number of factors was determined using the eigenvalue-one requirement [11]. Factors are interpreted by examining the varimax-rotated factor loadings, which are the correlations between each phenotype and the factor in question. Factor loadings greater than or equal to 0.40 in absolute value were used to interpret factors and to characterize the factor structures; this criterion ensures that the individual factor variables share at least 15% of their variance with the given factor [9]. The principle components were obtained by calculating the eigenvalues of the sample covariance matrix, which represent the amount of variance contributed by each factor. Only factors with eigenvalues higher than 1 were considered for linkage analysis.

Integrating data from linkage analysis for gene co-expression

Linkage signals of individual transcript expressions were recorded and the location of QTLs was compared to the location of the transcript gene in order to identify trans-regulatory sites. In addition, the location and LOD scores of QTLs identified in single individual transcript analysis (univariate analysis) were compared with the location of the QTLs identified using bivariate analysis or factors of the PCA. This allowed a determination of whether the bivariate analysis or PCA data reduction analysis improved our evidence for linkage, and if so, a more in-depth examination of the transcripts included in the principal components needs to be examined for biologic interactions on complex disorders.

Chromosome 4q28.2-q31.1 region

We then performed bivariate analysis of all pairwise co-localized transcripts on 4q28.2-q31.1 and found evidence for genetic correlation of co-localized genes, although without much increase in the magnitude of the LOD score (Figure 1). This analysis identified two networks of gene expressions (Figure 1). We obtained two factors using PCA of the first network (Group 1, SNX4, YWHAQ, ASS, MX2, and ISGF3G gene transcripts). Both factors had significant heritability; however, only Factor 1, loading heavily on the MX2 gene, localized to the 4q28.2-q31.1 region (Table 1), and the magnitude of the LOD score was lower than that of the univariate MX2 gene transcript analysis (LOD = 2.28). The heritability of one factor obtained using PCA for Group 2 transcripts was not significant and further analysis was not performed.

Chromosome 13q34 region

We performed analysis in an additional chromosome region of co-localized transcripts, 13q34 region, and noted similar results. Using univariate analysis, 12 transcripts co-localized in this region; and bivariate analysis revealed an intricate network of correlated traits (Table 2 and Figure 2). Using PCA, we obtained five factors, three of them with significant heritability. Similar to our previous findings on chomosome 4, PCA factors did not improve the magnitude of the LOD scores when compared to univariate analysis.

Transcript	Gene Name	Gene Symbol	H2 (SE)a	H2 p-value	LOD score	Trait locus location	Transcript gene locus
200805_at	lectin, mannose-binding 2	LMAN2	0.25(0.10)	0.0003	1.3	13q33.2-q34	5q35.3
209375_at	xeroderma pigmentosum, complementation group C	XPC	0.23 (0.13)	0.007	1.8	13q33.2-q34	3p25
211564_s_at	PDZ and LIM domain 4	PDLIM4	0.20 (0.13)	0.01	2.0	13q34	5q31.1
203366_at	polymerase (DNA directed), gamma	POLG	0.21 (0.11)	0.002	2.0	13q34	15q25
210502_s_at	peptidylprolyl isomerase E (cyclophilin E)	PPIE	0.40 (0.13)	<0.001	1.8	13q34	1p32
217922_at	Mannosidase, alpha, class 1A, member 2	MAN1A2	0.21 (0.10)	0.003	1.6	13q34	1p13
209715_at	chromobox homolog 5 (HP1 alpha homolog, Drosophila)	CBX5	0.28 (0.11)	<0.001	2.4	13q34	12q13.13
203880_at	COX17 homolog, cytochrome c oxidase assembly protein (yeast)	COX17	0.28 (0.11)	0.0004	1.3	13q34	3q13.33
201145_at	HCLS1 associated protein X-1	HAX1	0.33 (0.14)	0.0004	2.0	13q34	1q21.3
201157_s_at	N-myristoyltransferase 1	NMT1	0.31 (0.11)	0.0002	2.0	13q34	17q21.31
209219_at	RD RNA binding protein	RDBP	0.36 (0.11)	<0.001	1.6	13q34	6p21.3
217932_at	mitochondrial ribosomal protein S7	MRPS7	0.11 (0.08)	0.05	1.0	13q34	17q25
Principal component analysis, factor 2 (loading HAX1)b		N/A	0.28 (0.13)	0.001	1.6	13q34	N/A
factor 3 (loading MRPS7)		N/A	0.21 (0.10)	0.003	1.1	13q33.1-33.2	N/A
factor 5 (loading NMT1)		N/A	0.20 (0.12)	0.02	1.9	13q34	N/A

^aH², heritability; se, standard error; N/A, not apply.

^b For principal component analysis, only factors with significant heritability (alpha = 0.05) are shown.