Select500Accs

SNPデータの読み込みにgastonを用いる。

require(gaston) # vcfファイルの読み込み、GRMの計算 
require(cluster) # pam（k-medoids）の実行

データの読み込み。

# vcfファイルの読み込み
bm <- read.vcf("data/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.

bm # データのsummaryが表示される

## A bed.matrix with 1568 individuals and 38769 markers.
## snps stats are set
## ped stats are set

標準化された内積行列（ゲノム関係行列）の計算

standardize(bm) <- "mu_sigma"
B <- GRM(bm) * (ncol(bm) - 1)

関係行列からの距離を計算 \[b_{ij} = -\frac{1}{2}(d_{ij}^2 - d_{i.}^2 - d_{.j}^2 + d_{..}^2)\] \[-2 b_{ij} = d_{ij}^2 - d_{i.}^2 - d_{.j}^2 + d_{..}^2\] \[d_{ij}^2 = -2 b_{ij} + d_{i.}^2 + d_{.j}^2 - d_{..}^2\]

n <- nrow(B)
sum_bii.div.n <- sum(diag(B)) / n
di.2 <- diag(B) + sum_bii.div.n
d..2 <- 2 * sum_bii.div.n
D2 = - 2 * B + matrix(rep(di.2, n), nrow = n, byrow = T) + 
  matrix(rep(di.2, n), nrow = n, byrow = F) - d..2
D2[D2 < 0] <- 0

k-medoidsを実行して、500の代表を選ぶ

k <- 500
res <- pam(sqrt(D2), k, diss = T)

df <- data.frame(id.med = sort(res$id.med))
write.csv(df, "data/id_med.csv", quote = F)