주어지는 자료는 6열까지이다. fullsib 효과야 2열과 3열을 붙이면 되는 거니까 어렵지 않다. 문제는 동기축을 자료로 붙이는 일이다. 그리고 리넘된 자료로 해야 하므로 renumf90을 적당히 돌리고, 같은 펜에 있는 동기축의 renumber된 번호를 붙여야 한다. 여기서는 일단 수동으로 자료를 준비하여 분석한다.
sex 효과, direct genetic 효과는 일반적이다. 그러나 associative genetic 효과의 design matrix는 일반적인 design matrix와는 다르다. 일반적인 design matrix는 한 행에 1이 하나만 들어가나, associative genetic 효과의 design matrix는 한 행에 같은 pen(group)에 있는 동기축들을 1로 해야 한다. 그래서 7열은 레벨을 만들지 말고, 8열로만 레벨을 만들어 한 효과로 처리하고 동기축을 1로 한다는 의미이다.
실행화면
blupf90 blupf90_social.par | tee blupf90_social_01.log
로그
blupf90_social.par
BLUPF90 ver. 1.68
Parameter file: blupf90_social.par
Data file: social2_data.txt
Number of Traits 1
Number of Effects 6
Position of Observations 6
Position of Weight (1) 0
Value of Missing Trait/Observation 0
EFFECTS
# type position (2) levels [positions for nested]
1 cross-classified 5 2
2 cross-classified 1 15
3 cross-classified 7 0
4 cross-classified 8 15
5 cross-classified 9 3
6 cross-classified 4 3
Residual (co)variance Matrix
48.480
correlated random effects 2 3
Type of Random Effect: additive animal
Pedigree File: social_pedi.txt
trait effect (CO)VARIANCES
1 2 25.70 2.250
1 3 2.250 3.600
Random Effect(s) 5
Type of Random Effect: diagonal
trait effect (CO)VARIANCES
1 5 12.50
Random Effect(s) 6
Type of Random Effect: diagonal
trait effect (CO)VARIANCES
1 6 12.12
REMARKS
(1) Weight position 0 means no weights utilized
(2) Effect positions of 0 for some effects and traits means that such
effects are missing for specified traits
* The limited number of OpenMP threads = 4
* solving method (default=PCG):FSPAK
Data record length = 9
# equations = 38
G
25.700 2.2500
2.2500 3.6000
G
12.500
G
12.120
read 9 records in 0.1406250 s, 139
nonzeroes
read 15 additive pedigrees
finished peds in 0.1406250 s, 250 nonzeroes
solutions stored in file: "solutions"
여기서는 associative effect 는 없고, sex, animal, pen, common environment effect(full-sib) 효과만 있는 것으로 한다.
# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
COMBINE 7 2 3
DATAFILE
social1_data.txt
TRAITS
6
FIELDS_PASSED TO OUTPUT
WEIGHT(S)
RESIDUAL_VARIANCE
48.48
EFFECT
5 cross numer
EFFECT
1 cross alpha
RANDOM
animal
FILE
social1_pedi.txt
FILE_POS
1 2 3
PED_DEPTH
0
(CO)VARIANCES
25.7
EFFECT
4 cross numer
RANDOM
diagonal
(CO)VARIANCES
12.12
EFFECT
7 cross alpha
RANDOM
diagonal
(CO)VARIANCES
12.5
OPTION sol se
Effect group 1 of column 1 with 2 levels, effect # 1
Value # consecutive number
1 4 1
2 5 2
Effect group 3 of column 1 with 3 levels, effect # 3
Value # consecutive number
1 3 1
2 3 2
3 3 3
Effect group 4 of column 1 with 3 levels, effect # 4
Value # consecutive number
a1a4 3 1
a2a5 3 2
a3a6 3 3
single-stepGBLUP의 MME에는 H 행렬의 역행렬이 등장한다. 마치 NRM과 GRM의 적당한 결합처럼 보인다. 아니 무식해서 그렇게 생각했다. 그러나 NRM과 GRM의 적당한 결합이 아니다.
다시 GBLUP을 생각해 보자. GBLUP은 SNP genotyping한 개체들의 GRM을 만들어서 개체들의 GEBV를 추정한다. 집단의 모든 개체들의 GEBV를 추정하려면 모든 개체를 SNP genotyping 하여야 한다. 그러나 한정된 예산 때문에 또는 이미 개체들이 죽어서 genotyping을 할 수 없다. 그러면 죽은 선조까지 포함하여 전체 집단에 대해서 GBLUP을 할 수는 없는 걸까? 이런 질문에서 출발하여 일부 유전체 분석을 한 개체들의 gene content를 이용하여 유전체 분석을 하지 않은 개체(non genotyped animals)들의 gene content를 추정하여 GBLUP을 하자는 것이 ssGBLUP이고 그때 H matrix를 이용하는 것이다.
유전체 분석을 한 개체는 유전체 분석 결과를 이용하고, 유전체 분석을 안 한 개체는 유전체 분석을 한 개체로부터 gene content를 추정하여 만든 GRM이 바로 H matrix이다.
genotyped animals가 주어졌을 때 non genotyped animals의 gene contents의 분포를 이용하는 것이므로 정규 분포의 조건부 분포에 대한 이해가 필요할 수도 있다. 키워드로 conditional normal distribution, mean vector of conditional normal distribution, covariance matrix of conditional normal distribution 등을 이용할 수 있다.
H matrix는 매우 복잡하나, H 역행렬은 매우 간단한 형태를 취하고 있다. 놀라운 자연의 법칙이고, 이걸 알아내는 과학자들도 대단하다. Simple is beautiful.
# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
DATAFILE
ssgblup_data.txt
TRAITS
6
FIELDS_PASSED TO OUTPUT
WEIGHT(S)
7
RESIDUAL_VARIANCE
245
EFFECT
4 cross alpha
EFFECT
1 cross alpha
RANDOM
animal
FILE
ssgblup_pedi.txt
FILE_POS
1 2 3
SNP_FILE
ssgblup_snpgenotype.txt
PED_DEPTH
0
(CO)VARIANCES
35.241
OPTION no_quality_control
OPTION AlphaBeta 0.95 0.05
OPTION tunedG 0
OPTION thrStopCorAG 0.10
OPTION solv_method FSPAK
Renumf90 실행
renumf90 renumf90_ssgblup.par | tee renumf90_ssgblup.log
Renumf90 실행 화면
Renumf90 로그
RENUMF90 version 1.145
renumf90_ssgblup.par
datafile:ssgblup_data.txt
traits: 6
R
245.0
Processing effect 1 of type cross
item_kind=alpha
Processing effect 2 of type cross
item_kind=alpha
pedigree file name "ssgblup_pedi.txt"
positions of animal, sire, dam, alternate dam, yob, and group 1 2 3 0 0 0 0
SNP file name "ssgblup_snpgenotype.txt"
all pedigrees to be included
Reading (CO)VARIANCES: 1 x 1
Maximum size of character fields: 20
Maximum size of record (max_string_readline): 800
Maximum number of fields for input file (max_field_readline): 100
Pedigree search method (ped_search): convention
Order of pedigree animals (animal_order): default
Order of UPG (upg_order): default
Missing observation code (missing): 0
hash tables for effects set up
first 3 lines of the data file (up to 70 characters)
13 0 0 1 558 9 0.001792115
14 0 0 1 722 13.4 0.001385042
15 13 4 1 300 12.7 0.003333333
read 10 records
table with 1 elements sorted
added count
Effect group 1 of column 1 with 1 levels
table expanded from 10000 to 10000 records
added count
Effect group 2 of column 1 with 10 levels
wrote statistics in file "renf90.tables"
Basic statistics for input data (missing value code is '0')
Pos Min Max Mean SD N
6 4.8000 15.400 9.5500 3.3890 10
random effect with SNPs 2
type: animal
file: ssgblup_snpgenotype.txt
read SNPs 9 records
Effect group 2 of column 1 with 14 levels
random effect 2
type:animal
opened output pedigree file "renadd02.ped"
read 26 pedigree records
loaded 12 parent(s) in round 0
Pedigree checks
Number of animals with records = 10
Number of animals with genotypes = 9
Number of animals with records or genotypes = 14
Number of animals with genotypes and no records = 4
Number of parents without records or genotypes = 12
Total number of animals = 26
Wrote cross reference IDs for SNP file "ssgblup_snpgenotype.txt_XrefID"
Wrote parameter file "renf90.par"
Wrote renumbered data "renf90.dat" 10 records
본 예제의 genotype은 제대로 된 genotype이 아니다. 그래서 다음 프로그램이 정상적으로 동작하지 않는다. 그래서 no_quality_control, thrStopCorAG 등의 옵션을 넣어서 강제로 실행을 하는 것이다. 실제 분석에서는 세심한 주의를 기울여서 옵션을 주고 문제점이 있나 없나 확인을 해야 한다. 그렇지 않으면 유전체 자료를 이용하는 것이 독이 될 수도 있다.
중간 과정을 점검하기 위하여 H 역행렬, A 역행렬, GimA22i 등등이 필요하면 위의 "OPTION solv_method FSPAK"을 지우고 다음을 추가하여 pregsf90을 실행한다.
실제에서는 pregsf90을 실행하여 유전체 자료의 문제점이 있는지 없는지 반드시 확인하여야 한다. pregsf90은 문제점을 찾아내고 수정하기 위한 많은 옵션을 제공하는데 반드시 의미하는 바가 무엇인지 확인하고 진행하기를 바란다. 그리고 pregsf90, blupf90은 문제가 있어도 대충 덮고 실행을 계속하게 한다. 그러나 그렇게 하지 말기를 바란다. 문제가 발견되면 문제를 일으킨 유전체 자료가 무엇인지 확인하고 분석비용이 아깝긴 하지만 넣는 것이 바람직한 것인지 아니면 버리는 것이 나은 것인지 판단을 하고 진행해야 한다. 보통 버리는 것이 정신 건강상 좋다. 유전체 자료를 넣는다고 무조건 정확도가 올라가서 개량량이 증가하는 것은 아니다. 양질의 유전체 자료를 넣을 때만 그렇게 된다. 당연한 얘기가 실제에선 무시되고 마구 분석하는 현실이 안타깝기만 하다. 모든 것이 그렇듯이 유전체 자료를 이용한 유전체 선발도 마법의 도구가 아니다.
blupf90 실행하기 - renumf90이 만들어낸 파라미터 파일을 이용하여 blupf90 실행
다음 명령어로 실행
blupf90 renf90.par | tee blupf90_ssgblup.log
실행 화면
실행 로그
renf90.par
BLUPF90 ver. 1.68
Parameter file: renf90.par
Data file: renf90.dat
Number of Traits 1
Number of Effects 2
Position of Observations 1
Position of Weight (1) 2
Value of Missing Trait/Observation 0
EFFECTS
# type position (2) levels [positions for nested]
1 cross-classified 3 1
2 cross-classified 4 26
Residual (co)variance Matrix
245.00
Random Effect(s) 2
Type of Random Effect: additive animal
Pedigree File: renadd02.ped
trait effect (CO)VARIANCES
1 2 35.24
REMARKS
(1) Weight position 0 means no weights utilized
(2) Effect positions of 0 for some effects and traits means that such
effects are missing for specified traits
* The limited number of OpenMP threads = 4
* solving method (default=PCG):FSPAK
Options read from parameter file for genomic
* SNP format: BLUPF90 standard (text)
* SNP file: ssgblup_snpgenotype.txt
* SNP Xref file: ssgblup_snpgenotype.txt_XrefID
* No Quality Control Checks !!!!! (default .false.): T
* Set the threshold to Stop the analysis if low cor(A22,G) (default 0.3): 0.1
* Create a tuned G (default = 2): 0
* AlphaBeta defaults alpha=0.95, beta=0.05) : 0.95 0.05
Data record length = 4
# equations = 27
G
35.241
read 10 records in 0.9531250 s, 21
nonzeroes
read 26 additive pedigrees
*--------------------------------------------------------------*
* Genomic Library: Dist Version 1.281 *
* *
* Optimized OpenMP Version - 4 threads *
* *
* Modified relationship matrix (H) created for effect: 2 *
*--------------------------------------------------------------*
Read 26 animals from pedigree file: "D:\users\bhpark\2021\job\공부하기\01_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode\26_ssGBLUP_blupf90\renadd02.ped"
Number of Genotyped Animals: 9
Creating A22
Extracting subset of: 19 pedigrees from: 26 elapsed time: 0.0000
Calculating A22 Matrix by Colleau OpenMP...elapsed time: .0000
Numbers of threads=1 4
Reading SNP file
Column position in file for the first marker: 4
Format to read SNP file: (3x,400000i1)
Number of SNPs: 50
Format: integer genotypes (0 to 5) to double-precision array
Number of Genotyped animals: 9
Reading SNP file elapsed time: .00
Statistics of alleles frequencies in the current population
N: 50
Mean: 0.074
Min: 0.000
Max: 0.944
Var: 0.038
Quality Control - Monomorphic SNPs Exist - NOT REMOVED: 40
Genotypes missings (%): 0.000
Calculating G Matrix
Dgemm MKL #threads= 1 4 Elapsed omp_get_time: 0.0000
Scale by Sum(2pq). Average: 3.19135802469136
Detecting samples with similar genotypes
elapsed time= 0.0
Blend G as alpha*G + beta*A22: (alpha,beta) 0.950 0.050
Frequency - Diagonal of G
N: 9
Mean: 1.057
Min: 0.565
Max: 2.648
Range: 0.104
Class: 20
#Class Class Count
1 0.5645 1
2 0.6687 1
3 0.7729 1
4 0.8771 4
5 0.9813 1
6 1.085 0
7 1.190 0
8 1.294 0
9 1.398 0
10 1.502 0
11 1.606 0
12 1.711 0
13 1.815 0
14 1.919 0
15 2.023 0
16 2.127 0
17 2.232 0
18 2.336 0
19 2.440 0
20 2.544 1
21 2.648 0
Check for diagonal of genomic relationship matrix
** High Diagonal of genotype 4 2.65 Not Removed
** Low Diagonal of genotype 8 0.56 Not Removed
Check for diagonal of genomic relationship matrix, genotypes not removed: 2
------------------------------
Final Pedigree-Based Matrix
------------------------------
Statistic of Rel. Matrix A22
N Mean Min Max Var
Diagonal 9 1.000 1.000 1.000 0.000
Off-diagonal 72 0.201 0.000 0.500 0.013
----------------------
Final Genomic Matrix
----------------------
Statistic of Genomic Matrix
N Mean Min Max Var
Diagonal 9 1.057 0.565 2.648 0.377
Off-diagonal 72 -0.116 -0.742 0.725 0.142
Correlation of Genomic Inbreeding and Pedigree Inbreeding
Variance of Y is Zero !!
Correlation: 0.0000
Diagonal elements
Estimating Regression Coefficients G = b0 11' + b1 A + e
Regression coefficients b0 b1 = NaN NaN
Correlation diagonal elements G & A NaN
All elements - Diagonal / Off-Diagonal
Estimating Regression Coefficients G = b0 11' + b1 A + e
Regression coefficients b0 b1 = -0.395 1.413
Correlation all elements G & A 0.708
Off-Diagonal
Using 56 elements from A22 >= .02000
Estimating Regression Coefficients G = b0 11' + b1 A + e
Regression coefficients b0 b1 = -0.483 1.647
Correlation Off-Diagonal elements G & A 0.212
***********************************************************************
* CORRELATION FOR OFF-DIAGONALS G & A22 IS LOW THAN 0.50 !!!!! *
* MISIDENTIFIED GENOMIC SAMPLES OR POOR QUALITY GENOMIC DATA *
***********************************************************************
Creating A22-inverse
Inverse LAPACK MKL dpotrf/i #threads= 1 4 Elapsed omp_get_time: 0.0000
----------------------
Final A22 Inv Matrix
----------------------
Statistic of Inv. Rel. Matrix A22
N Mean Min Max Var
Diagonal 9 1.233 1.000 1.429 0.017
Off-diagonal 72 -0.101 -0.571 0.000 0.009
Creating G-inverse
Inverse LAPACK MKL dpotrf/i #threads= 1 4 Elapsed omp_get_time: 0.0000
--------------------------
Final Genomic Inv Matrix
--------------------------
Statistic of Inv. Genomic Matrix
N Mean Min Max Var
Diagonal 9 5.227 1.593 11.227 9.161
Off-diagonal 72 0.308 -5.482 3.166 3.142
Check for diagonal of Inverse Genomic - Inverse of pedigree relationship matrix
------------------------------
Final G Inv - A22 Inv Matrix
------------------------------
Statistic of Inv. Genomic - A22 Matrix
N Mean Min Max Var
Diagonal 9 3.994 0.593 10.021 8.878
Off-diagonal 72 0.408 -5.355 3.293 3.069
*--------------------------------------------------*
* Setup Genomic Done !!!, elapsed time: 0.049 *
*--------------------------------------------------*
finished peds in 1.156250 s, 108 nonzeroes
solutions stored in file: "solutions"
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.6 p192
single step GBLUP(ssGBLUP)은 유전체 자료를 가지고 있는 개체와 가지고 있지 않은 개체를 동시에 분석할 수 있는 방법이다. 또한 SNP efffect를 구한 다음 DGV 또는 GEBV를 구하는 것이 아니라 SNP effect를 구하는 단계를 생략하고 한 번에 유전체 자료를 가지고 있는 개체는 GEBV, 없는 개체는 EBV를 구한다.
ssGBLUP에 포함하는 개체는 다음의 세 가지로 분류할 수 있다.
1. 능력검정 기록만 가지고 있는 개체 - 유전체 자료를 생성하기 전에, 즉 과거에 유전평가에 포함된 개체
2. 능력검정 기록과 유전체 자료를 가지고 있는 개체 - 소위 참조집단(reference population)이라고 하는 개체들로 GEBV를 추정한다. 참조집단이 많아야 유전체 선발(genomic selection)의 효율(아래 3번 개체들의 정확도 향상)이 높아진다.
3. 유전체 자료만 가지고 있는 개체 - 보통 어린 개체들로 아직 능력검정을 할 단계에 이르지 못한 개체. 이들에 대해 GEBV를 계산하고 조기 선발함으로써 세대간격(generation interval)을 줄여 유전적 개량량을 높일 수 있다. 이들은 시간이 지나 능력검정할 단계에 다다르면 능력검정을 하여 표현형 자료를 확보하여 참조집단이 된다.
결과가 책과 다르다. 책 187쪽에는 5세대의 혈통을 이용하여 구한 A(NRM)가 나오는데 책에서 제시한 혈통을 가지고는 이걸 구할 수가 없다. 본 예제에서 책 187쪽의 혈통을 사용한다고 하였는데 프로그램의 흐름상 사용할 수가 없었다. 이런 이유로 해서 같은 해가 나올 수 없는 것으로 보인다. 본 예제를 위의 참고 포스팅에서는 blupf90으로 푸는데 역시 해가 다르다. 그 이유는 참고 포스팅에서는 EDC를 가중치로 줄 때 그대로 사용하였다. EDC가 아니라 EDC의 역수를 가중치로 주어야 한다. EDC의 역수를 가중치로 주었을 때 동일한 해를 구할 수 있었다.
