2017年我们在《华南农业大学学报》上发表了名为《基于空间效应与竞争效应的林木遗传分析模型》(http://xuebao.scau.edu.cn/zr/hnny_zr/ch/reader/view_abstract.aspx?file_no=20170513&flag=1)的论文。该文中提出了一种新的个体(或动物)模型,该模型涵盖了空间残差效应和遗传竞争效应。我们知道,如果一棵树周围都存在树木,则会有东西南北向和东南、东北、西南和西北向等8棵树,那么理论上,我们也可以拟合这8个方位的遗传竞争效应。这在ASReml中可以实现。本文将简单示范如何拟合东向和西向的遗传竞争效应,并假设它们相等。
简单示范的代码如下:
library(asreml)
library(AAfun4)
read.example(package = "RTDat1",setpath=TRUE)
sp.comp<-read.file(file='sp.comp.csv',header=TRUE)
数据集sp.comp的结构如下:
> str(sp.comp)
'data.frame': 1200 obs. of 25 variables:
$ ID : Factor w/ 1200 levels "10001","10002",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Dad : Factor w/ 10 levels "1","2","3","4",..: 5 3 2 1 6 5 8 4 3 10 ...
$ Mum : Factor w/ 10 levels "11","12","13",..: 4 8 6 3 2 8 1 3 6 2 ...
$ Row : Factor w/ 30 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Col : Factor w/ 40 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Fam : Factor w/ 100 levels "F1","F10","F100",..: 31 73 50 15 10 75 90 19 51 5 ...
$ Rep : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
$ Plot : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
$ NN : Factor w/ 1160 levels "10001","10002",..: NA NA NA NA NA NA NA NA NA NA ...
$ SN : Factor w/ 1160 levels "10041","10042",..: 1 2 3 4 5 6 7 8 9 10 ...
$ WN : Factor w/ 1170 levels "10001","10002",..: NA 1 2 3 4 5 6 7 8 9 ...
$ EN : Factor w/ 1170 levels "10002","10003",..: 1 2 3 4 5 6 7 8 9 10 ...
$ WSN : Factor w/ 1131 levels "10041","10042",..: NA 1 2 3 4 5 6 7 8 9 ...
$ ESN : Factor w/ 1131 levels "10042","10043",..: 1 2 3 4 5 6 7 8 9 10 ...
$ WNN : Factor w/ 1131 levels "10001","10002",..: NA NA NA NA NA NA NA NA NA NA ...
$ ENN : Factor w/ 1131 levels "10002","10003",..: NA NA NA NA NA NA NA NA NA NA ...
$ y : num 13.7 11.3 13.1 16.8 15.5 ...
$ y.NN : num NA NA NA NA NA NA NA NA NA NA ...
$ y.SN : num 16.89 12.71 11.63 9.92 17.31 ...
$ y.WN : num NA 13.7 11.3 13.1 16.8 ...
$ y.EN : num 11.3 13.1 16.8 15.5 11.8 ...
$ y.WSN: num NA 16.89 12.71 11.63 9.92 ...
$ y.ESN: num 12.71 11.63 9.92 17.31 12.17 ...
$ y.WNN: num NA NA NA NA NA NA NA NA NA NA ...
$ y.ENN: num NA NA NA NA NA NA NA NA NA NA ...
需要处理谱系矩阵以及东西向与谱系矩阵的关联:
scped<-sp.comp[,1:3] # pedigree
scpedinv<-ainverse(scped)
temp<-as.character(sp.comp$EN) # eastern direction
sp.comp$EN <- factor(temp,levels=attr(scpedinv,'rowNames'))
temp<-as.character(sp.comp$WN) # western direction
sp.comp$WN <- factor(temp,levels=attr(scpedinv,'rowNames'))
空间效应和东西向遗传效应的联合模型的拟合代码如下:
## model4: saptial +competition--SCM
m4<-asreml(y~1,random=~vm(ID,scpedinv)+vm(EN,scpedinv)+and(vm(WN,scpedinv)),
na.action=na.method(x='include'),
residual=~ar1(Row):ar1(Col),data=sp.comp)
> summary(m4)$varcomp
component std.error z.ratio bound %ch
vm(ID, scpedinv) 1.71379455 0.66108857 2.5923827 P 0.2
vm(EN, scpedinv) 0.08159341 0.07141173 1.1425771 P 0.2
Row:Col!R 6.76225030 0.48157228 14.0420256 P 0.0
Row:Col!Row!cor 0.20658643 0.03437369 6.0100157 U 0.0
Row:Col!Col!cor 0.02313905 0.03321758 0.6965906 U 0.1
还可以使用另一种拟合方法,代码如下:
m4b<-update(m4,random=~str(~vm(ID,scpedinv)+vm(EN,scpedinv)+and(vm(WN,scpedinv)),
~us(2):vm(ID,scpedinv))+units)
> summary(m4b)$varcomp[,1:3]
component std.error z.ratio
vm(ID, scpedinv)+vm(EN, scpedinv)+and(vm(WN, scpedinv))!us(2)_1:1 1.93997475 0.62897649 3.0843359
vm(ID, scpedinv)+vm(EN, scpedinv)+and(vm(WN, scpedinv))!us(2)_2:1 0.11937849 0.14176756 0.8420720
vm(ID, scpedinv)+vm(EN, scpedinv)+and(vm(WN, scpedinv))!us(2)_2:2 0.09149699 0.06768029 1.3519001
units 3.85835010 1.61900221 2.3831654
Row:Col!R 4.45572072 2.38348663 1.8694129
Row:Col!Row!cor 0.37574039 0.16017661 2.3457881
Row:Col!Col!cor 0.03315040 0.06447171 0.5141852
想要拟合空间效应与竞争效应的个体模型,最重要的是需要先处理数据集的结构,它不同于一般的数据。