2011-05-07

家系を扱おう

R DNA鑑定

家系情報を作る
- 家系情報は行列にする
- 行の順番は時系列を守ること
- 後述するように、尤度・確率の計算では、核家族を分離して計算することで計算量を減らすように作っている。したがって、家系図にループがある場合には、後半が対応していないことに留意すること
- 第１カラムは個人ID(上から順に1,2,...の連番)
- 第２カラムは母の個人ID。母が家系に入らないときは0
- 第３カラムは父。0は同様の定義
- 第４カラムは性別は男が0、女が1。性別不明は(今は不可)
- 第５カラムは、ディプロタイプ情報が得られるかどうかの情報。1は得られる。3は得られない。2は得られなかったり、得られたり。2は、その個人かどうかが不明な候補者のディプロタイプがある場合を想定している

p<-matrix(
c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
  0, 0, 0, 0, 2, 2, 4, 4, 6,  6,  0,  0, 12, 13,
  0, 0, 0, 0, 1, 1, 3, 3, 7,  7,  0,  0, 11, 10,
  0, 1, 0, 1, 0, 1, 0, 1, 0,  0,  0,  1,  1,  1,
  3,1,1,1,3,1,1,1,1,2,1,1,2,1),
  ncol=5)

> p
      [,1] [,2] [,3] [,4] [,5]
 [1,]    1    0    0    0    3
 [2,]    2    0    0    1    1
 [3,]    3    0    0    0    1
 [4,]    4    0    0    1    1
 [5,]    5    2    1    0    3
 [6,]    6    2    1    1    1
 [7,]    7    4    3    0    1
 [8,]    8    4    3    1    1
 [9,]    9    6    7    0    1
[10,]   10    6    7    0    2
[11,]   11    0    0    0    1
[12,]   12    0    0    1    1
[13,]   13   12   11    1    2
[14,]   14   13   10    1    1

家系図を描く→掲載図
- ディプロタイプが得られるかどうかの情報(第５カラムの1,2,3)を描き分ける
- 黒は1(ディプロタイプが確定的に得られる)
- 白は2か3
- 白は斜線があるかどうかでわかれる
  - 斜線ありは、3、斜線なしは2。斜線なしの白が探索対象
- 家系情報行列からkinshipパッケージのpedigreeオブジェクトを作って描図する

MakePedigreeFromFamilyInfo<-function(p){
	ns<-length(p[,1])
	affected<-status<-rep(1,ns)
	affected[which(p[,5]==2)]<-0
	affected[which(p[,5]==3)]<-0
	status[which(p[,5]==1)]<-0
	status[which(p[,5]==2)]<-0
	ptemp<-pedigree(id=p[,1],dadid=p[,3],momid=p[,2],sex=p[,4],affected=affected,status=status)
	if(sum(ptemp$affected)==0)ptemp$affected<-affected
	ptemp
}

library(kinship)
ptemp<-MakePedigreeFromFamilyInfo(p)
plot(ptemp)

家系図から核家族を抜き出す
- 核家族は、父母と１人以上の子からなる
- 核家族は、家系図と１人以上のメンバーで連結する
  - 特に、ただ一人のメンバーで家系図に連結してぶら下がっている核家族は「家系図の辺縁の核家族」である
    - この連結メンバーをピボットと呼ぶ
  - それ以外の核家族は家系図の中で二つ以上の核家族に連結している
  - 辺縁から核家族を選び進め、辺縁の核家族のピボットメンバー以外をこの核家族のピボットとすることで、すべての核家族は異なるメンバーをピボットとして持ち、すべての核家族には０ないし１人のピボットメンバーを持たせることができる。連結した家系図の核家族の中で、ピボットメンバーが０の核家族はただ１つ存在する
- paramlinkパッケージのlinkdat()関数内に書かれた、関数内関数を借用する

subnucs <- function(ped) {# subfunction in linkdat() in package "paramlink"
        parents <- unique(ped[, 2:3])
        parents = parents[-match(0, parents[, 1]), , drop = FALSE]
        list1 <- lapply(nrow(parents):1, function(i) {
            par = parents[i, ]
            c(fa = par[[1]], mo = par[[2]], offs = as.vector(ped[, 
                1])[which(ped[, 2] == par[[1]] & ped[, 3] == 
                par[[2]], useNames = FALSE)])
        })
        res = list()
        i = 1
        k = 1
        while (length(list1) > 1) {
            if (i > length(list1)) {
                warning("Loop detected, likelihood calculations will not work.")
                return(NULL)
            }
            link = ((sub <- list1[[i]]) %in% unlist(list1[-i]))
            if (sum(link) == 1) {
                res[[k]] <- c(pivot = sub[[which(link)]], sub)
                list1 <- list1[-i]
                k <- k + 1
                i <- 1
            }
            else i <- i + 1
        }
        res[[k]] <- c(pivot = 0, list1[[1]])
        res
    }

nucs<-subnucs(p)

> nucs
[[1]]
pivot    fa    mo  offs 
   13    12    11    13 

[[2]]
pivot    fa    mo  offs 
   10    13    10    14 

[[3]]
pivot    fa    mo offs1 offs2 
    7     4     3     7     8 

[[4]]
pivot    fa    mo offs1 offs2 
    6     6     7     9    10 

[[5]]
pivot    fa    mo offs1 offs2 
    0     2     1     5     6