Heritability Analysis
As an application of the variance component model, this note demonstrates the workflow for heritability analysis in genetics, using a sample data set cg10k with 6,670 individuals and 630,860 SNPs. Person IDs and phenotype names are masked for privacy. cg10k.bed, cg10k.bim, and cg10k.fam is a set of Plink files in binary format. cg10k_traits.txt contains 13 phenotypes of the 6,670 individuals.
;ls cg10k.bed cg10k.bim cg10k.fam cg10k_traits.txt
cg10k.bed cg10k.bim cg10k.fam cg10k_traits.txt
Machine information:
versioninfo()
Julia Version 0.6.0 Commit 903644385b (2017-06-19 13:05 UTC) Platform Info: OS: macOS (x86_64-apple-darwin13.4.0) CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz WORD_SIZE: 64 BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell) LAPACK: libopenblas64_ LIBM: libopenlibm LLVM: libLLVM-3.9.1 (ORCJIT, haswell)
Read in binary SNP data
We will use the SnpArrays.jl package to read in binary SNP data and compute the empirical kinship matrix. Issue
Pkg.clone("https://github.com/OpenMendel/SnpArrays.jl.git")
within Julia to install the SnpArrays package.
using SnpArrays
# read in genotype data from Plink binary file (~50 secs on my laptop) @time cg10k = SnpArray("cg10k")
22.902730 seconds (51.62 k allocations: 1005.845 MiB, 0.11% gc time)
6670×630860 SnpArrays.SnpArray{2}:
(false, true) (false, true) … (true, true) (true, true)
(true, true) (true, true) (false, true) (true, false)
(true, true) (true, true) (true, true) (true, true)
(true, true) (true, true) (false, true) (true, true)
(true, true) (true, true) (true, true) (false, true)
(false, true) (false, true) … (true, true) (true, true)
(false, false) (false, false) (true, true) (true, true)
(true, true) (true, true) (true, true) (false, true)
(true, true) (true, true) (true, true) (true, true)
(true, true) (true, true) (false, true) (true, true)
(true, true) (true, true) … (true, true) (true, true)
(false, true) (false, true) (true, true) (false, true)
(true, true) (true, true) (true, true) (false, true)
⋮ ⋱
(false, true) (false, true) (false, true) (false, true)
(false, true) (false, true) (false, true) (true, true)
(true, true) (true, true) … (false, true) (true, true)
(false, true) (false, true) (true, true) (false, true)
(true, true) (true, true) (false, true) (true, true)
(true, true) (true, true) (false, false) (false, true)
(true, true) (true, true) (true, true) (false, true)
(true, true) (true, true) … (true, true) (true, true)
(true, true) (true, true) (false, true) (true, true)
(true, true) (true, true) (true, true) (false, true)
(false, true) (false, true) (true, true) (true, true)
(true, true) (true, true) (true, true) (true, true)
Summary statistics of SNP data
people, snps = size(cg10k)
(6670, 630860)
# summary statistics (~50 secs on my laptop) @time maf, _, missings_by_snp, = summarize(cg10k);
24
# 5 number summary and average MAF (minor allele frequencies) quantile(maf, [0.0 .25 .5 .75 1.0]), mean(maf)
([0.00841726 0.124063 … 0.364253 0.5], 0.24536516625042462)
# Pkg.add("Plots") # Pkg.add("PyPlot") using Plots pyplot() histogram(maf, xlab = "Minor Allele Frequency (MAF)", label = "MAF")
# proportion of missing genotypes sum(missings_by_snp) / length(cg10k)
0.0013128198764010824
# proportion of rare SNPs with maf < 0.05 countnz(maf .< 0.05) / length(maf)
0.07228069619249913
Empirical kinship matrix
We estimate empirical kinship based on all SNPs by the genetic relation matrix (GRM). Missing genotypes are imputed on the fly by drawing according to the minor allele frequencies.
# GRM using SNPs with maf > 0.01 (default) (~10 mins on my laptop) srand(123) @time Φgrm = grm(cg10k; method = :GRM)
396.943890 seconds (8.43 G allocations: 127.378 GiB, 4.38% gc time)
6670×6670 Array{Float64,2}:
0.503024 0.00335505 -0.000120075 … -5.45185e-5 -0.00278072
0.00335505 0.498958 -0.00195952 0.000868471 0.0034285
-0.000120075 -0.00195952 0.493828 0.000174648 -0.000381467
0.000923828 -0.00329169 -0.00194166 -0.00223595 -0.00123508
-8.39649e-5 -0.00353358 0.0018709 0.00222858 -0.00171176
0.00204208 0.000572952 0.00254025 … 0.000861385 2.99785e-5
0.000569323 0.0024786 -0.00185743 0.00117649 -0.00118027
-0.000642144 0.00317992 -0.00099777 0.00354182 -0.000260645
-0.00102913 -0.00123475 -0.00061138 0.00173885 0.00177727
-0.00139442 0.00208423 0.000124525 -0.00145156 -0.001011
-0.00204555 0.00011055 -0.000419398 … -0.000198235 -0.00110353
0.000947587 0.00167346 0.00184451 -0.000690143 -0.00304087
0.000322759 -0.000899805 0.00303981 0.000739331 -0.00118835
⋮ ⋱
0.00298012 0.00130003 0.000998861 4.18454e-6 0.00303991
-0.00207748 0.00274717 -0.00191741 -0.00107073 0.00368267
0.000545569 -0.00244439 -0.00299578 … -0.000669885 0.00221027
-0.00423186 -0.00208514 -0.00108833 -0.000622127 -0.000567483
-0.00325644 -0.000781353 0.0030423 0.000501423 -0.00010267
0.00041055 -0.00200772 0.00274867 -0.00624933 -0.00521365
0.00210519 0.000879889 -0.00107817 -0.000797878 -0.000557352
-0.00230058 -0.000119132 0.000116817 … 0.000867087 -0.00233512
-0.0020119 0.00230772 -0.00128837 0.00194798 -0.00048733
-0.000944942 -0.000928073 -0.000175096 0.00126911 -0.00303766
-5.45185e-5 0.000868471 0.000174648 0.500829 0.000469478
-0.00278072 0.0034285 -0.000381467 0.000469478 0.500627
Phenotypes
Read in the phenotype data and compute descriptive statistics.
# Pkg.add("DataFrames") using DataFrames cg10k_trait = readtable( "cg10k_traits.txt"; separator = ' ', names = [:FID; :IID; :Trait1; :Trait2; :Trait3; :Trait4; :Trait5; :Trait6; :Trait7; :Trait8; :Trait9; :Trait10; :Trait11; :Trait12; :Trait13], eltypes = [String; String; Float64; Float64; Float64; Float64; Float64; Float64; Float64; Float64; Float64; Float64; Float64; Float64; Float64] ) # do not display FID and IID for privacy cg10k_trait[:, 3:end]
| Trait1 | Trait2 | Trait3 | Trait4 | Trait5 | Trait6 | Trait7 | Trait8 | Trait9 | Trait10 | Trait11 | Trait12 | Trait13 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | -1.81573145026234 | -0.94615046147283 | 1.11363077580442 | -2.09867121119159 | 0.744416614111748 | 0.00139171884080131 | 0.934732480409667 | -1.22677315418103 | 1.1160784277875 | -0.4436280335029 | 0.824465656443384 | -1.02852542216546 | -0.394049201727681 |
| 2 | -1.24440094378729 | 0.109659992547179 | 0.467119394241789 | -1.62131304097589 | 1.0566758355683 | 0.978946979419181 | 1.00014633946047 | 0.32487427140228 | 1.16232175219696 | 2.6922706948705 | 3.08263672461047 | 1.09064954786013 | 0.0256616415357438 |
| 3 | 1.45566914502305 | 1.53866932923243 | 1.09402959376555 | 0.586655272226893 | -0.32796454430367 | -0.30337709778827 | -0.0334354881314741 | -0.464463064285437 | -0.3319396273436 | -0.486839089635991 | -1.10648681564373 | -1.42015780427231 | -0.687463456644413 |
| 4 | -0.768809276698548 | 0.513490885514249 | 0.244263028382142 | -1.31740254475691 | 1.19393774326845 | 1.17344127734288 | 1.08737426675232 | 0.536022583732261 | 0.802759240762068 | 0.234159411749815 | 0.394174866891074 | -0.767365892476029 | 0.0635385761884935 |
| 5 | -0.264415132547719 | -0.348240421825694 | -0.0239065083413606 | 0.00473915802244948 | 1.25619191712193 | 1.2038883667631 | 1.29800739042627 | 0.310113660247311 | 0.626159861059352 | 0.899289129831224 | 0.54996783350812 | 0.540687809542048 | 0.179675416046033 |
| 6 | -1.37617270917293 | -1.47191967744564 | 0.291179894254146 | -0.803110740704731 | -0.264239977442213 | -0.260573027836772 | -0.165372266287781 | -0.219257294118362 | 1.04702422290318 | -0.0985815534616482 | 0.947393438068448 | 0.594014812031438 | 0.245407436348479 |
| 7 | 0.1009416296374 | -0.191615722103455 | -0.567421321596677 | 0.378571487240382 | -0.246656179817904 | -0.608810750053858 | 0.189081058215596 | -1.27077787326519 | -0.452476199143965 | 0.702562877297724 | 0.332636218957179 | 0.0026916503626181 | 0.317117176705358 |
| 8 | -0.319818276367464 | 1.35774480657283 | 0.818689545938528 | -1.15565531644352 | 0.63448368102259 | 0.291461908634679 | 0.933323714954726 | -0.741083289682492 | 0.647477683507572 | -0.970877627077966 | 0.220861165411304 | 0.852512250237764 | -0.225904624283945 |
| 9 | -0.288334173342032 | 0.566082538090752 | 0.254958336116175 | -0.652578302869714 | 0.668921559277347 | 0.978309199170558 | 0.122862966041938 | 1.4790926378214 | 0.0672132424173449 | 0.0795903917527827 | 0.167532455243232 | 0.246915579442139 | 0.539932616458363 |
| 10 | -1.15759732583991 | -0.781198583545165 | -0.595807759833517 | -1.00554980260402 | 0.789828885933321 | 0.571058413379044 | 0.951304176233755 | -0.295962982984816 | 0.99042002479707 | 0.561309366988983 | 0.733100030623233 | -1.73467772245684 | -1.35278484330654 |
| 11 | 0.740569150459031 | 1.40873846755415 | 0.734689999440088 | 0.0208322841295094 | -0.337440968561619 | -0.458304040611395 | -0.142582512772326 | -0.580392297464107 | -0.684684998101516 | -0.00785381461893456 | -0.712244337518008 | -0.313345561230878 | -0.345419463162219 |
| 12 | -0.675892486454995 | 0.279892613829682 | 0.267915996308248 | -1.04103665392985 | 0.910741715645888 | 0.866027618513171 | 1.07414431702005 | 0.0381751003538302 | 0.766355377018601 | -0.340118016143495 | -0.809013958505059 | 0.548521663785885 | -0.0201828675962336 |
| 13 | -0.795410435603455 | -0.699989939762738 | 0.3991295030063 | -0.510476261900736 | 1.51552245416844 | 1.28743032939467 | 1.53772393250903 | 0.133989160117702 | 1.02025736886037 | 0.499018733899186 | -0.36948273277931 | -1.10153460436318 | -0.598132438886619 |
| 14 | -0.193483122930324 | -0.286021160323518 | -0.691494225262995 | 0.0131581678700699 | 1.52337470686782 | 1.4010638072262 | 1.53114620451896 | 0.333066483478075 | 1.04372480381099 | 0.163206783570466 | -0.422883765001728 | -0.383527976713573 | -0.489221907788158 |
| 15 | 0.151246203379718 | 2.09185108993614 | 2.03800472474384 | -1.12474717143531 | 1.66557024390713 | 1.62535675109576 | 1.58751070483655 | 0.635852186043776 | 0.842577784605979 | 0.450761870778952 | -1.39479033623028 | -0.560984107567768 | 0.289349776549287 |
| 16 | -0.464608740812712 | 0.36127694772303 | 1.2327673928287 | -0.826033731086383 | 1.43475224709983 | 1.74451823818846 | 0.211096887484638 | 2.64816425140548 | 1.02511433146096 | 0.11975731603184 | 0.0596832073448267 | -0.631231612661616 | -0.207878671782927 |
| 17 | -0.732977488012215 | -0.526223425889779 | 0.61657871336593 | -0.55447974332593 | 0.947484859025104 | 0.936833214138173 | 0.972516806335524 | 0.290251013865227 | 1.01285359725723 | 0.516207422283291 | -0.0300689171988194 | 0.8787322524583 | 0.450254629309513 |
| 18 | -0.167326459622119 | 0.175327165487237 | 0.287467725892572 | -0.402652532084246 | 0.551181509418056 | 0.522204743290975 | 0.436837660094653 | 0.299564933845579 | 0.583109520896067 | -0.704415820005353 | -0.730810367994577 | -1.95140580379896 | -0.933504665700164 |
| 19 | 1.41159485787418 | 1.78722407901017 | 0.84397639585364 | 0.481278083772991 | -0.0887673728508268 | -0.49957757426858 | 0.304195897924847 | -1.23884208383369 | -0.153475724036624 | -0.870486102788329 | 0.0955473331150403 | -0.983708050882817 | -0.3563445644514 |
| 20 | -1.42997091652825 | -0.490147045034213 | 0.272730237607695 | -1.61029992954153 | 0.990787817197748 | 0.711687532608184 | 1.1885836012715 | -0.371229188075638 | 1.24703459239952 | -0.0389162332271516 | 0.883495749072872 | 2.58988026321017 | 3.33539552370368 |
| 21 | -0.147247288176765 | 0.12328430415652 | 0.617549051912237 | -0.18713077178262 | 0.256438107586694 | 0.17794983735083 | 0.412611806463263 | -0.244809124559737 | 0.0947624806136492 | 0.723017223849532 | -0.683948354633436 | 0.0873751276309269 | -0.262209652750371 |
| 22 | -0.187112676773894 | -0.270777264595619 | -1.01556818551606 | 0.0602850568600233 | 0.272419757757978 | 0.869133161879197 | -0.657519461414234 | 2.32388522018189 | -0.999936011525034 | 1.44671844178306 | 0.971157886040772 | -0.358747904241515 | -0.439657942096136 |
| 23 | -1.82434047163768 | -0.933480446068067 | 1.29474003766977 | -1.94545221151036 | 0.33584651189654 | 0.359201654302844 | 0.513652924365886 | -0.073197696696958 | 1.57139042812005 | 1.53329371326728 | 1.82076821859528 | 2.22740301867829 | 1.50063347195857 |
| 24 | -2.29344084351335 | -2.49161842344418 | 0.40383988742336 | -2.36488074752948 | 1.4105254831956 | 1.42244117147792 | 1.17024166272172 | 0.84476650176855 | 1.79026875432495 | 0.648181858970515 | -0.0857231057403538 | -1.02789535292617 | 0.491288088952859 |
| 25 | -0.434135932888305 | 0.740881989034652 | 0.699576357578518 | -1.02405543187775 | 0.759529223983713 | 0.956656110895288 | 0.633299568656589 | 0.770733932268516 | 0.824988511714526 | 1.84287437634769 | 1.91045942063443 | -0.502317207869366 | 0.132670133448219 |
| 26 | -2.1920969546557 | -2.49465664272271 | 0.354854763893431 | -1.93155848635714 | 0.941979400289938 | 0.978917101414106 | 0.894860097289736 | 0.463239402831873 | 1.12537133317163 | 1.70528446191955 | 0.717792714479123 | 0.645888049108261 | 0.783968250169388 |
| 27 | -1.46602269088422 | -1.24921677101897 | 0.307977693653039 | -1.55097364660989 | 0.618908494474798 | 0.662508171662042 | 0.475957173906078 | 0.484718674597707 | 0.401564892028249 | 0.55987973254026 | -0.376938143754217 | -0.933982629228218 | 0.390013151672955 |
| 28 | -1.83317744236881 | -1.53268787828701 | 2.55674262685865 | -1.51827745783835 | 0.789409601746455 | 0.908747799728588 | 0.649971922941479 | 0.668373649931667 | 1.20058303519903 | 0.277963256075637 | 1.2504953198275 | 3.31370445071638 | 2.22035828885342 |
| 29 | -0.784546628243178 | 0.276582579543931 | 3.01104958800057 | -1.11978843206758 | 0.920823858422707 | 0.750217689886151 | 1.26153730009639 | -0.403363882922417 | 0.400667296857811 | -0.217597941303479 | -0.724669537565068 | -0.391945338467193 | -0.650023936358253 |
| 30 | 0.464455916345135 | 1.3326356122229 | -1.23059563374303 | -0.357975958937414 | 1.18249746977104 | 1.54315938069757 | -0.60339041154062 | 3.38308845958422 | 0.823740765148641 | -0.129951318508883 | -0.657979878422938 | -0.499534924074273 | -0.414476569095651 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
describe(cg10k_trait[:, 3:end])
Trait1 Summary Stats: Mean: 0.002211 Minimum: -3.204128 1st Quartile: -0.645771 Median: 0.125010 3rd Quartile: 0.723315 Maximum: 3.479398 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait2 Summary Stats: Mean: 0.001353 Minimum: -3.511659 1st Quartile: -0.642621 Median: 0.033517 3rd Quartile: 0.657467 Maximum: 4.913423 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait3 Summary Stats: Mean: -0.001296 Minimum: -3.938436 1st Quartile: -0.640907 Median: -0.000782 3rd Quartile: 0.637108 Maximum: 7.916299 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait4 Summary Stats: Mean: 0.002309 Minimum: -3.608403 1st Quartile: -0.546086 Median: 0.228165 3rd Quartile: 0.715291 Maximum: 3.127688 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait5 Summary Stats: Mean: -0.001790 Minimum: -4.148749 1st Quartile: -0.690765 Median: 0.031034 3rd Quartile: 0.734916 Maximum: 2.717184 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait6 Summary Stats: Mean: -0.001196 Minimum: -3.824792 1st Quartile: -0.662796 Median: 0.036242 3rd Quartile: 0.741176 Maximum: 2.589728 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait7 Summary Stats: Mean: -0.001989 Minimum: -4.272455 1st Quartile: -0.638923 Median: 0.069801 3rd Quartile: 0.710423 Maximum: 2.653779 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait8 Summary Stats: Mean: 0.000614 Minimum: -5.625488 1st Quartile: -0.601575 Median: -0.038630 3rd Quartile: 0.527342 Maximum: 5.805702 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait9 Summary Stats: Mean: -0.001810 Minimum: -5.381968 1st Quartile: -0.601429 Median: 0.106571 3rd Quartile: 0.698567 Maximum: 2.571936 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait10 Summary Stats: Mean: -0.000437 Minimum: -3.548506 1st Quartile: -0.633641 Median: -0.096651 3rd Quartile: 0.498610 Maximum: 6.537820 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait11 Summary Stats: Mean: -0.000616 Minimum: -3.264910 1st Quartile: -0.673685 Median: -0.068044 3rd Quartile: 0.655486 Maximum: 4.262410 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait12 Summary Stats: Mean: -0.000589 Minimum: -8.851909 1st Quartile: -0.539686 Median: -0.141099 3rd Quartile: 0.350779 Maximum: 13.211402 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000 Trait13 Summary Stats: Mean: -0.000151 Minimum: -5.592104 1st Quartile: -0.492289 Median: -0.141022 3rd Quartile: 0.324804 Maximum: 24.174436 Length: 6670 Type: Float64 Number Missing: 0 % Missing: 0.000000
Y = convert(Matrix{Float64}, cg10k_trait[:, 3:15]) histogram(Y, layout = 13)
Pre-processing data for heritability analysis
To prepare variance component model fitting, we form an instance of VarianceComponentVariate. The two variance components are $(2\Phi, I)$.
using VarianceComponentModels # form data as VarianceComponentVariate cg10kdata = VarianceComponentVariate(Y, (2Φgrm, eye(size(Y, 1)))) fieldnames(cg10kdata)
3-element Array{Symbol,1}:
:Y
:X
:V
cg10kdata
VarianceComponentModels.VarianceComponentVariate{Float64,2,Array{Float64,2},Array{Float64,2},Array{Float64,2}}([-1.81573 -0.94615 … -1.02853 -0.394049; -1.2444 0.10966 … 1.09065 0.0256616; … ; 0.886626 0.487408 … -0.636874 -0.439825; -1.24394 0.213697 … 0.299931 0.392809], Array{Float64}(6670,0), ([1.00605 0.00671009 … -0.000109037 -0.00556144; 0.00671009 0.997916 … 0.00173694 0.00685701; … ; -0.000109037 0.00173694 … 1.00166 0.000938955; -0.00556144 0.00685701 … 0.000938955 1.00125], [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0]))
Before fitting the variance component model, we pre-compute the eigen-decomposition of $2\Phi_{\text{GRM}}$, the rotated responses, and the constant part in log-likelihood, and store them as a TwoVarCompVariateRotate instance, which is re-used in various variane component estimation procedures.
# pre-compute eigen-decomposition (~50 secs on my laptop) @time cg10kdata_rotated = TwoVarCompVariateRotate(cg10kdata) fieldnames(cg10kdata_rotated)
48.837361 seconds (39 allocations: 1021.427 MiB, 0.57% gc time)
5-element Array{Symbol,1}:
:Yrot
:Xrot
:eigval
:eigvec
:logdetV2
Save intermediate results
We don't want to re-compute SnpArray and empirical kinship matrices again and again for heritibility analysis.
# # Pkg.add("JLD") # using JLD # @save "cg10k.jld" # whos()
To load workspace
using SnpArrays, JLD, DataFrames, VarianceComponentModels, Plots pyplot() @load "cg10k.jld" whos()
Base Module
BinDeps 41348 KB Module
Blosc 41202 KB Module
ColorTypes 41457 KB Module
Colors 41480 KB Module
Compat 41196 KB Module
Conda 41205 KB Module
Core Module
DataArrays 41456 KB Module
DataFrames 41684 KB Module
DataStructures 41356 KB Module
FileIO 41310 KB Module
FixedPointNumbers 41695 KB Module
GZip 41181 KB Module
HDF5 41403 KB Module
IJulia 4185781 KB Module
Ipopt 41172 KB Module
JLD 41376 KB Module
JSON 41245 KB Module
LaTeXStrings 4058 bytes Module
LegacyStrings 41212 KB Module
LinearMaps 22 KB Module
MacroTools 41606 KB Module
Main Module
MathProgBase 41353 KB Module
MbedTLS 41269 KB Module
Measures 41175 KB Module
NaNMath 41200 KB Module
PlotThemes 41167 KB Module
PlotUtils 41332 KB Module
Plots 42960 KB Module
PyCall 41711 KB Module
PyPlot 41771 KB Module
RecipesBase 41283 KB Module
Reexport 41160 KB Module
Requires 41172 KB Module
SHA 62 KB Module
Showoff 41163 KB Module
SnpArrays 41218 KB Module
SortingAlgorithms 41178 KB Module
SpecialFunctions 41252 KB Module
StaticArrays 41744 KB Module
StatsBase 41810 KB Module
URIParser 41171 KB Module
VarianceComponentModels 41278 KB Module
Y 677 KB 6670×13 Array{Float64,2}
ZMQ 41223 KB Module
_ 77 KB 630860-element BitArray{1}
cg10k 1027303 KB 6670×630860 SnpArrays.SnpArray{2}
cg10k_trait 978 KB 6670×15 DataFrames.DataFrame
cg10kdata 695816 KB VarianceComponentModels.VarianceCo…
cg10kdata_rotated 348299 KB VarianceComponentModels.TwoVarComp…
h 24 bytes 3-element Array{Float64,1}
hST 104 bytes 13-element Array{Float64,1}
hST_se 104 bytes 13-element Array{Float64,1}
hse 24 bytes 3-element Array{Float64,1}
maf 4928 KB 630860-element Array{Float64,1}
missings_by_snp 4928 KB 630860-element Array{Int64,1}
people 8 bytes Int64
snps 8 bytes Int64
trait57_data 347778 KB VarianceComponentModels.TwoVarComp…
trait57_model 232 bytes VarianceComponentModels.VarianceCo…
traitall_model 2792 bytes VarianceComponentModels.VarianceCo…
traitidx 16 bytes 3-element UnitRange{Int64}
Σa 3848 bytes 13×13 Array{Array{Float64,2},2}
Σcov 2592 bytes 18×18 Array{Float64,2}
Σe 3848 bytes 13×13 Array{Array{Float64,2},2}
Φgrm 347569 KB 6670×6670 Array{Float64,2}
σ2a 104 bytes 13-element Array{Float64,1}
σ2e 104 bytes 13-element Array{Float64,1}
Heritability of single traits
We use Fisher scoring algorithm to fit variance component model for each single trait.
# heritability from single trait analysis hST = zeros(13) # standard errors of estimated heritability hST_se = zeros(13) # additive genetic effects σ2a = zeros(13) # enviromental effects σ2e = zeros(13) tic() for trait in 1:13 println(names(cg10k_trait)[trait + 2]) # form data set for trait j traitj_data = TwoVarCompVariateRotate(cg10kdata_rotated.Yrot[:, trait], cg10kdata_rotated.Xrot, cg10kdata_rotated.eigval, cg10kdata_rotated.eigvec, cg10kdata_rotated.logdetV2) # initialize model parameters traitj_model = VarianceComponentModel(traitj_data) # estimate variance components _, _, _, Σcov, _, _ = mle_fs!(traitj_model, traitj_data; solver=:Ipopt, verbose=false) σ2a[trait] = traitj_model.Σ[1][1] σ2e[trait] = traitj_model.Σ[2][1] @show σ2a[trait], σ2e[trait] h, hse = heritability(traitj_model.Σ, Σcov) hST[trait] = h[1] hST_se[trait] = hse[1] end toc()
Trait1 (σ2a[trait], σ2e[trait]) = (0.25978160614793233, 0.7369535197912689) Trait2 (σ2a[trait], σ2e[trait]) = (0.18647130348299173, 0.8129591079735827) Trait3 (σ2a[trait], σ2e[trait]) = (0.3188368159422607, 0.6798809726936244) Trait4 (σ2a[trait], σ2e[trait]) = (0.2651357653143703, 0.7308007669086968) Trait5 (σ2a[trait], σ2e[trait]) = (0.28083388108246, 0.7172036435586534) Trait6 (σ2a[trait], σ2e[trait]) = (0.2824159905728832, 0.7170988773569172) Trait7 (σ2a[trait], σ2e[trait]) = (0.2155274336968625, 0.7815346282986375) Trait8 (σ2a[trait], σ2e[trait]) = (0.194687807263945, 0.8049690651320599) Trait9 (σ2a[trait], σ2e[trait]) = (0.24706855916591713, 0.7512942998567308) Trait10 (σ2a[trait], σ2e[trait]) = (0.098712236297271, 0.9011756660217387) Trait11 (σ2a[trait], σ2e[trait]) = (0.1664264642608195, 0.8322427413046204) Trait12 (σ2a[trait], σ2e[trait]) = (0.0834296761650666, 0.9153609794266608) Trait13 (σ2a[trait], σ2e[trait]) = (0.05893968504298988, 0.940270012443928) elapsed time: 0.160999612 seconds 0.160999612
# heritability and standard errors [hST'; hST_se']
2×13 Array{Float64,2}:
0.260633 0.186578 0.319246 … 0.166648 0.0835307 0.0589863
0.0799732 0.0869002 0.0741007 0.08862 0.0944407 0.0953238
Pairwise traits
Joint analysis of multiple traits is subject to intensive research recently. Following code snippet does joint analysis of all pairs of traits, a total of 78 bivariate variane component models.
# additive genetic effects (2x2 psd matrices) from bavariate trait analysis; Σa = Array{Matrix{Float64}}(13, 13) # environmental effects (2x2 psd matrices) from bavariate trait analysis; Σe = Array{Matrix{Float64}}(13, 13) tic() for i in 1:13 for j in (i+1):13 println(names(cg10k_trait)[i + 2], names(cg10k_trait)[j + 2]) # form data set for (trait1, trait2) traitij_data = TwoVarCompVariateRotate(cg10kdata_rotated.Yrot[:, [i;j]], cg10kdata_rotated.Xrot, cg10kdata_rotated.eigval, cg10kdata_rotated.eigvec, cg10kdata_rotated.logdetV2) # initialize model parameters traitij_model = VarianceComponentModel(traitij_data) # estimate variance components mle_fs!(traitij_model, traitij_data; solver=:Ipopt, verbose=false) Σa[i, j] = traitij_model.Σ[1] Σe[i, j] = traitij_model.Σ[2] @show Σa[i, j], Σe[i, j] end end toc()
Trait1Trait2 (Σa[i, j], Σe[i, j]) = ([0.258822 0.174358; 0.174358 0.185108], [0.737892 0.585751; 0.585751 0.814301]) Trait1Trait3 (Σa[i, j], Σe[i, j]) = ([0.260236 -0.0144726; -0.0144726 0.319245], [0.736512 -0.11979; -0.11979 0.679488]) Trait1Trait4 (Σa[i, j], Σe[i, j]) = ([0.259615 0.222203; 0.222203 0.265149], [0.737116 0.599854; 0.599854 0.730788]) Trait1Trait5 (Σa[i, j], Σe[i, j]) = ([0.259574 -0.146827; -0.146827 0.28153], [0.737153 -0.254777; -0.254777 0.71653]) Trait1Trait6 (Σa[i, j], Σe[i, j]) = ([0.259476 -0.129115; -0.129115 0.282688], [0.73725 -0.23161; -0.23161 0.716837]) Trait1Trait7 (Σa[i, j], Σe[i, j]) = ([0.259115 -0.140455; -0.140455 0.215297], [0.737606 -0.197616; -0.197616 0.781774]) Trait1Trait8 (Σa[i, j], Σe[i, j]) = ([0.259778 -0.0327756; -0.0327756 0.194698], [0.736957 -0.127026; -0.127026 0.804959]) Trait1Trait9 (Σa[i, j], Σe[i, j]) = ([0.261858 -0.204589; -0.204589 0.246027], [0.734961 -0.307734; -0.307734 0.75232]) Trait1Trait10 (Σa[i, j], Σe[i, j]) = ([0.259649 -0.0994858; -0.0994858 0.0956585], [0.737083 -0.303942; -0.303942 0.904218]) Trait1Trait11 (Σa[i, j], Σe[i, j]) = ([0.25947 -0.138603; -0.138603 0.164709], [0.737257 -0.359557; -0.359557 0.83395]) Trait1Trait12 (Σa[i, j], Σe[i, j]) = ([0.261779 -0.145414; -0.145414 0.0807748], [0.735076 -0.041823; -0.041823 0.9181]) Trait1Trait13 (Σa[i, j], Σe[i, j]) = ([0.261125 -0.108774; -0.108774 0.0538214], [0.735674 -0.114123; -0.114123 0.945416]) Trait2Trait3 (Σa[i, j], Σe[i, j]) = ([0.186541 0.144056; 0.144056 0.320627], [0.812888 0.0995944; 0.0995944 0.678167]) Trait2Trait4 (Σa[i, j], Σe[i, j]) = ([0.186131 0.0746032; 0.0746032 0.265122], [0.813293 0.221109; 0.221109 0.730814]) Trait2Trait5 (Σa[i, j], Σe[i, j]) = ([0.186442 -0.0118093; -0.0118093 0.280842], [0.812987 -0.0365191; -0.0365191 0.717195]) Trait2Trait6 (Σa[i, j], Σe[i, j]) = ([0.18649 -0.00366533; -0.00366533 0.282471], [0.812941 -0.0206271; -0.0206271 0.717046]) Trait2Trait7 (Σa[i, j], Σe[i, j]) = ([0.186104 -0.030665; -0.030665 0.215304], [0.81332 -0.000667009; -0.000667009 0.781755]) Trait2Trait8 (Σa[i, j], Σe[i, j]) = ([0.187023 0.0331783; 0.0331783 0.195259], [0.812421 -0.0326343; -0.0326343 0.804415]) Trait2Trait9 (Σa[i, j], Σe[i, j]) = ([0.185032 -0.085334; -0.085334 0.245909], [0.814386 -0.0809638; -0.0809638 0.752433]) Trait2Trait10 (Σa[i, j], Σe[i, j]) = ([0.186587 -0.123303; -0.123303 0.0987387], [0.812872 -0.273083; -0.273083 0.901229]) Trait2Trait11 (Σa[i, j], Σe[i, j]) = ([0.185484 -0.117256; -0.117256 0.167776], [0.81393 -0.296772; -0.296772 0.830934]) Trait2Trait12 (Σa[i, j], Σe[i, j]) = ([0.185907 -0.0909104; -0.0909104 0.0827171], [0.813555 0.0457924; 0.0457924 0.916135]) Trait2Trait13 (Σa[i, j], Σe[i, j]) = ([0.185979 -0.0720811; -0.0720811 0.0568238], [0.8135 0.0751703; 0.0751703 0.942424]) Trait3Trait4 (Σa[i, j], Σe[i, j]) = ([0.3188 -0.154562; -0.154562 0.264323], [0.679917 -0.303223; -0.303223 0.731591]) Trait3Trait5 (Σa[i, j], Σe[i, j]) = ([0.319216 0.183527; 0.183527 0.282063], [0.679514 0.33724; 0.33724 0.716008]) Trait3Trait6 (Σa[i, j], Σe[i, j]) = ([0.319776 0.165672; 0.165672 0.284448], [0.678972 0.298667; 0.298667 0.715124]) Trait3Trait7 (Σa[i, j], Σe[i, j]) = ([0.318838 0.166283; 0.166283 0.215261], [0.67988 0.347706; 0.347706 0.781796]) Trait3Trait8 (Σa[i, j], Σe[i, j]) = ([0.320718 0.0566397; 0.0566397 0.197764], [0.678063 0.0451569; 0.0451569 0.801956]) Trait3Trait9 (Σa[i, j], Σe[i, j]) = ([0.319001 0.137699; 0.137699 0.246142], [0.679722 0.266704; 0.266704 0.752197]) Trait3Trait10 (Σa[i, j], Σe[i, j]) = ([0.31908 -0.076513; -0.076513 0.0996001], [0.679646 -0.142905; -0.142905 0.900298]) Trait3Trait11 (Σa[i, j], Σe[i, j]) = ([0.318094 -0.0177494; -0.0177494 0.16629], [0.6806 -0.1144; -0.1144 0.832376]) Trait3Trait12 (Σa[i, j], Σe[i, j]) = ([0.321164 0.0843842; 0.0843842 0.0874609], [0.677639 0.0341558; 0.0341558 0.911368]) Trait3Trait13 (Σa[i, j], Σe[i, j]) = ([0.323273 0.109443; 0.109443 0.0634295], [0.675635 -0.0060525; -0.0060525 0.935819]) Trait4Trait5 (Σa[i, j], Σe[i, j]) = ([0.26525 -0.215125; -0.215125 0.282572], [0.73068 -0.377406; -0.377406 0.715518]) Trait4Trait6 (Σa[i, j], Σe[i, j]) = ([0.265715 -0.199714; -0.199714 0.283942], [0.730231 -0.347732; -0.347732 0.715619]) Trait4Trait7 (Σa[i, j], Σe[i, j]) = ([0.26407 -0.18238; -0.18238 0.214324], [0.731843 -0.32655; -0.32655 0.782733]) Trait4Trait8 (Σa[i, j], Σe[i, j]) = ([0.266229 -0.0965381; -0.0965381 0.196655], [0.729739 -0.151461; -0.151461 0.803044]) Trait4Trait9 (Σa[i, j], Σe[i, j]) = ([0.269627 -0.226931; -0.226931 0.247265], [0.726443 -0.416085; -0.416085 0.751086]) Trait4Trait10 (Σa[i, j], Σe[i, j]) = ([0.265098 -0.0352926; -0.0352926 0.0981462], [0.730847 -0.226248; -0.226248 0.901736]) Trait4Trait11 (Σa[i, j], Σe[i, j]) = ([0.265178 -0.0970634; -0.0970634 0.164885], [0.73076 -0.272291; -0.272291 0.833762]) Trait4Trait12 (Σa[i, j], Σe[i, j]) = ([0.267732 -0.140985; -0.140985 0.081029], [0.728323 -0.0834791; -0.0834791 0.917815]) Trait4Trait13 (Σa[i, j], Σe[i, j]) = ([0.265695 -0.0970238; -0.0970238 0.0564809], [0.730259 -0.226115; -0.226115 0.942736]) Trait5Trait6 (Σa[i, j], Σe[i, j]) = ([0.281198 0.280259; 0.280259 0.281764], [0.716855 0.661013; 0.661013 0.717735]) Trait5Trait7 (Σa[i, j], Σe[i, j]) = ([0.280442 0.231918; 0.231918 0.211837], [0.717597 0.674491; 0.674491 0.785172]) Trait5Trait8 (Σa[i, j], Σe[i, j]) = ([0.280958 0.163168; 0.163168 0.193315], [0.717089 0.221817; 0.221817 0.806314]) Trait5Trait9 (Σa[i, j], Σe[i, j]) = ([0.283544 0.243884; 0.243884 0.240564], [0.714585 0.509072; 0.509072 0.757631]) Trait5Trait10 (Σa[i, j], Σe[i, j]) = ([0.281378 -0.0454427; -0.0454427 0.100081], [0.716678 -0.0579778; -0.0579778 0.899822]) Trait5Trait11 (Σa[i, j], Σe[i, j]) = ([0.280066 0.0195669; 0.0195669 0.165607], [0.71795 -0.0345589; -0.0345589 0.833047]) Trait5Trait12 (Σa[i, j], Σe[i, j]) = ([0.28101 0.0592641; 0.0592641 0.0831831], [0.717036 0.0552788; 0.0552788 0.915608]) Trait5Trait13 (Σa[i, j], Σe[i, j]) = ([0.281854 0.0680641; 0.0680641 0.0591899], [0.716223 0.0551992; 0.0551992 0.940027]) Trait6Trait7 (Σa[i, j], Σe[i, j]) = ([0.282435 0.220236; 0.220236 0.213997], [0.71708 0.581507; 0.581507 0.783041]) Trait6Trait8 (Σa[i, j], Σe[i, j]) = ([0.282435 0.18375; 0.18375 0.192999], [0.717081 0.436932; 0.436932 0.80663]) Trait6Trait9 (Σa[i, j], Σe[i, j]) = ([0.284516 0.233768; 0.233768 0.242478], [0.715071 0.477502; 0.477502 0.755765]) Trait6Trait10 (Σa[i, j], Σe[i, j]) = ([0.283087 -0.0427658; -0.0427658 0.100634], [0.716449 -0.0599491; -0.0599491 0.899275]) Trait6Trait11 (Σa[i, j], Σe[i, j]) = ([0.281046 0.0272144; 0.0272144 0.165044], [0.71843 -0.0516242; -0.0516242 0.833601]) Trait6Trait12 (Σa[i, j], Σe[i, j]) = ([0.28256 0.0548537; 0.0548537 0.083133], [0.716961 0.0502064; 0.0502064 0.915658]) Trait6Trait13 (Σa[i, j], Σe[i, j]) = ([0.283231 0.0585667; 0.0585667 0.0592752], [0.716314 0.055827; 0.055827 0.939942]) Trait7Trait8 (Σa[i, j], Σe[i, j]) = ([0.213998 0.0875641; 0.0875641 0.192993], [0.78304 -0.055939; -0.055939 0.806635]) Trait7Trait9 (Σa[i, j], Σe[i, j]) = ([0.219039 0.216925; 0.216925 0.243338], [0.778156 0.463024; 0.463024 0.754935]) Trait7Trait10 (Σa[i, j], Σe[i, j]) = ([0.216296 -0.0412106; -0.0412106 0.100663], [0.780785 -0.0868086; -0.0868086 0.899246]) Trait7Trait11 (Σa[i, j], Σe[i, j]) = ([0.2142 0.0204227; 0.0204227 0.165077], [0.782833 -0.0478727; -0.0478727 0.833569]) Trait7Trait12 (Σa[i, j], Σe[i, j]) = ([0.215054 0.0738562; 0.0738562 0.0814228], [0.782012 0.0366272; 0.0366272 0.917365]) Trait7Trait13 (Σa[i, j], Σe[i, j]) = ([0.216093 0.0728515; 0.0728515 0.0570272], [0.781006 0.0409945; 0.0409945 0.942189]) Trait8Trait9 (Σa[i, j], Σe[i, j]) = ([0.195154 0.111756; 0.111756 0.246453], [0.804528 0.185842; 0.185842 0.751896]) Trait8Trait10 (Σa[i, j], Σe[i, j]) = ([0.195015 -0.015506; -0.015506 0.0990776], [0.804651 0.0118538; 0.0118538 0.900815]) Trait8Trait11 (Σa[i, j], Σe[i, j]) = ([0.194421 0.0215044; 0.0215044 0.166226], [0.805231 -0.026247; -0.026247 0.83244]) Trait8Trait12 (Σa[i, j], Σe[i, j]) = ([0.194491 -0.00425152; -0.00425152 0.0832711], [0.805162 0.0349872; 0.0349872 0.915518]) Trait8Trait13 (Σa[i, j], Σe[i, j]) = ([0.19448 0.00235501; 0.00235501 0.0589351], [0.805173 0.0396048; 0.0396048 0.940275]) Trait9Trait10 (Σa[i, j], Σe[i, j]) = ([0.246455 -0.00257997; -0.00257997 0.0984563], [0.751895 0.0743439; 0.0743439 0.901429]) Trait9Trait11 (Σa[i, j], Σe[i, j]) = ([0.247001 0.0303415; 0.0303415 0.166421], [0.75136 0.153765; 0.153765 0.832248]) Trait9Trait12 (Σa[i, j], Σe[i, j]) = ([0.249421 0.0829968; 0.0829968 0.0890874], [0.749007 0.109331; 0.109331 0.909778]) Trait9Trait13 (Σa[i, j], Σe[i, j]) = ([0.24861 0.0916799; 0.0916799 0.0602352], [0.749811 0.100027; 0.100027 0.939032]) Trait10Trait11 (Σa[i, j], Σe[i, j]) = ([0.0914658 0.100613; 0.100613 0.166501], [0.908376 0.473847; 0.473847 0.83217]) Trait10Trait12 (Σa[i, j], Σe[i, j]) = ([0.0951392 0.0588424; 0.0588424 0.0796744], [0.904735 0.0828862; 0.0828862 0.919115]) Trait10Trait13 (Σa[i, j], Σe[i, j]) = ([0.0995192 -0.0257171; -0.0257171 0.0598595], [0.900397 0.163778; 0.163778 0.939368]) Trait11Trait12 (Σa[i, j], Σe[i, j]) = ([0.165386 0.0579914; 0.0579914 0.0796005], [0.83327 0.144637; 0.144637 0.919166]) Trait11Trait13 (Σa[i, j], Σe[i, j]) = ([0.166417 -0.000985185; -0.000985185 0.0595681], [0.832265 0.200012; 0.200012 0.939646]) Trait12Trait13 (Σa[i, j], Σe[i, j]) = ([0.085082 0.0696185; 0.0696185 0.0569655], [0.913729 0.572041; 0.572041 0.942247]) elapsed time: 3.587337102 seconds 3.587337102
3-trait analysis
Researchers want to jointly analyze traits 5-7. Our strategy is to try both Fisher scoring and MM algorithm with different starting point, and choose the best local optimum. We first form the data set and run Fisher scoring, which yields a final objective value -1.4700991+04.
traitidx = 5:7 # form data set trait57_data = TwoVarCompVariateRotate(cg10kdata_rotated.Yrot[:, traitidx], cg10kdata_rotated.Xrot, cg10kdata_rotated.eigval, cg10kdata_rotated.eigvec, cg10kdata_rotated.logdetV2) # initialize model parameters trait57_model = VarianceComponentModel(trait57_data) # estimate variance components @time mle_fs!(trait57_model, trait57_data; solver=:Ipopt, verbose=true) trait57_model
This is Ipopt version 3.12.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 78
Total number of variables............................: 12
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 3.0247565e+04 0.00e+00 1.00e+02 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
5 1.6835078e+04 0.00e+00 4.08e+02 -11.0 3.64e-01 - 1.00e+00 1.00e+00f 1 MaxS
10 1.4742941e+04 0.00e+00 1.10e+02 -11.0 2.35e-01 - 1.00e+00 1.00e+00f 1 MaxS
15 1.4701394e+04 0.00e+00 1.16e+01 -11.0 7.78e-02 -4.5 1.00e+00 1.00e+00f 1 MaxS
20 1.4701019e+04 0.00e+00 5.75e-01 -11.0 1.51e-04 -6.9 1.00e+00 1.00e+00f 1 MaxS
25 1.4701018e+04 0.00e+00 2.40e-02 -11.0 6.38e-06 -9.2 1.00e+00 1.00e+00f 1 MaxS
30 1.4701018e+04 0.00e+00 9.98e-04 -11.0 2.66e-07 -11.6 1.00e+00 1.00e+00f 1 MaxS
35 1.4701018e+04 0.00e+00 4.15e-05 -11.0 1.10e-08 -14.0 1.00e+00 1.00e+00h 1 MaxS
40 1.4701018e+04 0.00e+00 1.72e-06 -11.0 4.59e-10 -16.4 1.00e+00 1.00e+00f 1 MaxSA
45 1.4701018e+04 0.00e+00 7.17e-08 -11.0 1.91e-11 -18.8 1.00e+00 1.00e+00h 1 MaxSA
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
Number of Iterations....: 49
(scaled) (unscaled)
Objective...............: 4.4720359684330265e+02 1.4701017692082147e+04
Dual infeasibility......: 5.6081357364421780e-09 1.8435742302386474e-07
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 5.6081357364421780e-09 1.8435742302386474e-07
Number of objective function evaluations = 50
Number of objective gradient evaluations = 50
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 49
Total CPU secs in IPOPT (w/o function evaluations) = 0.014
Total CPU secs in NLP function evaluations = 0.076
EXIT: Optimal Solution Found.
0.097955 seconds (55.15 k allocations: 5.632 MiB)
VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,3), ([0.280777 0.279441 0.232208; 0.279441 0.28422 0.219831; 0.232208 0.219831 0.212832], [0.717266 0.66183 0.674206; 0.66183 0.715287 0.581891; 0.674206 0.581891 0.784183]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf)
We then run the MM algorithm, starting from the Fisher scoring answer. MM finds an improved solution with objective value 8.955397e+03.
# trait59_model contains the fitted model by Fisher scoring now @time mle_mm!(trait57_model, trait57_data; verbose=true) trait57_model
MM Algorithm
Iter Objective
-------- -------------
0 -1.470102e+04
1 -1.470102e+04
0.003006 seconds (21.01 k allocations: 1.551 MiB)
VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,3), ([0.280777 0.279441 0.232208; 0.279441 0.28422 0.219831; 0.232208 0.219831 0.212832], [0.717266 0.66183 0.674206; 0.66183 0.715287 0.581891; 0.674206 0.581891 0.784183]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf)
Do another run of MM algorithm from default starting point. It leads to a slightly better local optimum -1.470104e+04, slighly worse than the Fisher scoring result. Follow up anlaysis should use the Fisher scoring result.
# default starting point trait57_model = VarianceComponentModel(trait57_data) @time _, _, _, Σcov, = mle_mm!(trait57_model, trait57_data; verbose=true) trait57_model
MM Algorithm
Iter Objective
-------- -------------
0 -3.024757e+04
1 -2.040300e+04
2 -1.656070e+04
3 -1.528529e+04
4 -1.490986e+04
5 -1.480638e+04
6 -1.477811e+04
7 -1.476968e+04
8 -1.476639e+04
9 -1.476444e+04
10 -1.476286e+04
20 -1.475000e+04
30 -1.474011e+04
40 -1.473248e+04
50 -1.472658e+04
60 -1.472200e+04
70 -1.471840e+04
80 -1.471555e+04
90 -1.471328e+04
100 -1.471145e+04
110 -1.470997e+04
120 -1.470875e+04
130 -1.470775e+04
140 -1.470691e+04
150 -1.470621e+04
160 -1.470562e+04
170 -1.470511e+04
180 -1.470469e+04
190 -1.470432e+04
200 -1.470400e+04
210 -1.470372e+04
220 -1.470348e+04
230 -1.470326e+04
240 -1.470308e+04
250 -1.470291e+04
260 -1.470276e+04
270 -1.470263e+04
280 -1.470251e+04
290 -1.470241e+04
300 -1.470231e+04
310 -1.470223e+04
320 -1.470215e+04
330 -1.470208e+04
340 -1.470201e+04
350 -1.470195e+04
360 -1.470190e+04
370 -1.470185e+04
380 -1.470180e+04
390 -1.470176e+04
400 -1.470172e+04
410 -1.470168e+04
420 -1.470165e+04
430 -1.470162e+04
440 -1.470159e+04
450 -1.470156e+04
460 -1.470153e+04
470 -1.470151e+04
480 -1.470149e+04
490 -1.470147e+04
500 -1.470145e+04
510 -1.470143e+04
520 -1.470141e+04
530 -1.470139e+04
540 -1.470138e+04
550 -1.470136e+04
560 -1.470135e+04
570 -1.470133e+04
580 -1.470132e+04
590 -1.470131e+04
600 -1.470130e+04
610 -1.470129e+04
620 -1.470128e+04
630 -1.470127e+04
640 -1.470126e+04
650 -1.470125e+04
660 -1.470124e+04
670 -1.470123e+04
680 -1.470122e+04
690 -1.470122e+04
700 -1.470121e+04
710 -1.470120e+04
720 -1.470120e+04
730 -1.470119e+04
740 -1.470118e+04
750 -1.470118e+04
760 -1.470117e+04
770 -1.470117e+04
780 -1.470116e+04
790 -1.470116e+04
800 -1.470115e+04
810 -1.470115e+04
820 -1.470114e+04
830 -1.470114e+04
840 -1.470114e+04
850 -1.470113e+04
860 -1.470113e+04
870 -1.470112e+04
880 -1.470112e+04
890 -1.470112e+04
900 -1.470111e+04
910 -1.470111e+04
920 -1.470111e+04
930 -1.470111e+04
940 -1.470110e+04
950 -1.470110e+04
960 -1.470110e+04
970 -1.470109e+04
980 -1.470109e+04
990 -1.470109e+04
1000 -1.470109e+04
1010 -1.470109e+04
1020 -1.470108e+04
1030 -1.470108e+04
1040 -1.470108e+04
1050 -1.470108e+04
1060 -1.470108e+04
1070 -1.470107e+04
1080 -1.470107e+04
1090 -1.470107e+04
1100 -1.470107e+04
1110 -1.470107e+04
1120 -1.470107e+04
0.794377 seconds (168.12 k allocations: 15.640 MiB, 0.80% gc time)
VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,3), ([0.2808 0.279454 0.232256; 0.279454 0.284312 0.219977; 0.232256 0.219977 0.213052], [0.717243 0.661816 0.674158; 0.661816 0.715193 0.581746; 0.674158 0.581746 0.783965]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf)
Heritability from 3-variate estimate and their standard errors.
h, hse = heritability(trait57_model.Σ, Σcov) [h'; hse']
2×3 Array{Float64,2}:
0.281351 0.284453 0.213689
0.0778252 0.077378 0.084084
13-trait joint analysis
In some situations, such as studying the genetic covariance, we need to jointly analyze 13 traits. We first try the Fisher scoring algorithm.
# initialize model parameters traitall_model = VarianceComponentModel(cg10kdata_rotated) # estimate variance components using Fisher scoring algorithm @time mle_fs!(traitall_model, cg10kdata_rotated; solver=:Ipopt, verbose=true)
This is Ipopt version 3.12.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 16653
Total number of variables............................: 182
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.3113371e+05 0.00e+00 1.00e+02 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
5 8.2233766e+04 0.00e+00 6.03e+02 -11.0 2.32e+00 - 1.00e+00 1.00e+00f 1 MaxS
10 1.1960260e+05 0.00e+00 8.76e+02 -11.0 6.20e+01 -5.4 1.00e+00 1.00e+00h 1 MaxS
15 2.4416551e+05 0.00e+00 2.50e+02 -11.0 8.69e+02 -7.8 1.00e+00 1.00e+00f 1 MaxS
DomainError:
log will only return a complex result if called with a complex argument. Try log(complex(x)).
Stacktrace:
[1] nan_dom_err at ./math.jl:300 [inlined]
[2] log at ./math.jl:419 [inlined]
[3] logdet(::Array{Float64,2}) at ./linalg/generic.jl:1244
[4] VarianceComponentModels.TwoVarCompModelRotate(::VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}) at /Users/huazhou/.julia/v0.6/VarianceComponentModels/src/VarianceComponentModels.jl:127
[5] eval_f(::VarianceComponentModels.TwoVarCompOptProb{VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}},VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}},Array{Float64,2},Array{Float64,1},VarianceComponentModels.VarianceComponentAuxData{Array{Float64,2},Array{Float64,1}}}, ::Array{Float64,1}) at /Users/huazhou/.julia/v0.6/VarianceComponentModels/src/two_variance_component.jl:683
[6] (::Ipopt.#eval_f_cb#4{VarianceComponentModels.TwoVarCompOptProb{VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}},VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}},Array{Float64,2},Array{Float64,1},VarianceComponentModels.VarianceComponentAuxData{Array{Float64,2},Array{Float64,1}}}})(::Array{Float64,1}) at /Users/huazhou/.julia/v0.6/Ipopt/src/IpoptSolverInterface.jl:53
[7] eval_f_wrapper(::Int32, ::Ptr{Float64}, ::Int32, ::Ptr{Float64}, ::Ptr{Void}) at /Users/huazhou/.julia/v0.6/Ipopt/src/Ipopt.jl:89
[8] solveProblem(::Ipopt.IpoptProblem) at /Users/huazhou/.julia/v0.6/Ipopt/src/Ipopt.jl:304
[9] optimize!(::Ipopt.IpoptMathProgModel) at /Users/huazhou/.julia/v0.6/Ipopt/src/IpoptSolverInterface.jl:120
[10] #mle_fs!#29(::Int64, ::Symbol, ::Symbol, ::Bool, ::Function, ::VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}, ::VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}}) at /Users/huazhou/.julia/v0.6/VarianceComponentModels/src/two_variance_component.jl:893
[11] (::VarianceComponentModels.#kw##mle_fs!)(::Array{Any,1}, ::VarianceComponentModels.#mle_fs!, ::VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}, ::VarianceComponentModels.TwoVarCompVariateRotate{Float64,Array{Float64,2},Array{Float64,2}}) at ./<missing>:0
[12] include_string(::String, ::String) at ./loading.jl:515
From the output we can see the Fisher scoring algorithm ran into some numerical issues. Let's try the MM algorithm.
# reset model parameters traitall_model = VarianceComponentModel(cg10kdata_rotated) # estimate variance components using Fisher scoring algorithm @time mle_mm!(traitall_model, cg10kdata_rotated; verbose=true)
MM Algorithm
Iter Objective
-------- -------------
0 -1.311337e+05
1 -8.002108e+04
2 -5.806935e+04
3 -4.926111e+04
4 -4.611059e+04
5 -4.511606e+04
6 -4.482679e+04
7 -4.474294e+04
8 -4.471496e+04
9 -4.470174e+04
10 -4.469246e+04
20 -4.462243e+04
30 -4.456888e+04
40 -4.452774e+04
50 -4.449601e+04
60 -4.447134e+04
70 -4.445199e+04
80 -4.443665e+04
90 -4.442436e+04
100 -4.441442e+04
110 -4.440630e+04
120 -4.439961e+04
130 -4.439405e+04
140 -4.438938e+04
150 -4.438544e+04
160 -4.438210e+04
170 -4.437923e+04
180 -4.437676e+04
190 -4.437463e+04
200 -4.437277e+04
210 -4.437115e+04
220 -4.436972e+04
230 -4.436846e+04
240 -4.436735e+04
250 -4.436636e+04
260 -4.436548e+04
270 -4.436469e+04
280 -4.436399e+04
290 -4.436335e+04
300 -4.436278e+04
310 -4.436226e+04
320 -4.436179e+04
330 -4.436137e+04
340 -4.436098e+04
350 -4.436063e+04
360 -4.436030e+04
370 -4.436001e+04
380 -4.435974e+04
390 -4.435949e+04
400 -4.435926e+04
410 -4.435905e+04
420 -4.435886e+04
430 -4.435868e+04
440 -4.435851e+04
450 -4.435836e+04
460 -4.435822e+04
470 -4.435809e+04
480 -4.435797e+04
490 -4.435785e+04
500 -4.435775e+04
510 -4.435765e+04
520 -4.435756e+04
530 -4.435747e+04
540 -4.435739e+04
550 -4.435732e+04
560 -4.435725e+04
570 -4.435718e+04
580 -4.435712e+04
590 -4.435706e+04
600 -4.435701e+04
610 -4.435696e+04
620 -4.435691e+04
630 -4.435687e+04
640 -4.435683e+04
650 -4.435679e+04
660 -4.435675e+04
670 -4.435671e+04
680 -4.435668e+04
690 -4.435665e+04
700 -4.435662e+04
710 -4.435659e+04
720 -4.435657e+04
730 -4.435654e+04
740 -4.435652e+04
750 -4.435649e+04
760 -4.435647e+04
770 -4.435645e+04
780 -4.435643e+04
790 -4.435642e+04
800 -4.435640e+04
810 -4.435638e+04
820 -4.435637e+04
830 -4.435635e+04
840 -4.435634e+04
850 -4.435633e+04
860 -4.435631e+04
870 -4.435630e+04
880 -4.435629e+04
890 -4.435628e+04
900 -4.435627e+04
910 -4.435626e+04
920 -4.435625e+04
930 -4.435624e+04
940 -4.435623e+04
950 -4.435622e+04
960 -4.435621e+04
970 -4.435621e+04
980 -4.435620e+04
990 -4.435619e+04
1000 -4.435619e+04
1010 -4.435618e+04
1020 -4.435617e+04
1030 -4.435617e+04
1040 -4.435616e+04
1050 -4.435616e+04
1060 -4.435615e+04
1070 -4.435615e+04
1080 -4.435614e+04
1090 -4.435614e+04
3.551301 seconds (178.42 k allocations: 70.115 MiB, 0.42% gc time)
(-44356.138529861186, VarianceComponentModels.VarianceComponentModel{Float64,2,Array{Float64,2},Array{Float64,2}}(Array{Float64}(0,13), ([0.272384 0.190358 … -0.128222 -0.0980655; 0.190358 0.21692 … -0.0689912 -0.0444349; … ; -0.128222 -0.0689912 … 0.118227 0.0909188; -0.0980655 -0.0444349 … 0.0909188 0.107456], [0.724562 0.56992 … -0.0590518 -0.124939; 0.56992 0.782639 … 0.0238629 0.0475408; … ; -0.0590518 0.0238629 … 0.880671 0.550889; -0.124939 0.0475408 … 0.550889 0.891929]), Array{Float64}(0,0), Char[], Float64[], -Inf, Inf), ([0.0111619 0.0131088 … 0.0128956 0.0127641; 0.0131091 0.0151759 … 0.017162 0.0171466; … ; 0.0128956 0.017162 … 0.0173994 0.0182002; 0.0127643 0.0171461 … 0.0182003 0.0187848], [0.0112235 0.0133094 … 0.0130111 0.0127861; 0.01331 0.0158262 … 0.017867 0.017798; … ; 0.013011 0.0178666 … 0.0179487 0.0187579; 0.012786 0.0177975 … 0.0187578 0.0193328]), [0.000124587 7.24074e-5 … -3.35716e-7 -1.40982e-5; 7.24411e-5 0.000171849 … -2.05381e-5 -3.17975e-6; … ; -3.60221e-7 -2.05683e-5 … 0.000351859 -1.5168e-5; -1.40799e-5 -3.16738e-6 … -1.51641e-5 0.000373756], Array{Float64}(0,13), Array{Float64}(0,0))
It converges after ~1000 iterations.
Save analysis results
#using JLD #@save "copd.jld" #whos()