Computational protocol: An Immune Response Network Associated with Blood Lipid Levels

Similar protocols

Protocol publication

[…] After each expression array was scanned, background corrected probe signal intensities and bead counts were outputted from Illumina's BeadStudio software in order to undergo further processing. Strip-level quantile normalization was then used to force probe intensity distributions for all samples on all arrays to be the same . Since each sample was technically replicated, the normalized values were then used to measure their correlation via Pearson's product moment correlation coefficient (Ρ) and Spearman's rank correlation coefficient (ρ). Generally, reproducibility was high (). To further assess data quality, we also generated MA plots between replicate arrays after normalization . We manually inspected each sample's MA plot for curvature or overt deviation from the M = 0 axis, none exhibited these characteristics. A sample was removed from further analysis if its Ρ was <0.94 or ρ was <0.60 (9 samples fail).To combine raw signal intensities from corresponding replicates, the signals (S) were weighted by the number of beads (b) contributing to each signal and summed to arrive at one measure of signal intensity (δ) for each sample at each probe:Probes that did not meet certain criteria were removed from further analysis: (a) non-autosomal (b) complementary to cDNA from erythrocyte globin components (c) map to more than one genomic position.For each genotyping array, Cy3 and Cy5 signal intensities were exported from BeadStudio and pooled together for clustering with the Illuminus genotype calling algorithm . Samples were removed from further analysis if they showed low quality (call rate <0.95, 19 samples removed), failed to match Sequenom genotype fingerprinting (concordance <0.90 for at least 10 genotypes, 0 samples removed), or were a previously unknown close relation or duplication (pairwise identity by descent pi-hat >0.10, 1 sample removed). SNPs failing to meet the following quality thresholds were also removed from further analysis: call rate >0.95, minor allele frequency >0.01, and Hardy-Weinberg equilibrium P value >1.0×10−6. 37,558 SNPs were removed in total.Un-observed SNPs were imputed with the software IMPUTE version 5 using phased HapMap release 22 haplotypes from the CEU panel . A genotype was assigned if its posterior probability was >0.95 or missing if not, and all SNPs underwent the same filtering as those above. 249,345 SNPs were removed in total, leaving 2,061,516 SNPs for further analysis. [...] To control for structure in the Finnish population, we used principal components analysis (PCA) on the genotypic data in order to identify outliers who descend from outside the Helsinki region (). All SNPs underwent PCA with the EIGENSOFT software ; regression residuals involving the 2 previous SNPs were used as inputs to correct for SNP linkage disequilibrium. Samples exceeding eight standard deviations along any statistically significant principal component were removed from further analysis (, 17 samples removed). A principal component was considered significant if its Tracy-Widom P value was <0.01. [...] Network analysis was done using the R packages, WGCNA , , and NEO .The undirected transcription network considered the top 10% of expression signals for meta-lipids (3,520 unique signals). The correlation matrix was constructed via all against all Pearson correlation coefficient calculations and the adjacency matrix was calculated with a soft threshold power of nine (). Genes were then hierarchically clustered and visualized in a dendrogram (), where a ‘leaf’ constitutes an individual gene and ‘branches’ are clusters of tightly correlated genes. The dynamic tree cut function in WGCNA with a minimum module size of 10 genes was used to determine initial modules. Individual module expression profiles underwent singular value decomposition and the summary module profiles from the vector corresponding to the first singular value were clustered to identify modules that were highly correlated (those less than a dendrogram height of 0.20). These modules were then merged.To correlate module summary profiles with lipid traits, a t-test of Spearman's rank correlation was used. The corresponding Spearman correlation coefficients and P values can be observed in . Statistical significance was determined by estimating the number of co-expression modules in the entire dataset. Given the 23 modules calculated from 1000 expression probes, we estimated the total number of modules to be (23×35419/1000) = 814.637. Therefore, the appropriate alpha level was determined to be (0.05/814.637) = 6.14×10−5. Calculations of module membership and individual gene significance () have been previously defined . Only module K (the Lipid Leukocyte, LL, module) was used in further analyses.NEO was used to predict the directedness of the network using causal SNPs as anchors. Of the lipid traits associated with LL module expression, HDL and TG were selected because the genetic variation underlying them has been studied extensively. Since the choice of SNPs can have a large impact on the directedness of the network (non-causal SNPs can introduce noise) and the DILGOM dataset (N = 518) is not sufficiently powered to significantly detect many of the known variants, we use only the strongest signals from recent genome-wide association studies , ; rs3764261 (CETP) was included for HDL and rs1260326 (GCKR) was included for TG. In our dataset, the strongest signal previously found for TG, rs964184 (APOA1-C3-A4-A5), did not pass quality control filters. Since rs2251746 has been shown to be an eQTL for FCER1A and LL module expression, we also include it as a causal anchor. To further verify that these loci can be considered causal anchors in the DILGOM dataset, we adopt the automatic SNP selection approach in NEO using both a greedy method and forward-stepwise regression . We observed that all SNPs were correctly assigned to their respective nodes. An edge exists if the edge score (the absolute value of the Pearson correlation between nodes A and B) exceeds a threshold of 0.3. Since all nodes have a causal anchor, the NEO score (the log10 ratio of a fitted causal model P value to the next best causal model P value) defined in the main text is the NEO.NB.OCA score. An edge is considered significantly oriented if the NEO score exceeds a threshold of 0.3. Simulation studies have shown that a NEO.NB.OCA score of 0.3 or more corresponds to a false positive rate of 5% or less (cite NEO). We further considered the path coefficient for A→B (Z test statistic >1.96 or <−1.96) and, since the NEO score is a ratio of model P values, the fit of the primary model MA→A→B←MB (P value should be >0.05). See for directed network edge statistics. […]

Pipeline specifications

Software tools Illuminus, IMPUTE, EIGENSOFT, WGCNA, Dynamic Tree Cut
Applications Population genetic analysis, GWAS
Organisms Homo sapiens
Diseases Drug Hypersensitivity, Hypersensitivity, Atherosclerosis