パッケージの読み込み

library(gaston)

データの準備

McCouchら（2016; Nature Communication7:10532） が用いたイネのゲノムデータを用いて解析を行う。

Geno <- read.vcf( "/home/kimura-biomet/Desktop/Trials/data/Data_Rice/HDRA-G6-4-RDP1-RDP2-NIAS.AGCT_MAF005MIS001NR_imputed.vcf.gz" )
## ped stats and snps stats have been set. 
## 'p' has been set. 
## 'mu' and 'sigma' have been set.
G_mat <- as.matrix(Geno) - 1
Pheno <- read.csv( "/home/kimura-biomet/Desktop/Trials/data/Data_Rice/RiceDiversityPheno4GWASGS.csv", row.names = "X" )
dim(G_mat)
## [1]  1568 38769
colnames(Pheno)
##  [1] "Flowering.time.at.Arkansas"       "Flowering.time.at.Faridpur"      
##  [3] "Flowering.time.at.Aberdeen"       "FT.ratio.of.Arkansas.Aberdeen"   
##  [5] "FT.ratio.of.Faridpur.Aberdeen"    "Culm.habit"                      
##  [7] "Leaf.pubescence"                  "Flag.leaf.length"                
##  [9] "Flag.leaf.width"                  "Awn.presence"                    
## [11] "Panicle.number.per.plant"         "Plant.height"                    
## [13] "Panicle.length"                   "Primary.panicle.branch.number"   
## [15] "Seed.number.per.panicle"          "Florets.per.panicle"             
## [17] "Panicle.fertility"                "Seed.length"                     
## [19] "Seed.width"                       "Seed.volume"                     
## [21] "Seed.surface.area"                "Brown.rice.seed.length"          
## [23] "Brown.rice.seed.width"            "Brown.rice.surface.area"         
## [25] "Brown.rice.volume"                "Seed.length.width.ratio"         
## [27] "Brown.rice.length.width.ratio"    "Seed.color"                      
## [29] "Pericarp.color"                   "Straighthead.suseptability"      
## [31] "Blast.resistance"                 "Amylose.content"                 
## [33] "Alkali.spreading.value"           "Protein.content"                 
## [35] "Year07Flowering.time.at.Arkansas" "Year06Flowering.time.at.Arkansas"

両者に共通する個体を取り出す。

Common <- intersect( rownames(G_mat), rownames(Pheno) )
G_mat <- G_mat[ match( Common, rownames(G_mat)),]
Pheno <- Pheno[ match(Common, rownames(Pheno)),]
length(Common)
## [1] 388

マーカー行列のランク確認

今回はマーカー回帰をRidgeなどのスパースモデリングを用いずに直接行うが、そのためにはマーカーの行列のランクがマーカー数と一致している必要がある。そもそもこの仮定が成り立つかを確認する。行列のランクはqr関数を用いることで調べることができる。

d <- dim(G_mat)[2]
QR <- qr(G_mat)
print(c(d,QR$rank))
## [1] 38769   388

圧倒的にマーカー数のほうが多い。ここからは合理的な方法でマーカー数を削減し、行列のランクとマーカー数を一致させる。しかし、流石に数が多すぎるので、まずは偶数番のマーカーを抜くことにする。

G_mat <- G_mat[ , 1:ncol( G_mat ) %% 2 + 1 == 1 ]
dim(G_mat)
## [1]   388 19384

マーカーの削減

マーカーをどのように減らすかが問題になるが、似たマーカー同士はランク落ちを引き起こすため、どちらか片方を残すべきである。“似ている”をどう評価するかであるが、各マーカーをベクトルとみなし、その内積で考えるのが良いように思われる。以下では、各マーカーの長さを1にし、マーカー間の内積を行列で表した

G_scale <- apply( G_mat, 2, function( x ) sqrt( sum( x* x ) ) ) 
G_new <- sweep(G_mat, 2, G_scale, FUN = "/")
GRM <- t(G_new)%*%G_new
hist(GRM)

ヒストグラムでみると様々な分布を持っていることがわかる。 ではこの情報をもとに早速マーカー数を減らしていく。以下のコードでは、内積の絶対値が大きいマーカーのペアを見つけ、片方をなくす作業をランクとマーカー数が一致するまで繰り返している。

n <- 0.25
while( d != QR$rank ){
  
  Dupl <- which(GRM>=(1-n),arr.ind = TRUE)
  Dupl2 <- Dupl[which( Dupl[,1]!=Dupl[,2] ),]
  Dupl3 <- Dupl2[  ! which( Dupl2[,1]>Dupl2[,2])]
  G_mat2 <- G_mat[,!colnames(G_mat)%in%rownames(Dupl2)]
  
  X <- G_mat2
  X <- cbind( rep(1,nrow(X)), X)
  d <- dim(X)[2]
  QR <- qr(X)
  
  n <- n + 0.001
  print( 1 - n )
}
## [1] 0.749
## [1] 0.748
## [1] 0.747
## [1] 0.746
## [1] 0.745
## [1] 0.744
## [1] 0.743
## [1] 0.742
## [1] 0.741
## [1] 0.74
## [1] 0.739
## [1] 0.738
## [1] 0.737
## [1] 0.736
## [1] 0.735
## [1] 0.734
## [1] 0.733
## [1] 0.732
## [1] 0.731
## [1] 0.73
## [1] 0.729
## [1] 0.728
## [1] 0.727
## [1] 0.726
## [1] 0.725
## [1] 0.724
## [1] 0.723
## [1] 0.722
## [1] 0.721
## [1] 0.72
## [1] 0.719
## [1] 0.718
## [1] 0.717
## [1] 0.716
## [1] 0.715
## [1] 0.714
## [1] 0.713
## [1] 0.712
## [1] 0.711
print(d)
## [1] 375

最終的にマーカー数は350程度となったが、これは表現型の数とほぼ一致する。このことからも合理的な方法でマーカーを減らすことができたことがわかる。

マーカー回帰とその効果の分散

最後にマーカー回帰を行い、その効果の分散を求める。求め方はGBLUPの解を参考にしている。絶対値こそ小さいがやはりマーカーによってかなり分散が違うことがわかった。しかし、表現型の違いは誤差分散にしか影響を与えず、それはマーカー効果の誤差分散に対しては定数倍の違いしか生み出さないため、マーカー効果の分散は表現型間で定数倍を除き変化しなかった。

n_pheno <- ncol(Pheno)

for( j in 1:n_pheno ){
  y <- Pheno[,j]
  y[which(is.na(y))] <- mean(y,na.rm = TRUE)
  b <- solve( t(X) %*% X ) %*% t( X ) %*% matrix(y)
  barplot(b[-1,1],main=colnames(Pheno)[j], ylab= "marker effects")
}

for( i in 1:n_pheno ){
  y <- Pheno[,i]
  y[which(is.na(y))] <- mean(y,na.rm = TRUE)
  b <- solve( t(X) %*% X ) %*% t( X ) %*% matrix(y)
  e <- y - X %*% b
  vare <- (t(e)%*%e)/length(y)
  varb <- solve( t(X) %*% X ) * as.numeric(vare)
  barplot(diag(varb)[-1],main=colnames(Pheno)[i], ylab= "variance")
}