Single-step SNP-BLUP with on-the-fly imputed genotypes and residual polygenic effects

Taskinen e al. Gene Sel E ol (2017) 49:36 DOI 10.1186/s12711-017-0310-9 RESEARCH ARTICLE Single-s ep SNP-BLUP wi hon- he- ly impu ed geno ypes and esidual polygenic e ec s Ma i Taskinen* , Esa A. Män ysaa i and Ismo S andén Abs ac Backg ound: Single-s ep genomic bes linea unbiased p edic ion (BLUP) e alua ion combines ela ionship in o ma- ion om pedig ee and genomic ma ke da a. The inclusion o he genomic in o ma ion in o mixed model equa ions equi es he in e se o he combined ela ionship ma ix H , which has a dense ma ix block o geno yped animals. Me hods: To a oid in e sion o dense ma ices, single-s ep genomic BLUP can be ans o med o single-s ep single nucleo ide polymo phism BLUP (SNP-BLUP) which ha e obse ed and impu ed ma ke coe icien s. Simple block LDL ype decomposi ions o he single-s ep ela ionship ma ix H we e de i ed o ob ain di e en ypes o linea ly equi alen single-s ep genomic mixed model equa ions wi h di e en se s o epa ame ized andom e ec s. Fo non-geno yped animals, he impu ed ma ke coe icien e ms in he single-s ep SNP-BLUP we e calcula ed on- he- ly du ing he i e a i e solu ion using spa se ma ix decomposi ions wi hou s o ing he impu ed geno ypes. Residual polygenic e ec s we e added o geno yped animals and ansmi ed o non-geno yped animals using ela ionship coe icien s ha a e simila o impu ed geno ypes. The ela ionships we e u he o hogonalized o imp o e con e - gence o i e a i e me hods. Resul s: All p esen ed single-s ep SNP-BLUP models can be sol ed e icien ly using i e a i e me hods ha ely on i e a ion on da a and spa se ma ix app oaches. The e iciency, accu acy and i e a ion con e gence o he de i ed mixed model equa ions we e es ed wi h a small da ase ha included 73,579 animals o which 2885 we e geno yped wi h 37,526 SNPs. Conclusions: In e sion o he la ge and dense genomic ela ionship ma ix was a oided in single-s ep e alua ion by using ully o hogonalized single-s ep SNP-BLUP o mula ions. The numbe o i e a ions un il con e gence was smalle in single-s ep SNP-BLUP o mula ions han in he o iginal single-s ep GBLUP when he i abili y was low, bu inc eased abo e ha o he o iginal single-s ep when he i abili y was high. © The Au ho (s) 2017. This a icle is dis ibu ed unde he e ms o he C ea i e Commons A ibu ion 4.0 In e na ional License (h p://c ea i ecommons.o g/licenses/by/4.0/), which pe mi s un es ic ed use, dis ibu ion, and ep oduc ion in any medium, p o ided you gi e app op ia e c edi o he o iginal au ho (s) and he sou ce, p o ide a link o he C ea i e Commons license, and indica e i changes we e made. The C ea i e Commons Public Domain Dedica ion wai e (h p://c ea i ecommons.o g/ publicdomain/ze o/1.0/) applies o he da a made a ailable in his a icle, unless o he wise s a ed. Backg ound The i s model o simul aneously combine genomic in o ma ion wi h non-geno yped animal in o ma ion was single-s ep bes linea unbiased p edic ion (BLUP) [1, 2] o ssGBLUP. When he numbe o geno yped animals is la ge, ssGBLUP may become compu a ion- ally in easible o p ac ical pu poses because i equi es he in e ses o dense ma ices o sizeequal o he num- be o geno yped animals, pa icula ly he in e se o he genomic ela ionship ma ix G−1 g . In addi ion, ma ix Gg can be singula when he numbe o geno yped indi idu- als exceeds he numbe o ma ke s. Compu a ional chal- lenges may ha e been a eason o he slow adop ion o ssGBLUP ins ead o a mul i-s ep app oach. A compu a- ionally scalable al e na i e, he algo i hm o p o en and young (APY), has been sugges ed [3]. In APY, a spa se G−1 APY app oxima ion o he G−1 g ma ix is c ea ed by se - ing a diagonal ma ix o a g oup o indi iduals. In p ac- ice, APY has been able o educe compu a ional cos s signi ican ly when he numbe o geno yped animals is e y la ge [4, 5]. Howe e , di e en se s o co e animals Open Access G ene ics S elec ion E olu ion *Co espondence: ma i. askinen@luke. i Na u al Resou ces Ins i u e Finland (Luke), Mylly ie 1, Jokioinen, Finland Page 2 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 in APY will gi e di e en e alua ions, which may a ec selec ion decisions. An al e na i e o mula ion called he eina e single- s ep single nucleo ide polymo phism BLUP (ssSNP- BLUP) [6] o e comes some o he majo compu a ional challenges in ssGBLUP. In pa icula , he e is no need o cons uc o in e he genomic ela ionship ma ix Gg . The o iginal idea in ssSNP-BLUP ci cum en s he p ob- lems in ssGBLUP by p edic ing o impu ing geno ypes o non-geno yped animals, and elying on compu a- ionally less demanding SNP-based p edic ion ins ead o b eeding alue based p edic ion. An addi ional ad an age is ha he ma ke e ec solu ions a e easie o use o in e im p edic ions. The ssSNP-BLUP has some compu a ional challenges as well. A simple implemen a ion o ssSNP-BLUP gen- e a es and s o es geno ypes o all SNPs o all animals. This will lead o e y la ge disk s o age and as eading equi emen s ha will p ohibi use o he app oach o la ge popula ions. An al e na i e is o make he equi ed geno ype impu a ions “on- he- ly” ins ead o s o ing he e y la ge amoun o impu ed geno ypes o ile. Fo p ac- ical pu poses, he on- he- ly impu a ion equi es a as compu ing app oach o he impu a ion s ep and/o as con e gence o he i e a i e me hod. The o iginal desc ip ion o ssSNP-BLUP app oach p esen s a wide ange o models han ssGBLUP such as he use o a numbe o di e en Bayesian SNP model o - mula ions [6]. In ssGBLUP, in con as o ssSNP-BLUP, i is ypical o include esidual polygenic (RPG) in o ma- ion o enhance genomic in o ma ion by including ped- ig ee-based ela ionships in o he genomic ela ionship ma ix. Thus, he genomic ela ionship ma ix is usually “adjus ed” wi h pa o he pedig ee ela ionship ma ix ei he o supply mo e addi i e ela ionship in o ma- ion o simply o make he genomic ela ionship ma ix in e ible. In he o iginal ssSNP-BLUP o mula ion, his in o ma ion has no been included. The compu a ional challenges o ssGBLUP has led o he in oduc ion o se e al equi alen models (e.g., [7–9]). Howe e , hese al e na i e app oaches ha e had poo con e gence by i e a i e me hods [7, 9]. One eason is ha he co a i- ance s uc u es ha e a poo e condi ion numbe which is a a io o he la ges and smalles eigen alues and is used o measu e nume ical s abili y [9]. Some o he al e na- i e e sions o ssGBLUP ha e had SNP e ec s. Because ssSNP-BLUP is equi alen o ssGBLUP and ssGBLUP has many equi alen o ms, al e na i e o mula ions o mixed model equa ions (MME) can be de i ed o ssSNP-BLUP as well. In his pape , simple block LDL ype decomposi ions [10] o he ssGBLUP ela ionship ma ix a e de i ed o ob ain se e al linea ly equi alen MME. This allows de i a ion and es ing o se e al equi alen ssSNP-BLUP MME ha a oid making and s o ing he impu ed geno- ypes. In his pape , an explici impu a ion o geno ypes is no needed bu ins ead we apply spa se ma ix decom- posi ions o a ain pedig ee-based eg essions o e alu- a ions o geno yped animals on non-geno yped animals using spa se ma ix decomposi ions. Accu acy and i e a- ion con e gence o he de i ed equi alen MME a e es ed on a small No dic dai y ca le da ase . Me hods Fo any model, an in ini e numbe o equi alen mod- els exis . In he ollowing, i s we ecall he concep o equi alen models and how equi alen models can be made by a aching co a iance in o ma ion o he model design ma ix. As an example, equi alence o genomic BLUP (GBLUP) and SNP-BLUP is p esen ed. Then, ssG- BLUP is ecalled, and i s co a iance s uc u e o mula ed using LDL decomposi ion. Equi alen MME a e de i ed whe e in o ma ion om he genomic ela ionship ma ix a e a ached o he model design ma ix simila ly as shown o GBLUP and SNP-BLUP, esul ing in ssSNP- BLUP. Finally, MME ha ing o hogonalized andom e ec s o diagonal co a iance s uc u es a e de i ed in o de o imp o e he con e gence in i e a i e me hods. A small da ase is used o illus a e pe o mance o he de i ed equi alen models. Linea ly equi alen models Two mixed linea e ec s models, i.e. wi h he same obse a ions y bu di e en ixed ( b and  b ) and andom e ec s ( u and  u ), esidual e o s ( e and  e ), model ma ices ( X , Z ,  X , and  Z ) and a iance s uc- u es ( G , R ,  G , and  R ), a e said o be linea ly equi alen models [11–13] i he expec ed alues and he a iances o he obse a ions a e equal. Thus, models (1) and (2) a e equi alen i : Sepa a e equa ions o  ixed and andom e ec s Mixed model(1) can be sol ed om sepa a e equa ions o ixed and andom e ec s [14] as: whe e ma ix V needs o be in e ible. The size o ma ix V is he numbe o obse a ions, which can be e y (1) y=Xb +Zu +e,Va (u)=G,Va (e)=R (2) y=  X  b+  Z u+  e,Va ( u)=  G,Va ( e)=  R (3)  Xb =  X  b, V=ZGZ′+R=  Z  G  Z′+  R=  V . (4)  b =(X′V−1X) −1 X′V−1y  u=GZ′V− 1 (y−X  b) , Page 3 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 la ge, and, he e o e, sol ing he mixed model using his me hod is seldom easible in p ac ice. Mixed model equa ions In p ac ice, Eq.(4) can be sol ed by Hende son’s MME [14] as ollows: whe e he a iance ma ix o he andom e ec s G and he esidual a iance ma ix R need o be in e ible. MME(5) usually lead o mo e spa se ma ix sys ems han Eq.(4). Equi alen models byspli ing he a iance ma ix Suppose he a iance ma ix G can be exp essed as a ma ix p oduc as ollows: whe e M is ec angula and  G is an in e ible squa e ma ix. He e, he ma ix  G could be, o example, an iden i y ma ix and could ha e di e en dimensions han he G ma ix. Ma ices G and Z a e always oge he in Eq.(4). This allows us o e-pa ame ize he model as: and he ea e : whe e he equi alen model has he same esidual a i- ance ma ix (  R=R ), he new model ma ix  Z=ZM , and he andom e ec s  u a e de ined as: Now, acco ding o Eq. (3), he quan i ies wi h a ilde oge he wi h he same obse a ion ec o y , he o iginal ixed e ec s (  b=  b ), and design ma ix (  X=X ) o m a linea ly equi alen model(2). Linea ly equi alen MME O iginal ixed e ec s  b and he new andom e ec s  u o he linea ly equi alen model can be sol ed, simila ly as in Eq.(4), om: o om he co esponding MME(5): (5)  X ′ R −1 XX ′ R −1 Z Z ′ R − 1XZ ′ R − 1Z + G − 1  b  u = X ′ R −1 y Z ′ R − 1y , (6) G=M  GM′, (7) V =ZGZ ′ +R =ZM  GM′Z′+R =  Z  G  Z′+  R=  V, (8)  u =M  GM′Z′V−1(y−X  b) =M  G  Z′ V− 1 (y−X  b)=M u, (9)  u=  G  Z ′ V − 1(y − X  b) . (10)  b =(X′  V−1X) −1 X′  V−1y  u=  G  Z′ V− 1 (y−X  b), Because o he equi alence  u=M u in Eq.(8), he o igi- nal andom e ec s  u can be ob ained om he solu ion o he linea ly equi alen MME (11) by p e-mul iplying he new andom e ec s  u wi h ma ix M , i.e. whe e iden i y ma ix Ib has dimension o  b . No e ha he in e se o he MME ma ix in Eq.(12) is usually no e alua ed explici ly, bu a he , he co esponding linea ma ix equa ion is sol ed using ei he di ec o i e a i e solu ion me hods. The o iginal MME (5) can, hus, be sol ed om a modi ied linea ly equi alen MME (11) whe e he num- be o andom e ec s in  u could be smalle o la ge han in  u , he a iance s uc u e  G could be easie o ob ain o in e han he o iginal G , o he new ma ix sys em could be o he wise nume ically mo e e icien . P esen a ion o GBLUP andSNP‑BLUP As an example, conside single- ai GBLUP MME [15]. The a iance ma ix o andom e ec s u is based on a genomic ela ionship ma ix ( Gg ), which desc ibes he genomic ela ionships be ween indi iduals, i.e. Va (u) = σ 2 u G g . The genomic ela ionship ma ix is usu- ally ully dense and inc eases in o de as he numbe o geno yped animals inc eases. The in e ed genomic ela ionship ma ix Gg −1 is needed in he solu ion o he GBLUP MME: whe e a single ai case is assumed o simplici y, R= σ 2 eI , and = σ 2 e σ 2 u . The e a e many ways o cons uc Gg . Assume he genomic ela ionship ma ix Gg can be exp essed, in sim- pli ied o m, using (cen e ed and scaled) ma ke ma ix Zm [15] so ha : whe e and Im is an iden i y ma ix o size equal o he numbe o ma ke s. Now, a linea ly equi alen MME(11), al e - na i e o GBLUP MME(13), can be de i ed and o iginal e ec s sol ed simila ly as in(12) so ha : (11)  X′R−1XX ′R−1  Z  Z′R− 1 X Z′R− 1  Z+ G− 1  b  u = X′R−1y  Z ′ R − 1y . (12)  b  u  =  Ib0 0M  X′R−1XX ′R−1  X  Z′R−1X  Z′R−1  Z+  G−1 −1 X′R−1y  Z′R−1y , (13)  b  u = X′XX ′Z Z′XZ ′Z + G−1 g−1 X′y Z′y , (14) Gg= Z m Z ′ m= M  GM ′, (15) M=Zmand  G=Im, Page 4 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 whe e  Z=ZZm and  is he same a iance a io as in Eq. (13). This equi alen MME sys em, known as he SNP-BLUP [16], has ma ke s as andom e ec s ins ead o indi iduals. Random e ec s a e also “o hogonalized”, and, hus, in e sion o he dense genomic ela ionship ma ix Gg is a oided. No e ha , i he ma ke e ec s ha e unequal a iances, he ela ionship ma ix can be build as: whe e and he ma ix B is a diagonal co a iance ma ix desc ib- ing he a iances o di e en ma ke e ec s. Now he solu ion o SNP-BLUP becomes: Single‑s ep SNP‑BLUP In ssGBLUP [1, 2], some indi iduals ha e genomic in o - ma ion while some ha e only pedig ee in o ma ion. The model o ssGBLUP is a special case o mixed e ec mod- els, whe e he be ween-animal ela ionships a e modeled ia he agg ega ed ela ionship ma ix H [1, 2]. The ela- ionships in H a e desc ibed by he pedig ee-based ela- ionship ma ix A , and he genomic ela ionship ma ix among geno yped animals by Gg . Rela ionships among non-geno yped indi iduals a e cons uc ed om he ped- ig ee bu modi ied acco ding o he ela ionships among geno yped animals. Assuming ha he non-geno yped indi iduals a e deno ed wi h sub- and supe -sc ip s 1 and he geno yped indi iduals by sub- and supe -sc ip s 2, he pedig ee ela- ionship ma ix and i s in e se a e he ollowing: Assume, as be o e, ha he genomic ela ionship ma ix has he o m Gg= Z m Z ′ m . I he geno yped popula ion con ains iden ical wins, i.e. clones, o i he e a e mo e indi iduals han ma ke s, he genomic ela ionship ma ix Gg becomes singula and he usual MME(5) can- no be cons uc ed. Hence, Gg is commonly adjus ed by eg essing i owa ds he pedig ee ela ionship ma ix A22 wi h:  b  u = Ib0 0Zm  X′XX ′  Z  Z′X Z′ Z+  Im−1 X′y  Z ′ y , (16) Gg= Z m BZ ′ m= M  GM ′, (17) M=Zmand  G=B,  b  u = Ib0 0Zm  X′XX ′  Z  Z′X Z′ Z+  B− 1 −1 X′y  Z ′ y . (18) A= A11 A12 A21 A22  and A−1 = A 11 A 12 A21 A22 . (19) Gw=wA22 +(1−w)Gg, whe e w is a scala weigh be ween 0 and 1 and can be in e p e ed as he ela i e weigh on he polygenic e ec [2, 15]. Single‑s ep ela ionship ma ix The in e se o he ssGBLUP a iance ma ix H is [1, 2]: Using block LDL decomposi ion, i is equal o: whe e I1 and I2 a e iden i y ma ices o size equal o he numbe o non-geno yped and geno yped animals, espec i ely. Fo impu a ion o geno ypes: is a eg ession p edic ion o an impu a ion ope a o ha expands he genomic ela ionship in o ma ion om he geno yped o he non-geno yped indi iduals [2, 6, 17]. By in e ing he LDL decomposi ion o Eq.(22), he a iance ma ix H has a simila decomposi ion: By using block ma ix in e sion iden i ies (23) and o he “impu ed” A22 ma ix o Gw (19) he a iance ma ix H (24) can be al e na i ely exp essed as: whe e Gimp is an impu ed genomic ela ionship ma ix as ollows: No e ha his ope a es on Gg ins ead o Gw . Also, no e ha Gimp has a size equal o he numbe o all animals. Geno yped animals ha e obse ed ma ke da a bu non- geno yped animals ha e impu ed ma ke da a. (20) H −1 = A−1 +00 0G w− 1 − (A 22 ) − 1  (21) = A 11 A 12 A 21 A 22 +Gw− 1 −(A22)− 1 . (22) H −1 = I10 − A′ imp I2  A 11 0 0G w − 1  I1−Aimp 0I 2 , (23) Aimp = A12(A22) − 1 =− (A11) −1 A 12 (24) H= I1Aimp 0I 2  (A11) −1 0 0G w I10 A′ imp I2 . (25) (A 11 ) −1 =A11 −A12(A22)− 1 A21 (26)  Aimp I2  A22  A′ imp I2 = A − (A11) −1 0 00 , (27) H= (1 − w)  (A11) −1 0 00 + wA + (1 − w)Gimp , (28) G imp = Aimp I2  Gg  A′ imp I2 . Page 5 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 Equi alen single‑s ep MME In he ollowing, he single-s ep ela ionship ma ix H will be exp essed as six di e en decomposi ions equi a- len o Eq.(6): o di e en ma ices Mi and  Gi . All o hese lead o a linea ly equi alen ssGBLUP MME sys em o new se s o andom e ec s wi h, po en ially, di e en nume ical p ope ies. The linea ly equi alen MME o hese modi- ied se s o andom e ec s a e simila o Eq.(11) and he o iginal e ec s can be sol ed simila ly as in Eq.(12): whe e  Zi=ZMi and a single ai case is assumed. No e ha in hese equi alen MME, squa e ma ix  Gi ep esen s he co a iance s uc u e o he epa a- me ized andom e ec s, Mi hasa ow o each o iginal andom e ec in u o change he model ma ix Z , and  Zi is he ede ined model ma ix. In o de o de i e linea ly equi alen MME o mul iple ai cases: whe e ⊗ is he K onecke p oduc , he gene ic (co) a i- ance ma ix G0 o size numbe o ai s is assumed o ha e a decomposi ion: Fo example, a simple case would ha e M0 as iden i y ma ix and  G0 as G0 . Now he a iance ma ix G can be exp essed using he decomposi ions o he single-s ep ela ionship ma ix(29) simila ly as in Eq.(6): whe e E ec s can hen be sol ed om linea ly equi alen mul i- ple ai MME (12): whe e  Zi=Z(M0⊗Mi) . (29) H=Mi  GiM′ i ,i = 1, ..., 6,  b  u  =  Ib0 0M i  X′XX ′  Zi  Z′ iX Z′ i Zi+   G− 1 i−1 X′y  Z ′ i y , (30) G=G0⊗H, (31) G0=M0  G0M′ 0. (32) G=G0⊗H=(M0  G0M′ 0)⊗(Mi  GiM′ i) (33) =(M0⊗Mi)(  G0⊗  Gi)(M0⊗Mi)′ (34) =M  GM′, (35) M=M0⊗Miand  G=  G0⊗  Gi.   b  u  =  Ib0 0M 0⊗Mi  ×  X′R−1XX ′R−1  Zi  Z′ iR−1X  Z′ iR−1  Zi+  G−1 0⊗  G−1 i  −1  X′R−1y  Z′ iR−1y , Basic equi alen ssGBLUP MME The LDL decom- posi ion(24) can be used di ec ly o build he i s lin- ea ly equi alen ssGBLUP MME o his pape using H=M1  G1M′ 1 whe e: F om he modi ied ela ionship ma ix  G1 (36), i can be seen ha his basic equi alen ssGBLUP MME has andom e ec s o non-geno yped animals wi h a iances (A 11 ) −1 and o geno yped animals wi h a iances o he adjus ed genomic ela ionship ma ix Gw . The numbe o e ec s in he sys em is, hus, he same as in he o iginal ssGBLUP. Basic RPG ssSNP‑BLUP MME O he linea ly equi a- len ssGBLUP MME o he o m(29) can be de i ed as well. No e ha he adjus ed genomic ela ionship ma ix Gw in Eq.(19) can be exp essed by ma ix p oduc s as ollows: whe e Gg= Z m Z ′ m . The second linea ly equi alen ssGB- LUP MME can be buil by subs i u ing Gw in Eq. (38) o Eq. (36): In his o m, we a oid he in e se o Gg in he MME (11). The coe icien s w and (1−w) we e also spli using squa e oo s o ma ix M2 so ha he new a iance ma ix  G2 can be in e ed e en when w is 0 o 1. The i s g oup o new andom e ec s in Eq.(39), o he non-geno yped animals, is he same as in he i s , basic equi alen ssG- BLUP in Eq.(36). Howe e , he geno yped animals now ha e andom e ec s ela ed h ough he a iance ma ix A22 . E ec s in his second e ec g oup can be seen as esidual polygenic e ec s ha can desc ibe e ec s ha he ma ke e ec s a e unable o model [8]. The hi d g oup o andom e ec s a e he ma ke e ec s as in SNP- BLUP(15) and so, his decomposi ion(39) can be called basic RPG ssSNP‑BLUP MME. Compa ed o he o iginal ssGBLUP, he second equi alen MME has ma ke e ec s in addi ion o he animal e ec s. (36) M 1 = I1Aimp 0I 2  and  G1 = (A11) −1 0 0G w. (37) G w = I2I2  wA22 0 0(1 − w)Gg  I2 I2  (38) = I2Zm  wA22 0 0(1 − w)Im  I2 Z′ m, (39) M 2= � I1 √ wAimp √ 1−wAimpZm 0√wI2√1−wZm � � G2=  (A11)−100 0A 22 0 0 0I m   . Page 6 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 Expanded RPG ssSNP‑BLUP MME A hi d linea ly equi alen MME o o m(29) can be de i ed om he al e na i e exp ession o ma ix H in Eq. (27) and by spli ing Gg= Z m Z ′ m in Eq.(28) as: He e, ma ices E1 and E2 a e ec angula spa se inci- dence ma ices ha selec he subse s o non-geno yped and geno yped animals, espec i ely, om he A ma ix. Bo h E1 and E2 ha e he same numbe o columns, i.e. numbe o all animals. Ma ix E1 has a ow o each non- geno yped and ma ix E2 o each geno yped animal co - esponding o animal’s column in ma ices Z1 and Z2 o Each ow o bo h E1 and E2 has only one non-ze o ele- men , a alue one a he column co esponding o ha animal’s loca ion among all o he animals. Hence, when ows and columns o he ma ices o all animals a e in he same o de as in ma ix A in Eq.(18), ma ix  E1 E2  is an iden i y ma ix o he size o all animals. The hi d equi alen MME (40) has h ee g oups o e ec s simila o he second MME(39). The hi d e ec g oup has, again, he o hogonal ma ke e ec s, and so his o mula ion is a ssSNP-BLUP as well. The i s e ec g oup, o he non-geno yped animals has, howe e , a cons an mul iplie √1−w . Also, he second g oup, ela ed h ough he pedig ee ela ionship ma ix A , has now e ec s o all animals, and no jus o he geno yped animals. Thus, in his hi d expanded RPG ssSNP‑BLUP MME, he non-geno yped animals ha e wo se s o an- dom e ec s. Special cases o equi alen ssGBLUP MME These h ee equi alen MME, (36), (39) and (40), will app oach he usual animal model when w→1 . A he limi ( w=1 ), he expanded RPG ssSNP-BLUP MME 3 (40) has clea ly he ecognizable co a iance s uc u e o A . The basic equi - alen ssGBLUP MME1(36) and he basic RPG ssSNP- BLUP MME 2(39) a e models whe e he geno yped ani- mals ac as base animals and he non-geno yped animals a e eg essed on hem. Fo he o he di ec ion o w→0 , he basic equi alen ssGBLUP MME 1(36) con e ges o an al e na i e p es- en a ion o he s anda d ssGBLUP. Howe e , i di ides (40) M 3= �√ 1−wI1 √ wE1 √ 1−wAimpZm 0√wE2√1−wZm � � G3=  (A11)−100 0 A0 0 0I m   . (41) Z=[Z1Z2]. he b eeding alues o non-geno yped animals in o eg essions on geno yped animals and in o non-impu ed b eeding alues ha a e no condi ional on hem. Simi- la ly, a he limi w=0 he basic and he expanded RPG ssSNP-BLUP MME, 2(39) and 3(40), coincide wi h he simple ssSNP-BLUP wi hou esidual polygenic e ec s. Fo example, in he single ai case, he MME coe icien ma ix o he basic RPG ssSNP-BLUP MME2 (39) is whe e W=(Z 1 Aimp +Z 2 )Zm . In he case whe e w=0 , i.e. he e is no adjus men o he genomic ela ionship ma ix and, he e o e, no esid- ual polygenic e ec s, he basic and he expanded RPG ssSNP-BLUP MME, 2 (39) and 3(40), a e essen ially he same MME as was de i ed by Fe nando e al.[6]. The main di e ences a e ha hey ha e mo ed he cen e - ing e m o he ma ke ma ix Zm in o an addi ional ixed e ec , and hey p oposed o sol e he MME using Bayes- ian eg ession. E icien implemen a ion The h ee linea ly equi alen MME based on Eqs.(36), (39), and (40) con ain in e ed and non-in e ed e ms o he pedig ee ela ionship ma ix A . In an e icien se up o sol e hese MME, all hese e ms can be exp essed, o modi ied in o a o m ha can be exp essed, wi h spa se ma ices o spa se decomposi ions o spa se ma ices. This is expec ed o gi e h ee e icien implemen a ions o hese MME. In p ac ice, ma ix equa ions o hese MME a e assumed o be sol ed i e a i ely by he p econ- di ioned conjuga e g adien (PCG) algo i hm. Then, only a ma ix- ec o p oduc o he MME coe icien ma ix imes a ec o is pe o med once e e y i e a ion. Spa se ma ices anddecomposi ions The in e se o he pedig ee ela ionship ma ix A−1 can be exp essed e icien ly [18] as: whe e animals a e so ed in ( e e sed) age o de om he younges o he oldes using spa se pe mu a ion ma ix Q′ , so ha ma ix L= (I − 1 2 P) ′ D 1 2 becomes a lowe i- angula ma ix in Q′A−1Q=LL′ . The diagonal ma ix D has alues 4/(4−k−Fs) whe e k is he numbe o known pa en s and Fs is he sum o pa en inb eeding coe icien s. In he “pa en al ma ix” P on ow i , he e a e 1s in columns co esponding o pa en s o animal i . The pa en al ma ix can be in e p e ed, oge he wi h iden i y (42)     X ′ XX ′ Z1 0X ′ W Z′ 1XZ ′ 1Z1+A11 0Z ′ 1W 00A22−10 W′XW ′Z10W ′W + Im     , (43) A −1 = Q  I − 1 2P ′ D  I − 1 2P  Q′ = QLL′Q′ , Page 7 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 ma ix I , as a e y spa se lowe iangula “Cholesky” ma ix( L )[18]. The pedig ee ela ionship ma ix A can be exp essed as he in e se o i s in e se, and so he subma ices o he pedig ee ela ionship ma ix and i s in e se(18) can be ob ained by selec ing he app op ia e ows and columns as ollows: whe e i,j=1, 2 . Ma ix- ec o p oduc s Aijx and Aijx can be e icien ly compu ed using hese decomposi- ions [19, 20]. The subma ices Aij a e e y spa se, so hey could al e na i ely be exp essed as sepa a e spa se ma ices. The in e se o A−1 o pa icula pa s o i (e.g. (A 11 ) −1 ) a e, howe e , in gene al non-spa se and, hus, hese in e se ma ix e ms should ne e be compu ed explic- i ly. The spa se subma ix A11 can be exp essed using spa si y p ese ing Cholesky ac o iza ion so ha : whe e ma ix L1 is spa se lowe iangula and Q1 spa se pe mu a ion ma ix. No e ha ma ix L1 has o be compu ed explici ly, as opposed o ma ix L in Eqs. (43) o(46). Ma ix L1 has some mo e ill-ins compa ed o A11 bu is s ill e y spa se and e icien in use. In [4], i was demons a ed ha he compu a ions emain a o d- able e en when he da ase size g ows. Compu a ions in ol ing ma ix in e sions o A−1 o pa s o i can be ans o med in o solu ions o spa se ma ix equa ion sys ems [20]: whe e and 1 a e ec o s o app op ia e sizes in o which he in e se ma ix ope a ions a e pe o med. He e he backslash ( ) is an ope a o indica ing o wa d o back‑ wa d subs i u ions and emphasizes he impo ance o a oiding in e ing ma ices. In o he wo ds, L y is he solu ion x o equa ion Lx =y o can be exp essed as sol ing ( L , y ). No e ha he ma ix p oduc s a e ca e ully (44) A=(A− 1 )−1 =Q(L′)− 1 L− 1 Q′, (45) A ij = EiQ(L′) −1 L−1Q′E ′ j (46) Aij = EiQLL ′ Q ′ E ′ j, (47) A11 =Q1L1L′ 1Q′ 1, (48) A =Q(L′) −1 L−1Q′ = Q  L ′  L (Q ′ )  (49) ( A11) −1 1 = Q1  L ′ 1  L1 (Q ′ 1 1) , nes ed wi h pa en hesis so ha only ma ix- ec o ope a- ions a e pe o med. Fu he mo e, ollowing [4, 9] he in e se o he ma ix A22 , needed in he in e sion o he second modi ied ela ionship ma ix  G2 in Eq.(39), can be exp essed e i- cien ly using block ma ix in e sion iden i y simila o Eq.(25) as: whe e all e ms can be compu ed using Eqs. (45) o(49). On‑ he‑ ly impu a ion ope a ion The de i ed equi alen o mula ions in Eqs. (36), (39), and (40) con ain impu a ion ope a o Aimp =− (A11) −1 A 12 (23) in ma ices Mi ha a e needed when ope a ing wi h he modi ied model ma ix  Zi=ZMi and i s anspose  Z′ i= M ′ i Z ′ , and when calcu- la ing he o iginal andom e ec s in Eq.(12). When he MME a e sol ed by he PCG i e a ion algo i hm, he co e o he algo i hm is a mul iplica ion o he so-called di ec ion ec o by he le hand side o he MME(11). In his mul iplica ion, he impu a ion ope a o , as pa o he MME coe icien ma ix, ope a es ei he wi h a pa o he ec o pe aining o andom e ec s o he geno- yped animals (i.e. Aimp 2 ) o o a ec o o he ma ke e ec s m h ough he ma ke ma ix (i.e. AimpZm m ). Thus, in he anspose  Z′ i side, he impu a ion e m ope - a es on a ec o o size equal o he numbe o non-gen- o yped animals (i.e. A′ imp 1 ). In all cases, he size o he ec o e m ha ope a es on (A 11 ) −1 equals he numbe o non-geno yped animals, i.e. size o 1 . Fo example, he impu ed genomic ma ke da a e m, −(A 11 ) −1A 12 Zm m , ha expands he genomic in o ma ion om geno yped o non-geno yped animals, can be calcula ed using: and Eq.(49) so ha : Vec o 1 is calcula ed om Eq.(46), o wi hou con- s uc ing any o he ma ices by using ules o A−1 by pedig ee in o ma ion [18], o , al e na i ely, as a spa se ma ix- ec o p oduc o sepa a e spa se ma ix A12 . No e ha he ac ual impu a ion o he genomic ma ke in o ma ion is no needed. The impu a ion ope a ion is pe o med only implici ly, “on- he- ly” du ing he i e a- i e solu ion wi hou he need o use, o example, disk (50) (A22)− 1 =A 22 −A 21 (A 11 ) −1A 12 , (51)  2=Zm m (52) 1=A12 2 (53) − (A11) −1 A12Z m m=− Q1  L ′ 1  L1 (Q ′ 1 1) . Page 8 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 s o age. In he no mal impu a ion p ocess [6], he ma ke in o ma ion needs o be calcula ed o housands, o e en hund eds o housands ma ke ec o s o geno yped animals, i.e. columns o ma ke ma ix Zm . The p edic ed ma ke da a ma ix con ains eal numbe s and can be e y la ge, and, hus, akes a lo o ime and disk space o gene a e and use. In he on- he- ly impu a ion p ocess o gene ic e ec s, howe e , impu a ion is an ope a ion on a “p ojec ion ec o ” o he geno yped animals, i.e. a linea combina- ion o he ma ke ec o s. I needs o be pe o med only wice wi hin each i e a ion ound o each ai . Once o ma ix  Zi and ano he ime o he anspose  Z′ i in ma ix mul iplica ion o MME coe icien ma ix o Eq.(11). The impu a ion ope a ion is also needed once be o e he i e - a ion when calcula ing he igh -hand-side o he new andom e ec s (  Z′ i R − 1 y ) in Eq.(11) and once a he end o he i e a ion in o de o e ie e he o iginal andom e ec s (  u=Mi  u ) in Eq.(12). O hogonaliza ion o  andom e ec s In he SNP-BLUP e sions o he de i ed equi alen ssGBLUP, in Eqs. (39) and (40), he ma ke e ec s a e o hogonal, i.e. hei co a iance ma ix is diagonal. I u ns ou ha he PCG i e a ion numbe s o hese wo equi alen ssSNP-BLUP MME a e conside ably la ge han he o iginal ssGBLUP. In Eq. (39), he RPG e ec s and genomic alues p edic ed by SNPs ha e colinea i y, and in Eq.(40), he RPG and he ani- mal e ec s o non-geno yped animals a e di icul o sepa a e. The key o main aining good nume ical p ope ies o he o iginal ssGBLUP seems o be o “o hogonalize” he emaining new andom e ec s, oo. The emaining a i- ance s uc u es can be o hogonalized by spli ing he a iance ma ices and a aching he wo “hal s” in o he coe icien ma ices M as in Eqs.(6) and(29). The e m (A 11 ) −1 in Eqs.(39) and (40) can be o hog- onalized by using he spa se Cholesky ac o iza ion in Eq.(47) as ollows: whe e No e ha he pe mu a ion ope a o Q1 can be pe o med ou side he in e se ope a o . Howe e , he e m A22 in Eq.(39) seems o be much mo e di icul o decompose. S ill, i can be exp essed using he spa se decomposi ion o he ull ma ix A in Eq.(44) as: (54) (A 11 )−1 =M11  G11M′ 11, (55) M11 =Q1 ( L′ 1 ) −1 and  G11 =I1. (56) A22 =E2AE′ 2=M22  G22M′ 22, whe e and The ec angula ma ix  A 1 2 2 has dimensions numbe o geno yped indi iduals imes o al numbe o indi iduals, hence he ade-o he e is ha Eq.(56) will expand he second andom e ec g oup o Eq.(39) om geno yped o all indi iduals (size o I ). Using Eqs.(54), (56), and (58), he linea ly equi alen MME (39) and (40) can be “o hogonalized” in o ou h: and i h equi alen MME. Bo h o hese linea ly equi alen (29) ssSNP-BLUPs sha e he same o hogonal a iance s uc u e: and, hus, bo h also ha e he same numbe o new an- dom e ec s: andom e ec s o he geno yped indi idu- als, wo se s o andom e ec s o he non-geno yped indi iduals, and andom e ec s o he ma ke s. The di - e ence in equi alen MME 4 (59) and 5 (60) is on how hey di ide he RPG on non-geno yped animals. The ou h equi alen MME (59) can be called o hog‑ onal ssSNP‑BLUP MME, and he i h MME(60), o igi- na ing om he expanded RPG ssSNP-BLUP (40), o hogonal expanded ssSNP‑BLUP MME. Reduc ion o  he numbe o e ec s byusing ances o s o geno yped animals Ma ix A22 , as he co a iance s uc u e o he geno- yped animals, was epa ame ized in Eq.(56) using he ull pedig ee ela ionship ma ix A . This epa ame iza- ion inc eases he numbe o co esponding new an- dom e ec s om geno yped animals o all animals in he pedig ee. Howe e , compu a ions in ol ing A22 equi e only he geno yped indi iduals and hei ances o s. Thus, o educe he numbe o ex a new e ec s, A22 can be exp essed using a smalle pedig ee and ela ionship (57) M 22 =  A 1 2 2 and  G22 = I , (58)  A 1 2 k= E k Q(L′)−1,k =1, 2. (59) M 4=   M11 √wAimp � A 1 2 2 √ 1−wAimpZm 0√w � A 1 2 2 √1 − wZm  , (60) M 5=  √ 1−wM11 √w � A 1 2 1 √ 1−wAimpZm 0√w � A 1 2 2 √1 − wZm  , (61) � G 4= � G5=   I1 00 0I 0 0 0I m  , Page 9 o 15 Taskinen e al. Gene Sel E ol (2017) 49:36 ma ix  A con aining he geno yped animals and hei ances o s [4]. Le he in e se o he pedig ee ela ionship ma ix(43) o his smalle pedig ee be: whe e  Q and  L a e as be o e in Eq.(44) bu in ol e geno- yped animals and hei ances o s only. Ma ix A22 in Eq.(56) can hen be ep esen ed using he smalle pedi- g ee as: whe e  E2 selec s he geno yped indi iduals om he smalle pedig ee, and he size o iden i y ma ix Iganc is he num- be o geno yped animals and hei ances o s. The six h linea ly equi alen , educed o hogonal ssSNP‑BLUP MME, can be de i ed om Eqs.(59) and (61) as: He e only he non-geno yped ances o s o he geno- yped ha e wo se s o andom e ec s, all he o he non- geno yped and all geno yped animals ha e single se s o e ec s, in addi ion o he ma ke e ec s. Da a The de i ed MME we e es ed using a small No dic Red dai y ca le da ase and a simple model. The small da a- se and he model we e pa ially chosen in o de o be able o use di ec spa se ma ix solu ions o he o iginal ssGBLUP o ob ain accu a e “co ec solu ions”. The da a we e de eg essed p oo s o milk yield ha we e based on es ima ed b eeding alues om he No dic p oduc ion ai e alua ions by NAV (No dic E alua ions, Denma k). The e we e 73,579 animals in he pedig ee o which 2885 we e geno yped. Geno yped animals oge he wi h hei ances o s o m a smalle pedig ee o 6833 animals. The animals had been geno yped wi h he Illumina Bo ine SNP50 Bead Chip (Illumina, San Diego, USA). The analy- sis used 37,526 SNPs ha passed quali y con ol. The e we e 66,426 non-geno yped and 1222 geno yped animals wi h pheno ypes. Hence, 1663 animals had a geno ype bu no pheno ype. We conside ed a single ai model (62)  A− 1 =  Q  L  L′  Q′, (63) A22 =  E 2  A  E ′ 2=  M 22  G 22  M ′ 22, (64)  M 22 =  E2  Q(  L ′ ) −1 and  G22 = I ganc, (65) M 6= � M11 √ wAimp � M22 √ 1−wAimpZm 0√w� M22 √1−wZm � � G6=  I100 0I ganc 0 0 0I m   . and assumed a he i abili y o 0.5. The genomic da a con- ained one pai o animals wi h iden ical genomic ma ke da a and a couple o mo e nea iden ical pai s ha led o p oblems o he in e sion o he genomic ela ionship ma ix Gg wi hou A22 adjus men , i.e. in he case w=0 . Compa ison s a is ics The o iginal ssGBLUP wi h in e se a iance ma ix H−1 o Eq.(20) and he six linea ly equi alen o mula ions o he o m(29), om Eqs.(36), (39), (40), (59), (60), (61), and (65), we e implemen ed and es ed in an Oc a e[21] en i onmen . Spa se ma ix ac o iza ions we e based on CHOLMOD ou ines[22]. Six di e en weigh s w we e es ed: 0.00, 0.01, 0.10, 0.20, 0.30, and 1.00. Because o he singula i y in he in e se o he genomic ela ion- ship ma ix Gg −1 wi h he es da a, he o iginal ssGB- LUP ma ix (20) and he i s equi alen o mula ion(36) we e no calcula ed when w=0 . The de i ed new o - mula ions we e compa ed agains he o iginal ssGBLUP, mainly ocusing on e iciency, accu acy, and numbe o i e a ions. E iciency o he de i ed equi alen ssGBLUP o mu- la ions elies on he spa si y o he pedig ee ela ionship ma ices and hei decomposi ions. In e se ma ix ope a- ions o hese spa se ma ices we e ans o med in o sol - ing spa se lowe iangula ma ix sys ems. The e iciency o hese sol ing ope a ions depend on he spa si y s uc- u e o he ma ices, i.e. numbe o non-ze o elemen s. Accu acy o he o mula ions was es ed by sol ing he MME wi h he PCG me hod using Oc a e’s PCG ou ine (pcg) wi h he diagonal o he MME coe icien ma ix as he p econdi ione , o wi hou p econdi ioning. Con- e gence ole ance in pcg was ela i e esidual no m. The ole ance was chosen o be small ( 10−12 ) so ha all sol ed e ec s, wi hou doub , con e ged. Accu acies, o a he he di e ences om he “exac solu ion”, we e cal- cula ed as ela i e esidual e o s (  ei ) be ween i e a i ely ob ained MME solu ions ( si ) and he di ec solu ion ( sdi ec ) o he o iginal ssGBLUP: whe e subsc ip i=1, ...,6 is he o mula ion numbe . Implemen a ions o he de i ed equi alen ssGBLUP o mula ion we e no ye s eamlined o speed and, hus, he execu ion imes we e nei he op imal, no compa a- ble. The pe o mance o he o mula ions is, he e o e, es ed by compa ing he numbe o i e a ions o he i e a- i e solu ion. The pu pose was o demons a e ha he i e a ion coun s a e compa able o hose ob ained by he o iginal ssGBLUP. Wi h a la ge numbe o geno yped indi iduals, he in e sion o he genomic ela ionship ma ix in he o iginal ssGBLUP becomes a bo leneck (66)  ei = � s di ec − s i� �sdi ec �,