phasedされたゲノムデータから組み換えを考慮した上で後代のジェノタイプを生成する方法を軽くここにまとめる。なお、この方法はbreedSimulatRに実装されている方法を参考にしている。

配偶子の生成まで

まずはゲノムデータから組み換えを考慮して配偶子を生成する方法をまとめる。この配偶子からさらに後代の遺伝子型まで生成する方法はこの次にまとめる。

ゲノムデータの準備

まずはゲノム情報を人工的に生成する。ここではマーカーは1000個とし、Mというリストにphasedされたマーカー情報と遺伝子型をを格納してある。

m <- 1000 # number of markers
ChrName <- paste( rep( "Chr1", m), 1:m, sep = "_") # name of each marker

M <- vector( "list", 2 )
names( M ) <- c( "gamete", "genotype" )

# randomly generate gamete deriving from male and female 
M$gamete <- matrix( sample( 0 : 1, 2* m, replace = TRUE ), ncol = m ) 

colnames( M$gamete ) <- ChrName
rownames( M$gamete ) <- c( "male", "female" )

M$gamete[,1:5]
##        Chr1_1 Chr1_2 Chr1_3 Chr1_4 Chr1_5
## male        1      0      1      1      0
## female      1      0      0      0      0
# calculate genotype from gamete
M$genotype <- apply( M$gamete, 2, sum )
names( M$genotype ) <- ChrName

head( M$genotype )
## Chr1_1 Chr1_2 Chr1_3 Chr1_4 Chr1_5 Chr1_6 
##      2      0      1      1      0      1

次にマーカーの情報を格納したM_infというデータフレームを生成する。ここにはマーカーの名前とその物理位置、染色体地図上の位置の三つの情報が格納されている。図は物理位置と染色体地図の位置を比較したものである。本来は非線形の関係が見られるが、ここでは人工的に生成しているためそのような関係は見られない。なお、染色体の長さは1e6, cMでの長さは100としてある。

# marker information
LnChr <- 1e6
LnChrMorg <- 100

# sample markers based on the physical position
PhyPos <- sort( runif( m, 0, LnChr ) )

# sample markers based on the gene map
MorgPos <- sort( runif( m, 0, LnChrMorg ) )

# create data frame to store genome information
M_inf <- data.frame( name = ChrName, Pos = PhyPos, Pos_Chr = MorgPos )

plot( M_inf$Pos, M_inf$Pos_Chr )

組み換えと配偶子の生成

次に上で準備したマーカー情報とゲノムデータをもとに配偶子を生成する。生成する際はまずポワソン分布に基づいて組み換えの数を指定し、その数分染色体地図場から組み換えの位置をサンプルした。そのうえでその組み換え位置ごとに配偶子を切り替えこの配偶子を生成した。

# randomly take a gamete from tow gametes
i_gamete <- sample( 1:2 ,1 )
g <- M$gamete[ (i_gamete %% 2 + 1 ), ]

# calculate the number of cross over based on poisson distribution
nCO <- rpois( 1, LnChrMorg )

# sample the position of cross over based on gene map
posCO <- sort( runif( nCO, 0, LnChrMorg ) )
posCO <- c( 0, posCO, LnChrMorg )


g_ind <- rep( NA, m )
for( i in 1 : ( nCO + 1 ) ){
  g_ind[ posCO[ i ] < M_inf$Pos_Chr & M_inf$Pos_Chr < posCO[ (i + 1 ) ] ] <- 
    g[ posCO[ i ] < M_inf$Pos_Chr & M_inf$Pos_Chr < posCO[ (i + 1 ) ] ]
  i_gamete <- i_gamete + 1
  g <- M$gamete[ (i_gamete %% 2 + 1 ), ]
}

names( g_ind ) <- ChrName
head( g_ind )
## Chr1_1 Chr1_2 Chr1_3 Chr1_4 Chr1_5 Chr1_6 
##      1      0      0      0      0      0

関数化

後代の遺伝子型の生成に向けここではマーカー情報、ゲノムデータから配偶子を生成するまでを関数化する。

