AlphaSimRとは
こちらの論文で紹介されている、植物および動物の育種プログラムをシミュレートするためのRパッケージです。Faux et al., 2016で紹介された’AlphaSim’の後継にあたり、Rの関数を実行することで育種をシミュレートします。遺伝的な集団構造の設定、交配や表現型の取得などが関数になっていてそれらを組み合わせて一連の育種業務をシミュレートすることができるようです。
今回は初めてこのパッケージを使ってみようと思っているのでこちらにある”Introduction to AlphaSimR”を読みながら時々立ち止まりつつ試してみたいと思います。
パッケージのロード
まず、下準備としてパッケージのロードと、乱数シードの設定を行います。
library(AlphaSimR) ##> Loading required package: R6 library(ggplot2) set.seed(utf8ToInt("🐶"))
先述したIndroductionのページでは、AlphaSimRで育種プログラムのシミュレートをする上での一般的な手順として、以下の4ステップがあげられています。今回の記事でもこの手順で進めつつ、各ステップで登場する関数とオブジェクトについて深掘りしながら見ていきたいと思います。
- Creating Founder Haplotypes
- Setting simulation
- Modeling the Breeding Program
- Examining the Results
(https://cran.r-project.org/web/packages/AlphaSimR/vignettes/intro.html より引用)
1. 創始者集団のハプロタイプを構築する (Create Founder Haplotypes)
最初に登場する関数はquickHaplo()です。こちらは創始者集団のハプロタイプを定義し、MapPopというオブジェクトとして返します。この関数は完全にランダムにハプロタイプを割り当てるので、より現実的なハプロタイプを得るにはrunMacs()やnewMapPop()関数の利用が推奨されます3。
# 10染色体を持ち、各染色体に1000カ所の均等な分離するサイト(locus)を持つ仮想的な集団1000個体を生成 founderPop <- quickHaplo(nInd = 1000, nChr = 10, segSites = 1000) class(founderPop) # オプジェクトのクラス名を見る ##> [1] "MapPop" ##> attr(,"package") ##> [1] "AlphaSimR" sloop::otype(founderPop) # オブジェクトのタイプを見る ##> [1] "S4" print(founderPop) ##> An object of class "MapPop" ##> Ploidy: 2 ##> Individuals: 1000 ##> Chromosomes: 10 ##> Loci: 10000
上記の例ではオプションに、染色体数として nChr = 10
、1染色体あたりの分離するサイト(座)数として segSites = 1000
を指定していました。そのため、printの出力では Loci: 10000 (= 10 *
1000) となっています。
MapPopオブジェクトに含まれる情報の中で重要なものの1つが遺伝地図情報かと思います。遺伝地図の情報を見てみるには
founderPop@genMap
とするか、getGenMap()関数を用います。以下では後者を利用して遺伝地図情報をデータフレームとして得た結果です。
head(getGenMap(founderPop)) ##> id chr pos ##> 1 1_1 1 0.000000000 ##> 2 1_2 1 0.001001001 ##> 3 1_3 1 0.002002002 ##> 4 1_4 1 0.003003003 ##> 5 1_5 1 0.004004004 ##> 6 1_6 1 0.005005005
表で見てもイメージしづらいので、getGenMap()で得られる情報から遺伝地図を簡易的にプロットしてみたいと思います。ただし、10000
locusもあると図が見づらくなってしまうので、遺伝地図情報に関係するオプションを適宜変えながら、図を見やすくするついでに、オプションを変えたときの結果として得られる遺伝地図の変化も見ていきます。
simple_genmap_plot <- function(mappop) { getGenMap(mappop) |> ggplot(aes(chr, pos)) + geom_line() + geom_point() + ggrepel::geom_text_repel(aes(label = id)) + scale_y_reverse() } # 染色体数を3、染色体毎サイト数を3に下げる simple_genmap_plot(quickHaplo(nInd = 1, nChr = 3, segSites = 3))
# 染色体数を3、サイト数はベクトルで指定することで染色体ごとに変えられる simple_genmap_plot(quickHaplo(nInd = 1, nChr = 3, segSites = c(3, 5, 7)))
# genLenオプションで染色体の長さを変えられる（segSites同様ベクトルでの指定も可） simple_genmap_plot(quickHaplo(nInd = 1, nChr = 3, segSites = 10, genLen = 2.0))
# runMacsを使うと現実に即したより複雑な遺伝地図が構築される simple_genmap_plot(runMacs(nInd = 1, nChr = 3, segSites = 10, species = "WHEAT"))
ここで、genMapのpos列の単位について個人的に少し気になったのですが、マニュアルにその言及がある箇所を見つけられませんでした。おそらくはM(モルガン)だろうなと思いますが、あとで組換え価を計算することによって実験的にこの予想の間接的な証拠を見てみたいと思います。
2. シミュレーションの設定 (Setting Simulation Parameters)
シミュレーションを実行していく上で、形質値がどんなQTLで定義されるのか、性別はセットするのかなど、最初に設定しておくべき項目を定義することができます。そのためにはまずSimParam$new()関数でSimParamオブジェクトを生成します。ここで、変数名をSP(大文字)としておくことには大きな理由があります。AlphaSimRで使われる多くの関数は内部でグローバル変数をSPという名前で検索し利用します。Introductionでも”The
name SP should almost always be used,
“と書いてあるので、ここは従います3。
ここで注意なのが、SPという名前の変数が存在する中で、SPを上書き変更しようと思って間違えて、たとえばspという小文字の別の変数に格納した場合は以降の関数にその変更の影響を与えることができなくなってしまうことです。しかし、エラーは起きないのでこのバグに気づきづらいという状態になります。個人的にはこういう状況があるのはプログラミング上よくないような気がしますが、1スクリプトで1シミュレーションにとどめるなどしたら良いのかもしれません。
SP <- SimParam$new(founderPop) class(SP) ##> [1] "SimParam" "R6" sloop::otype(SP) ##> [1] "R6"
SimParamオブジェクトの実体はR6オブジェクトであり、Pythonといった他の多くのプログラミング言語に採用されているオブジェクト指向プログラミングの仕組みと似たものになっています。つまり、SimParamオブジェクトにはメソッドが実装されており、そのメソッドを呼び出すことで設定を実行していきます。
以下では1染色体につき、1000個の相加効果を持つQTLを設定します。また、性別については、シミュレーションの過程で集団に追加されていく順番に、システマチックに性別を割り当てるように設定しています。
SP$addTraitA(nQtlPerChr=1000) SP$setSexes("yes_sys")
SimParamクラスに実装されているメソッドは ?SimParam
で見ることができます。非常に多くのメソッドがあり、さまざまな設定が可能であることがわかります。
3. 育種業務をモデル化する (Modeling the Breeding Program)
今回の設定はこれで全て終わったので、ここで具体的な育種の過程をコードで表現していく段階になります。まずは創始者ハプロタイプと設定したパラメータから、初期育成集団を作成することになります。newPop()関数を使うことで表現します。ここで、先ほど作成したSPがコードには全く出てきませんが、裏でSPという名前の変数を読み込んでいるため、先ほどの設定はちゃんと読み込まれていることに注意です。
pop <- newPop(founderPop) class(pop) ##> [1] "Pop" ##> attr(,"package") ##> [1] "AlphaSimR" print(pop) ##> An object of class "Pop" ##> Ploidy: 2 ##> Individuals: 1000 ##> Chromosomes: 10 ##> Loci: 10000 ##> Traits: 1
newPop()関数で得られるオブジェクトはPopクラスです。シミュレーションの過程ではこのクラスのオブジェクトに対して各種操作(選抜や交配)を適用していくことになります。print(pop)の出力を見てみると、先ほどのMapPopオブジェクトの出力で見られていた情報に加え、Traitsという項目が増えていることがわかり、
先ほどのSP$addTraitA(nQtlPerChr=1000)
というコードがちゃんと機能していることがわかります。
実際にこのTraitsに関するGenetic
valueをgv()関数により取得し分布を見ると以下のようになります。
firstPopGv <- gv(pop) hist(firstPopGv, probability = TRUE) # 正規分布の密度関数も一緒にプロット curve(dnorm(x), from = min(firstPopGv), to = max(firstPopGv), add = TRUE)
SP$addTraitA()はデフォルトで平均0分散1の正規分布に従うような相加効果を与えるらしいです。
ちなみに、Additive effectによるGenetic
valueがどのように計算されているかをあえて手動で計算すると以下のようになるかと思います。遺伝子型はpullQtlGeno()で、各QTLの効果はSimParamオブジェクトのtraitsプロパティからアクセスできます。
firstPopQtlGeno <- pullQtlGeno(pop) unique(as.vector(firstPopQtlGeno)) ##> [1] 0 1 2 summary(SP$traits[[1]]@addEff) ##> Min. 1st Qu. Median Mean 3rd Qu. Max. ##> -5.272e-02 -9.382e-03 -2.498e-05 -2.858e-05 9.293e-03 4.918e-02 SP$traits[[1]]@intercept ##> [1] 0.02822522 # Additive EffectによるGenetic valueを遺伝子型と書くlocusの効果から手動で計算 firstPopGvManual <- (firstPopQtlGeno - 1) %*% SP$traits[[1]]@addEff + SP$traits[[1]]@intercept # gv()関数の結果と一致することを確認 all.equal(firstPopGv, firstPopGvManual, check.attributes = FALSE) ##> [1] TRUE
これをみると遺伝子型は0, 1, 2でコードされており、0のときはその遺伝子座の効果を引き、1のときは何もせず(中間の効果)、そして2のときは効果を足し合わせているようなイメージになっているかと思います。
それでは、実際に育種プログラム全体をコードで表現してみようと思います。以下のコードでは、狭義の遺伝率h2 = 0.2として表現型をとり(setPheno()関数)、20回の世代更新を行なっています。各世代では表現型値で上位のメス500個体とオス25個体を選抜してランダムに1000組合せの交配を行い(selectCross()関数)、後代を得ています。また、各世代での平均表現型値を取得して(meanP関数)、phenoMean変数に記録しています。
pop <- setPheno(pop, h2 = 0.2) phenoMean <- meanP(pop) for(generation in 1:20){ pop <- setPheno(pop, h2 = 0.2) pop <- selectCross(pop=pop, nFemale=500, nMale=25, use="pheno", nCrosses=1000) pop <- setPheno(pop, h2 = 0.2) phenoMean <- c(phenoMean, meanP(pop)) }
4. 結果を評価する (Examining the Results)
最後に先ほどの20世代の育種プログラムの結果を、世代を通しての平均表現型値の変遷の折れ線グラフとしてプロットしてみます。
plot(0:20, phenoMean, xlab="Generation", ylab="Mean Phenotype Value", type="l")
各世代で順調に平均表現型値が上がっていて育種の効果があることを示しています。実際の育種現場ではさらに何個もの形質を評価しながら進めていかなければならないのでこんな単純なケースはそうそう無いとは思いますが、ひとまず、単純な育種プロセスをシミュレートすることができました。
Breeder’s equation
Breeder’s equationをシミュレーションから見てみたいと思います。1世代経たときの親集団と子集団の形質平均値の差をR(引用記事中ではΔZ)、選択された集団と元の集団の形質平均値の差S(selection differential)、狭義の遺伝率をh2とするとBreeder’s equationは以下の様になります。
R = h2 * S
h2はsetPheno()関数から設定可能なので、いくつかのh2についてシミュレートしてみたいと思います。
親集団1000個体から形質値で上位10%を選抜して、ランダムに交配し後代についても1000個体得るということを5回試行し、結果を表示してみます。結果のプロットではh2をx軸にとり、R/Sの値をy軸にとって散布図をかきました。Breeder’s equationから、このプロットはy=xの直線上付近に見られると予想されます。
examine_breeders_equation <- function(h2) { pop <- newPop(founderPop) pop <- setPheno(pop, h2 = h2) n <- nInd(pop) popSel <- selectInd(pop, nInd = n %/% 10, use = "pheno") pop2 <- randCross(popSel, nCrosses = n) pop2 <- setPheno(pop2, h2 = h2) c( R = unname(meanP(pop2) - meanP(pop)), S = unname(meanP(popSel) - meanP(pop)), h2 = h2 ) } breeders_equation_result <- purrr::map( seq(0.1, 1.0, 0.2), ~as.data.frame(t(replicate(5, examine_breeders_equation(.x)))) ) |> dplyr::bind_rows() breeders_equation_result |> ggplot(aes(h2, R/S)) + geom_point() + geom_abline(slope = 1, intercept = 0) + coord_equal()
実際y=x付近にプロットされました。
世代を超えたIBDの追跡（組換えの可視化）
AlphaSimRではSP$setTrackRec(TRUE)とすることで組換えイベントを保持してシミュレーションしてくれます。これを行うと、pullIbdHaplo()関数で集団のハプロタイプをIBD(Identidy
by
descent)ベースで見ることができます。これは創始者集団の染色体断片がどのように後代に遺伝するかを見ることができるということです。そこで、簡単なグフラフィカルジェノタイプ的なプロットをする関数を作って実際に集団のいくつかの個体についてのIBDハプロタイプを見てみたいと思います。
描画を簡単にするため創始者集団は2個体の二倍体で、1染色体のみを持つとします。すると創始者集団にはハプロタイプとして合計4種類を含むことになり、これらを4色で色分けしてプロットします。また、データ数が多くなるとグラフが見にくくなるので、各世代10個体のみを得ることにしました。
simple_ibd_haplotype_plot <- function(pop, map) { pullIbdHaplo(pop) |> as.data.frame() |> tibble::rownames_to_column("haplo_id") |> tidyr::pivot_longer(-haplo_id, values_to = "haplotype") |> dplyr::mutate( haplotype = factor(haplotype) ) |> dplyr::full_join(map, by = c("name" = "id")) |> ggplot(aes(pos, haplo_id)) + geom_raster(aes(fill = haplotype)) + theme(axis.text.x = element_blank()) } founderPop <- quickHaplo(2, nChr = 1, segSites = 100) map <- getGenMap(founderPop) SP <- SimParam$new(founderPop) SP$setTrackRec(TRUE) pop <- newPop(founderPop) simple_ibd_haplotype_plot(pop, map) # 最初の集団
それでは、上記4つに色分けされたハプロタイプ（染色体）がどう後代に遺伝するかを追ってみます。
S1世代
組換えが起こり、色分けにスイッチングしている部分が見られます。二倍体が10個体なので、合計20本分のデータがプロットに見られます。
pop <- makeCross(pop, crossPlan = matrix(c(1, 2), ncol = 2), nProgeny = 10) simple_ibd_haplotype_plot(pop, map) # 第1世代
S2世代
pop <- randCross(pop, nCrosses = 10) simple_ibd_haplotype_plot(pop, map) # 第2世代
S10世代
世代を経るごとに連鎖ブロックが途切れ途切れになっている様子が見えるかと思います。また、集団サイズを10と小さく取っているからか2番目のハプロタイプがかなり除かれていることが見て取れます。
for (i in 1:8) { pop <- randCross(pop, nCrosses = 10) } simple_ibd_haplotype_plot(pop, map) # 第10世代
遺伝距離と組換え価の関係
先ほど、getGenMap()で得られるpos列の単位について、おそらくM(モルガン)であろうという話をしましたが、それを間接的に裏付けるデータを取ってみようと思います。具体的には親を1個体を用意して、その自殖第1世代(S1)1000個体を得るということをシミュレーションしてみます。このとき、親、およびS1のハプロタイプはpullSegSiteHaplo()で完全なものが得られるため、簡単に組換え価が計算できます。具体的な計算はコードを見ていただくとして、実際に組換え価とposの値の差の関係をプロットしてみます。
founderPop <- runMacs(nInd = 1, nChr = 1, segSites = 20, species = "MAIZE") pos <- getGenMap(founderPop)$pos SP <- SimParam$new(founderPop) # indexが1の個体を親として新しいPopオブジェクト作成 parent <- newPop(founderPop[1]) haplo_p <- pullSegSiteHaplo(parent) haplo_p_sub <- haplo_p[, haplo_p[1, ] != haplo_p[2, ]] # 後代で分離するサイトのみ保持 # 親個体を自殖して1000個体を得る s1 <- self(parent, nProgeny = 1000) haplo_s1_sub <- pullSegSiteHaplo(s1)[, haplo_p[1, ] != haplo_p[2, ]] # 各サイト間の遺伝距離(pos列の差分)と組換え価を計算する n <- ncol(haplo_p_sub) recomb_rate <- NULL pos_diff <- NULL for (i in 1:(n-1)) { for (j in (i+1):n) { match_to_haplo1 <- apply(haplo_s1_sub[, c(i, j)], 1, function(x) all(x == haplo_p_sub[1, c(i, j)])) match_to_haplo2 <- apply(haplo_s1_sub[, c(i, j)], 1, function(x) all(x == haplo_p_sub[2, c(i, j)])) recomb_rate <- c(recomb_rate, mean(!(match_to_haplo1 | match_to_haplo2))) pos_diff <- c(pos_diff, abs(pos[j] - pos[i])) } } plot(pos_diff, recomb_rate) curve((1-exp(-2 * x))/2, from = min(pos_diff), to = max(pos_diff), add = TRUE, col = "blue") curve(tanh(2 * x)/2, from = min(pos_diff), to = max(pos_diff), add = TRUE, col = "red") legend("bottomright", c("Haldane", "Kosambi"), col = c("blue", "red"), lty =c(1, 1))
するとこれらの関係は綺麗にKosambi mapping function上に乗ってきます。このことから、おそらくposはM単位の値であると考えられるかと思います。
あえてシミュレーションパラメータにSPという名前をつけない例
最後に、AlphaSimRでは独特の作法でSimParamオブジェクトにはSPという変数名をつけるということが推奨されています。ただ、これは一般的なプログラミングにおいては普通ではないので、なんとなく違和感を感じるという方のために、SPという変数名を使わない場合のコードを以下に置いておきます。
具体的には各関数にあるsimParam引数にSimParamオブジェクトを格納した変数を与えてやることで可能になります。以下の例では大文字のSP（推奨の変数名）の代わりに小文字のspという変数名を用いてみています。ちなみにこちらの例をもとにしています。
rm(list = ls()) # 先ほど定義したSPとその他もろもろをを消しておく founderPop <- quickHaplo(nInd=1000, nChr=10, segSites=1000) sp <- SimParam$new(founderPop) sp$addTraitA(nQtlPerChr=1000) sp$setSexes("yes_sys") #pop <- newPop(founderPop) # SPという変数がない場合、これではerrorが発生 pop <- newPop(founderPop, simParam = sp) genMean <- meanG(pop) for(generation in 1:20){ pop = selectCross(pop=pop, nFemale=500, nMale=25, use="gv", nCrosses=1000, simParam = sp) genMean <- c(genMean, meanG(pop)) } # Examining the Results plot(0:20, genMean, xlab="Generation", ylab="Mean Genetic Value", type="l")
参考文献
- Gaynor, R. Chris, Gregor Gorjanc, and John M. Hickey. 2021. AlphaSimR: an R package for breeding program simulations. G3 Gene|Genomes|Genetics 11(2):jkaa017. https://doi.org/10.1093/g3journal/jkaa017.
- Faux, A.-M., Gorjanc, G., Gaynor, R.C., Battagin, M., Edwards, S.M., Wilson, D.L., Hearne, S.J., Gonen, S. and Hickey, J.M. (2016), AlphaSim: Software for Breeding Program Simulation. The Plant Genome, 9: plantgenome2016.02.0013. https://doi.org/10.3835/plantgenome2016.02.0013
- http://cran.r-project.org/web/packages/AlphaSimR/vignettes/intro.html
- https://www.nature.com/scitable/knowledge/library/the-breeder-s-equation-24204828/
- https://en.wikipedia.org/wiki/Genetic_map_function
