(Sigma Chemical HPA014769, 1/20,000), Actin (Sigma Chemical A2066, 1/5,000), Acox1 (Abnova PAB4367, 1/2000). Data filtering For transcript data we applied three filtering steps based on 1) genetic heritability, 2) probeset annotation. We have profiled 3 mice per strain which allowed us to estimate the broad sense heritability for each transcript. Broad sense heritability for each transcript was measured using ANOVA where strain information was used as a grouping factor. The broad sense heritability which is defined as the ratio of genetic variance over total variance for the phenotype was estimated by dividing the sum of squares of the strain information factor over total sum of squares in the ANOVA. The significance of heritability was established if the p- value for the strain information term in ANOVA was below the nominal 0.05 threshold. The selection cutoff for including gene in the analysis based on heritability was set to heritability p-value of 0.05. From the 22700 probesets 10186 probesets did not meet this cutoff. For annotation filtering, we acquired the Ensembl Gene ID for each Affymetrix probeset and selected those probesets that were only annotated to only one Ensembl gene. From the initial 22700 probesets, 4401 probesets had ambiguous annotation (either did not map to a gene or mapped to more than 2 Ensembl genes). Overall, 9896 probesets met both filtering criteria (significant heritability and unique Ensembl annotation). For the protein data, the initial filtering steps were based on 1) eliminating peptides with excessive missing values which would otherwise have unreliable mapping information, 2) eliminating peptides with missed internal cleavage sites which cause unreliable measurement. To annotate peptides we utilized the SpliceCenter web-based tool [42] to obtain the location of the exon each peptide represents. Peptides which mapped to multiple exons of more than one gene (as determined by SpliceCenter) were excluded from the analysis because of ambiguous annota- tion. The genomic exon coordinates for each peptide was then used to query the Ensembl database to acquire the Ensembl Gene IDs. In the second step of filtering peptides that had more than the value of 2 for the ratio of total variance over technical variance were chosen. Overall, 1543 peptides met the two-stage filtering described above. Genotyping and genome-wide association mapping Inbred strains were previously genotyped by the Broad Institute (http://www.broadinstitute.org/mouse/hapmap ), and they were combined with the genotypes from Wellcome Trust Center for Human Genetics (WTCHG). Genotypes of RI strains at the Broad SNPs were inferred from WTCHG genotypes by interpolating alleles at polymorphic SNPs among parental strains, calling ambiguous genotypes missing. Of the 140,000 SNPs available, 95,854 were informative with an allele frequency greater than 10% and missing values in less than 10% of the strains. These SNPs were used for both protein and transcript genome-wide association analysis. We applied the following linear mixed model to account for the population structure and genetic relatedness among strains in the genome-wide association mapping [27]: y = m+xb+u+e. In the formula, m represents mean, x represents SNP effect, u represents random effects due to genetic relatedness with Var(u) = sg2K and Var(e) = se2, where K represents IBS (identi- ty-by-state) matrix across all genotypes in the HMDP panel. A restricted maximum likelihood (REML) estimate of sg2 and se2 are computed using EMMA, and the association mapping is performed based on the estimated variance component with a standard F-test to test b?0. We applied EMMA (Efficient Mixed Global Mouse mRNA and Protein Levels Comparison Model Association) as an R implementation of a linear mixed model. The percent of variance explained for each molecular phenotype was calculated using the SNP effect calculated from EMMA by defining it as 1-(variance of residuals/variance of original pheno- types). It should be noted that since EMMA is orders of magnitude faster than other implementations commonly used, we were able to perform statistical analyses for all pairs of transcripts and genome wide markers in a few hours using a cluster of 50 processors. Both pQTL and eQTL were defined as ‘‘local’’ if the peak association SNP position was within a 4 Mb interval, flanking 2 Mb on either side of the transcription start and end of the gene under regulation. Genome-wide cutoff: Genome-wide cutoffs were calculated as the false discovery rates using the ‘‘qvalue’’ package for FDR calculation in the R statistical software [43]. Due to the computational complexity associated with evaluating q-values for over 400 million p-values, we computed the FDRs by taking the average FDR for 100 samples each containing 5 million randomly selected p-values from the original calculated p-values. FDR calculation was carried out separately for the protein and transcript dataset. GO analysis and other statistical methods/software All statistical analyses and data visualizations were carried out using the R statistical software (available at http://cran.r-project. org/). Classification of proteins and transcripts to various GO categories was accomplished using Mouse Genome Informatics website at Jackson Laboratories and the GO ontology tool at Princeton. For each probe and each peptide, we first obtained the MGI IDs using the MGI batch query tool at http://www. informatics.jax.org/ [44]. Using MGI IDs we utilized the GO Term Mapper at http://go.princeton.edu, which is based on map2slim algorithm [45] to obtain the GO annotations and summary statistics. The background geneset used in this analysis was the list of all genes annotated by MGI. To assess the significance of the correlation coefficients observed in the PT-pair GO analysis, for each GO category we created 100,000 bootstrap datasets each equal in size to the number of genes assigned to the GO term. Bootstrapping was carried out randomly and without replacement from the pool of 584 original correlation p-values among the PT- pairs. The significance of the observed average p-values for each GO term is reported as the two-tailed test against the empirical distribution created by the corresponding 100,000 permutation set. All the correlation coefficients and corresponding p-values reported in the paper are calculated using the bicor function in the WGCNA R package [46]. The main advantage of using bicor, which performs biweight midcorrelation calculation, over Pearson’s correlation is based the robustness of the correlation coefficient measurement to the presence of outliers in the data. Supporting Information Dataset S1 Peptide expression values as measured by the LCMS. Found at: doi:10.1371/journal.pgen.1001393.s001 (TXT) Dataset S2 Peptide annotations. Found at: doi:10.1371/journal.pgen.1001393.s002 (TXT) Dataset S3 SNP information for the Affymetrix probes. Found at: doi:10.1371/journal.pgen.1001393.s003 (XLS) Figure S1 Immunoblotting validation results. Found at: doi:10.1371/journal.pgen.1001393.s004 (DOC) Figure S2 Cellular compartment representation of the measured proteins and transcripts. Found at: doi:10.1371/journal.pgen.1001393.s005 (TIF) PLoS Genetics | www.plosgenetics.org 15 June 2011 | Volume 7 | Issue 6 | e1001393