CreateGamete <- function( M, M_inf, LnChrMorg ){
  # randomly take a gamete from tow gametes
  i_gamete <- sample( 1:2 ,1 )
  g <- M$gamete[ (i_gamete %% 2 + 1 ), ]

  # calculate the number of cross over based on poisson distribution
  nCO <- rpois( 1, LnChrMorg )

  # sample the position of cross over based on gene map
  posCO <- sort( runif( nCO, 0, LnChrMorg ) )
  posCO <- c( 0, posCO, LnChrMorg )


  g_ind <- rep( NA, m )
  for( i in 1 : ( nCO + 1 ) ){
    g_ind[ posCO[ i ] < M_inf$Pos_Chr & M_inf$Pos_Chr < posCO[ (i + 1 ) ] ] <- 
      g[ posCO[ i ] < M_inf$Pos_Chr & M_inf$Pos_Chr < posCO[ (i + 1 ) ] ]
    i_gamete <- i_gamete + 1
    g <- M$gamete[ (i_gamete %% 2 + 1 ), ]
  }

  names( g_ind ) <- M_inf$name

  return( g_ind )

}

これを動かすと確かに配偶子が生成されていることがわかる。

head( CreateGamete( M, M_inf, LnChrMorg) )
## Chr1_1 Chr1_2 Chr1_3 Chr1_4 Chr1_5 Chr1_6 
##      1      0      1      1      0      1

後代の遺伝子型の生成

次に母親・父親両方を準備し、それぞれから配偶子を生成し、そのうえで後代の遺伝子型までの生成を試みる。

交配親の情報を生成

まずは両親のゲノムデータを生成する。ここではM1とM2に格納した

M1 <- vector( "list", 2 )
M2 <- vector( "list", 2 )
names( M1 ) <- c( "gamete", "genotype" )
names( M2 ) <- c( "gamete", "genotype" )

# randomly generate gamete deriving from male and female 
M1$gamete <- matrix( sample( 0 : 1, 2* m, replace = TRUE ), ncol = m ) 
M2$gamete <- matrix( sample( 0 : 1, 2* m, replace = TRUE ), ncol = m ) 

colnames( M1$gamete ) <- ChrName
colnames( M2$gamete ) <- ChrName
rownames( M1$gamete ) <- c( "male", "female" )
rownames( M2$gamete ) <- c( "male", "female" )

# calculate genotype from gamete
M1$genotype <- apply( M1$gamete, 2, sum )
M2$genotype <- apply( M2$gamete, 2, sum )
names( M1$genotype ) <- ChrName
names( M2$genotype ) <- ChrName

配偶子の生成と合成

次にそれぞれから配偶子を生成し、それを統合して後代の遺伝子型を決定する。決定した子のゲノムデータはM3に格納した。

gamete1 <- CreateGamete( M1, M_inf, LnChrMorg)
gamete2 <- CreateGamete( M2, M_inf, LnChrMorg)

M3 <- vector( "list", 2 )
names( M3 ) <- c( "gamete", "genotype" )

M3$gamete <- t( cbind( gamete1, gamete2 ) )
M3$genotype <- apply( M3$gamete, 2, sum )
names( M3$genotype ) <- ChrName

関数化

CreateGenotype <- function( M1, M2, M_inf, LnChrMorg ){
  
  gamete1 <- CreateGamete( M1, M_inf, LnChrMorg)
  gamete2 <- CreateGamete( M2, M_inf, LnChrMorg)

  M3 <- vector( "list", 2 )
  names( M3 ) <- c( "gamete", "genotype" )

  M3$gamete <- t( cbind( gamete1, gamete2 ) )
  M3$genotype <- apply( M3$gamete, 2, sum )
  names( M3$genotype ) <- ChrName
  
  return( M3 )
  
}

これで自由に子供の遺伝子型を生成することができるようになった

ind1 <- CreateGenotype( M1, M2, M_inf, LnChrMorg )
ind2 <- CreateGenotype( M1, M2, M_inf, LnChrMorg )
ind3 <- CreateGenotype( M1, M2, M_inf, LnChrMorg )
ind4 <- CreateGenotype( M1, M2, M_inf, LnChrMorg )