クラス・メソッド・総称的関数

クラス (class) 構造をもつオブジェクトで,クラス属性 (class attribute) という文字列がひとつまたは複数登録されている.Rの関数が返すオブジェクトはたいてい何らかのクラスに属し,特定の構造をもっている.また,特定のクラスオブジェクトを引数にとれる関数も数多くある.
異なるクラスオブジェクトに対する類似の処理を,ひとつの関数で代表させる機構が,総称的関数 (generic function) とメソッド選択適用 (method dispatch).総称的関数"function()"に,引数としてクラス"class"が与えられると,"function.class()"という名前の関数が探され適用される.

> # plot関数のメソッド関数一覧
> methods(plot)
 [1] plot.acf*           plot.data.frame*    plot.decomposed.ts*
 [4] plot.default        plot.dendrogram*    plot.density       
 [7] plot.ecdf           plot.factor*        plot.formula*      
[10] plot.function       plot.hclust*        plot.histogram*    
[13] plot.HoltWinters*   plot.isoreg*        plot.lm            
[16] plot.medpolish*     plot.mlm            plot.ppr*          
[19] plot.prcomp*        plot.princomp*      plot.profile.nls*  
[22] plot.spec           plot.stepfun        plot.stl*          
[25] plot.table*         plot.ts             plot.tskernel*     
[28] plot.TukeyHSD      

   Non-visible functions are asterisked
> 


S3クラス・メソッド

本質的には,クラス属性文字列の与えられたオブジェクトにすぎない.該当クラスにふさわしい構造をもつかどうかに関係なく与えることができてしまう.

> x <- runif(10)
> class(x)
[1] "numeric"
> class(x) <- "lm"
> class(x) # だませてしまう
[1] "lm"
> summary(x) # しかし無意味
 エラー: $ operator is invalid for atomic vectors
> 
y <- c("R", "language")
class(y)
class(y) <- "matrix" # 多少の整合性チェックはする

Rシステムに組み込みのS3総称的関数は以下で確認できる.

> names(.knownS3Generics)
 [1] "Math"           "Ops"            "Summary"        "Complex"       
 [5] "as.character"   "as.data.frame"  "as.environment" "as.matrix"     
 [9] "as.vector"      "cbind"          "labels"         "print"         
[13] "rbind"          "rep"            "seq"            "seq.int"       
[17] "solve"          "summary"        "t"              "edit"          
[21] "str"            "contour"        "hist"           "identify"      
[25] "image"          "lines"          "pairs"          "plot"          
[29] "points"         "text"           "add1"           "AIC"           
[33] "anova"          "biplot"         "coef"           "confint"       
[37] "deviance"       "df.residual"    "drop1"          "extractAIC"    
[41] "fitted"         "formula"        "logLik"         "model.frame"   
[45] "model.matrix"   "predict"        "profile"        "qqnorm"        
[49] "residuals"      "se.contrast"    "terms"          "update"        
[53] "vcov"          
> 


S4クラス・メソッド

含むべきデータであるスロット (slot) とその型 (type) が明示的に表現され,チェックがはたらくようになっている.また,他のクラスを継承 (inheritance) したり,クラス定義を変更されないように封印 (seal) することができる.クラスオブジェクトには原型 (prototype) を与えておくことができる.クラス定義はパッケージ (環境) に記録される.

S4クラスの定義

setClass()関数で定義,名前つきスロットはrepresentationに,名前なしスロットはcontainsに書く.そのクラス属性を持つ新しいオブジェクト (インスタンス) はnew()関数で生成する.

> setClass("ex", representation(id = "character"), contains = "numeric")
> getClass("ex")
Class "ex" [in ".GlobalEnv"]

Slots:
                          
Name:      .Data        id
Class:   numeric character

Extends: 
Class "numeric", from data part
Class "vector", by class "numeric", distance 2
> # 名前なしスロットは .Data として与えられる
> setClass("mat", representation(id = "character", num = "matrix"))
> getClass("mat")
Class "mat" [in ".GlobalEnv"]

Slots:
                          
Name:         id       num
Class: character    matrix
> 

クラス"ex"のインスタンスを作成

> (ex.1 <- new("ex", runif(3), id = "first 3"))
An object of class "ex"
[1] 0.3433830 0.8273612 0.7917083
Slot "id":
[1] "first 3"

> str(ex.1)
Formal class 'ex' [package ".GlobalEnv"] with 2 slots
  ..@ .Data: num [1:3] 0.343 0.827 0.792
  ..@ id   : chr "first 3"
> 

不適切なオブジェクトはスロットに代入できない.

> ex.2 <- new("ex", runif(3), id = numeric(2))
 以下にエラー validObject(.Object) : 
   不正なクラス “ex” オブジェクト:  invalid object for slot "id" in class "ex": got class "numeric", should be or extend class "character"
> 

クラススロットの操作

スロットの中身を抜き出すには @ 演算子を使う.または,slot(instance, ".Data") などのようにする.

> ex.1@id
[1] "first 3"
> ex.1@.Data
[1] 0.3433830 0.8273612 0.7917083
> slot(ex.1, "id")
[1] "first 3"
> 

スロットにオブジェクトを付値.整合性の検査機構がはたらく.

(> (ex.1@id <- "First 3")
[1] "First 3"
> ex.2 <- new("ex", representation(id = numeric(2), contains = 1:3))
 以下にエラー representation(id = numeric(2), contains = 1:3) : 
   representation の要素 1 は単一の文字列ではありません 
> 

クラスの拡張

containsは,すでにあるクラスから拡張 (etend) されたクラスと定義される.
既存のクラスは上位クラス (super-class),拡張されたクラスは下位クラス (sub-class) となる.

スロットのプロトタイプ

スロットにはプロトタイプ (既定値) を与えられる.

> setClass("prot", representation(id = "character", num = "numeric"),
+   prototype = list(id = "ND", num = 0))
> (prot.1 <- new("prot", id = "Carbon", num = 12))
An object of class "prot"
Slot "id":
[1] "Carbon"

Slot "num":
[1] 12

> (prot.2 <- new("prot"))
An object of class "prot"
Slot "id":
[1] "ND"

Slot "num":
[1] 0

> 

クラスの継承

既存のクラスをクラス定義に含めることができる.含まれたクラスは上位クラス,含むクラスは下位 (拡張) クラスとなる.

> setClass("prot2", representation(id2 = "character", num2 = "numeric"),
+   prototype = list(id2 = "YES", num2 = 1), contains = "prot")
> (prot2.1 <- new("prot2", id2 = "Nitrogen", num2 = 14, prot.1))
An object of class "prot2"
Slot "id2":
[1] "Nitrogen"

Slot "num2":
[1] 14

Slot "id":
[1] "Carbon"

Slot "num":
[1] 12

>

同じ名前のスロットがあると,同じものとされてしまう.

> setClass("prot3", representation(id = "character", num3 = "numeric"),
+   prototype = list(id = "YES", num3 = 1), contains = "prot")
> (prot3.1 <- new("prot3", id = "Oxigen", num3 = 16, prot.1))
An object of class "prot3"
Slot "id":
[1] "Oxigen"

Slot "num3":
[1] 16

Slot "num":
[1] 12

> 

クラス定義の封印

クラス定義を変更できないようにする.

> setClass("ex3", representation(x = "numeric", y = "numeric"))
> setClass("ex3", representation(x = "logical", z = "numeric"))
> sealClass("ex3", where = .GlobalEnv)
> setClass("ex3", representation(x = "numeric", y = "numeric"))
 以下にエラー setClass("ex3", representation(x = "numeric", y = "numeric")) : 
   "ex3" は封印されたクラス定義を持ち、再定義は出来ません 
> 

S4メソッド

総称的関数を指定して,それに対するメソッド関数を定義する.総称的関数は,引数オブジェクトのクラスの組み合わせに対応する最適なメソッド関数を探して,割り当てる.

> # クラス"track"を定義
> setClass("track", representation(x = "numeric", y = "numeric"))
> 
> # 組み込みの算術演算総称的関数Arithに対し,"track"クラスに対応するメソッドを定義
> setMethod("Arith", c("track", "numeric"),
+   function(e1, e2){e1@y <- callGeneric(e1@y, e2); e1})
[1] "Arith"
> # 逆も
> setMethod("Arith", c("numeric", "track"),
+   function(e1, e2){e2@y <- callGeneric(e1, e2@y); e2})
[1] "Arith"
> 
> # "track"クラスをもつインスタンス生成
> t1 <- new("track", x = 1:3, y = rnorm(1))
> t1 - 1 # 引き算メソッドが適用される
An object of class "track"
Slot "x":
[1] 1 2 3

Slot "y":
[1] -3.623696

> 1 / t1 # 割り算メソッドが適用される
An object of class "track"
Slot "x":
[1] 1 2 3

Slot "y":
[1] -0.3811417

> 

言語オブジェクト

Rの言語オブジェクトについてまとめてみました.
Rプログラミングマニュアル』(間瀬茂, 2007年, 数理工学社) を参考にしました.
言語としてのRの基礎や,豊富な役立つtipsが載っており,とても参考になる本です.


言語オブジェクト

Rの言語オブジェクトには3つの型がある.呼び出し call,表現式 expression,名前 name である.


呼び出しオブジェクトは関数名と引数を与えてつくる.
作成したオブジェクトは,eval()関数で評価実行できる.

> c1 <- call("round", 2.4)
> c1
round(2.4)
> eval(c1)
[1] 2


quote()関数によって表現式オブジェクトを呼び出しに変換できる.引数は,eval()関数を用いたとき,またはオブジェクトをデータとして扱ったときに,はじめて評価される.substitute関数と異なり,quote関数は含まれる引数を評価しないで構文木を返す.

> q1 <- quote(sum(2, 4))
> q1
sum(2, 4)
> is.call(q1)
[1] TRUE
> eval(q1)
[1] 6


as.call()関数はリストを呼び出しに変換する.as.list()関数はその逆.

> f <- sum
> A <- c(1, 3)
> (g <- as.call(list(f, quote(A))))
.Primitive("sum")(A)
> eval(g)
[1] 4
> as.list(g)
[[1]]
function (..., na.rm = FALSE)  .Primitive("sum")

[[2]]
A

リストにcallモードを与えても同じことができる.

> g <- list(f, quote(A))
> eval(g)
[[1]]
function (..., na.rm = FALSE)  .Primitive("sum")

[[2]]
A

> mode(g) <- "call"
> eval(g)
[1] 4


呼び出しはリスト風の構造をしている.

> q2 <- quote(sum(1, 2))
> q2[[1]]
sum
> q2[[2]]
[1] 1
> q2[[3]]
[1] 2
> q2[[1]] <- c
> eval(q2)
[1] 1 2


呼び出しの各成分は名前オブジェクトで置き換えられる.

> q3 <- quote(plot(x, y))
> q3[[1]]
plot
> q3[[1]] <- as.name("+")
> q3
x + y


代入

substitute(expr, env)関数は,環境envに拘束された変数を代入した表現式exprに対する構文木を返す.plot()の軸ラベルなどを生成する際に,評価された実引数の値ではなく,実引数を生成する式が表示されるのは,このはたらきによる.
結果は通常,逆構文解析関数であるdeparse()で文字列に変換され,使われる.

> f1 <- function(x){
+   cat(x, fill = TRUE)
+   cat(deparse(substitute(x)), fill = TRUE)
+ }
> f1(rnorm(2))
1.829917 -0.07375462
rnorm(2)

substitute()関数は,スロットに記録された予約オブジェクトの中身を取り出している (c.f. 遅延評価).


envirに環境やリストを指定すると,その中で変数が検索され,あれば代入される.

> # envirにリストを指定
> (A <- substitute(expression(a + b), list(a = 1, b = 2)))
expression(1 + 2)
> eval(eval(A))
[1] 3
> # シンボルを変更
> (B <- substitute(expression(a + b), list(a = 1, b = quote(c))))
expression(1 + c)
> eval(eval(B))
 以下にエラー 1 + c :  二項演算子の引数が数値ではありません 
> # cは関数c()とされる
> (C <- substitute(expression(a + b), list(a = 1, b = c)))
expression(1 + .Primitive("c"))
> eval(eval(C))
 以下にエラー 1 + .Primitive("c") :  二項演算子の引数が数値ではありません 
> # シンボルに値を代入
> d <- 3
> (D <- substitute(expression(a + b), list(a = d, b = 2)))
expression(3 + 2)
> eval(eval(D))
[1] 5

代入は形式的で,意味を持つかどうかはチェックされない.


遅延評価により,代入式が評価された時点で値が代入されるので,以下のような事態が起こる.

> f1 <- function(y, main = deparse(substitute(y))){
+   y <- log(y)
+   plot(y, main = main) # ここでmainが評価される
+ }
> f2 <- function(y, main = deparse(substitute(y))){
+   main # ここでmainが評価されてしまう
+   y <- log(y)
+   plot(y, main = main)
+ }
> f1(1:3)
> f2(1:3)


呼び出しの即実行

do.call(what, args)関数は,関数名whatと引数リストargsをもとに関数呼び出しを構成し,実行する.lapply()関数との違いは,lapplyはリストの各成分ごとに関数whatを適用して,リストとして結果を返す点.

> (x <- data.frame(A = 1:3, B = -1, C = 0))
  A  B C
1 1 -1 0
2 2 -1 0
3 3 -1 0
> lapply(x, "sum")
$A
[1] 6

$B
[1] -3

$C
[1] 0

> do.call(sum, x)
[1] 3


言語オブジェクトを操作する関数

表現式オブジェクトの作成

expression()関数は表現式オブジェクトを作成する.表現式は,呼び出しオブジェクトを並べたリスト風の構造で,リスト風の操作ができる.

> (ex1 <- expression(sum(1:3), "B", 0))
expression(sum(1:3), "B", 0)
> length(ex1)
[1] 3
> ex1[[2]]
[1] "B"


表現式を指定環境で評価

evalq()関数はeval(quote(expr))と同値で,式中の引数を,現在のスコープではなく,オプションで指定された環境で評価する.

> a <- -1
> b <- 1
> eval(a + b, list(a = 0))
[1] 0
> evalq(a + b, list(a = 0))
[1] 1


テキストから表現式 とその逆

parse()関数は文字列textの入力を構文解析し,表現式のリストにする.eval関数と組み合わせて使うと便利.
deparse()関数は表現式を逆構文解析し,文字列に変換する.fig用のラベル文字列を式から構成するときなどに用いる.

> a <- -1
> b <- 1
> txt1 <- parse(text = "a + b")
> eval(txt1)
[1] 0
> deparse(txt1)
[1] "structure(expression(a + b), srcfile = <environment>, wholeSrcref = structure(c(1L, "
[2] "0L, 2L, 0L, 0L, 0L, 1L, 2L), srcfile = <environment>, class = \"srcref\"))"   


Rの言語オブジェクトについてまとめてみました.
Rプログラミングマニュアル』(間瀬茂, 2007年, 数理工学社) を参考にしました.
言語としてのRの基礎や,豊富な役立つtipsが載っており,とても参考になる本です.


言語オブジェクト

Rの言語オブジェクトには3つの型がある.呼び出し call,表現式 expression,名前 name である.


呼び出しオブジェクトは関数名と引数を与えてつくる.
作成したオブジェクトは,eval()関数で評価実行できる.

> c1 <- call("round", 2.4)
> c1
round(2.4)
> eval(c1)
[1] 2


quote()関数によって表現式オブジェクトを呼び出しに変換できる.引数は,eval()関数を用いたとき,またはオブジェクトをデータとして扱ったときに,はじめて評価される.substitute関数と異なり,quote関数は含まれる引数を評価しないで構文木を返す.

> q1 <- quote(sum(2, 4))
> q1
sum(2, 4)
> is.call(q1)
[1] TRUE
> eval(q1)
[1] 6


as.call()関数はリストを呼び出しに変換する.as.list()関数はその逆.

> f <- sum
> A <- c(1, 3)
> (g <- as.call(list(f, quote(A))))
.Primitive("sum")(A)
> eval(g)
[1] 4
> as.list(g)
[[1]]
function (..., na.rm = FALSE)  .Primitive("sum")

[[2]]
A

リストにcallモードを与えても同じことができる.

> g <- list(f, quote(A))
> eval(g)
[[1]]
function (..., na.rm = FALSE)  .Primitive("sum")

[[2]]
A

> mode(g) <- "call"
> eval(g)
[1] 4


呼び出しはリスト風の構造をしている.

> q2 <- quote(sum(1, 2))
> q2[[1]]
sum
> q2[[2]]
[1] 1
> q2[[3]]
[1] 2
> q2[[1]] <- c
> eval(q2)
[1] 1 2


呼び出しの各成分は名前オブジェクトで置き換えられる.

> q3 <- quote(plot(x, y))
> q3[[1]]
plot
> q3[[1]] <- as.name("+")
> q3
x + y


代入

substitute(expr, env)関数は,環境envに拘束された変数を代入した表現式exprに対する構文木を返す.plot()の軸ラベルなどを生成する際に,評価された実引数の値ではなく,実引数を生成する式が表示されるのは,このはたらきによる.
結果は通常,逆構文解析関数であるdeparse()で文字列に変換され,使われる.

> f1 <- function(x){
+   cat(x, fill = TRUE)
+   cat(deparse(substitute(x)), fill = TRUE)
+ }
> f1(rnorm(2))
1.829917 -0.07375462
rnorm(2)

substitute()関数は,スロットに記録された予約オブジェクトの中身を取り出している (c.f. 遅延評価).


envirに環境やリストを指定すると,その中で変数が検索され,あれば代入される.

> # envirにリストを指定
> (A <- substitute(expression(a + b), list(a = 1, b = 2)))
expression(1 + 2)
> eval(eval(A))
[1] 3
> # シンボルを変更
> (B <- substitute(expression(a + b), list(a = 1, b = quote(c))))
expression(1 + c)
> eval(eval(B))
 以下にエラー 1 + c :  二項演算子の引数が数値ではありません 
> # cは関数c()とされる
> (C <- substitute(expression(a + b), list(a = 1, b = c)))
expression(1 + .Primitive("c"))
> eval(eval(C))
 以下にエラー 1 + .Primitive("c") :  二項演算子の引数が数値ではありません 
> # シンボルに値を代入
> d <- 3
> (D <- substitute(expression(a + b), list(a = d, b = 2)))
expression(3 + 2)
> eval(eval(D))
[1] 5

代入は形式的で,意味を持つかどうかはチェックされない.


遅延評価により,代入式が評価された時点で値が代入されるので,以下のような事態が起こる.

> f1 <- function(y, main = deparse(substitute(y))){
+   y <- log(y)
+   plot(y, main = main) # ここでmainが評価される
+ }
> f2 <- function(y, main = deparse(substitute(y))){
+   main # ここでmainが評価されてしまう
+   y <- log(y)
+   plot(y, main = main)
+ }
> f1(1:3)
> f2(1:3)


呼び出しの即実行

do.call(what, args)関数は,関数名whatと引数リストargsをもとに関数呼び出しを構成し,実行する.lapply()関数との違いは,lapplyはリストの各成分ごとに関数whatを適用して,リストとして結果を返す点.

> (x <- data.frame(A = 1:3, B = -1, C = 0))
  A  B C
1 1 -1 0
2 2 -1 0
3 3 -1 0
> lapply(x, "sum")
$A
[1] 6

$B
[1] -3

$C
[1] 0

> do.call(sum, x)
[1] 3


言語オブジェクトを操作する関数

表現式オブジェクトの作成

expression()関数は表現式オブジェクトを作成する.表現式は,呼び出しオブジェクトを並べたリスト風の構造で,リスト風の操作ができる.

> (ex1 <- expression(sum(1:3), "B", 0))
expression(sum(1:3), "B", 0)
> length(ex1)
[1] 3
> ex1[[2]]
[1] "B"


表現式を指定環境で評価

evalq()関数はeval(quote(expr))と同値で,式中の引数を,現在のスコープではなく,オプションで指定された環境で評価する.

> a <- -1
> b <- 1
> eval(a + b, list(a = 0))
[1] 0
> evalq(a + b, list(a = 0))
[1] 1


テキストから表現式 とその逆

parse()関数は文字列textの入力を構文解析し,表現式のリストにする.eval関数と組み合わせて使うと便利.
deparse()関数は表現式を逆構文解析し,文字列に変換する.fig用のラベル文字列を式から構成するときなどに用いる.

> a <- -1
> b <- 1
> txt1 <- parse(text = "a + b")
> eval(txt1)
[1] 0
> deparse(txt1)
[1] "structure(expression(a + b), srcfile = <environment>, wholeSrcref = structure(c(1L, "
[2] "0L, 2L, 0L, 0L, 0L, 1L, 2L), srcfile = <environment>, class = \"srcref\"))"   

環境 (その2)

Rの環境について,2回に分けてまとめました.
第1回目は,環境全般とスコープ規則についてでした.
第2回目は,環境を操作する関数について.
なお,『Rプログラミングマニュアル』(間瀬茂, 2007年, 数理工学社) を参考にしました.
言語としてのRの基礎や,豊富な役立つtipsが載っており,とても参考になる本です.


環境へのアクセス

sys.関数群は,関数にともなう環境へのアクセス手段を,フレームスタックをたどって提供する.
.GlobalEnvはフレームリストのうちで番号0をもつ.関数が評価されると,フレームスタックが1ずつ増加していく.この評価環境にはsys.関数でアクセスできる.
一方,関数評価の親フレームは,かならずしも現在の評価フレームから1をひいたものではなく,関数が定義された環境でもなく,関数が呼び出された環境になる.

ff <- function(x) gg(x)
gg <- function(y) sys.status()
str(ff(1)) # objectの構造を表示する

> str(ff(1))
List of 3
 $ sys.calls  :Dotted pair list of 4
  ..$ : language str(ff(1))
  ..$ : language ff(1)
  ..$ : language gg(x)
  ..$ : language sys.status()
 $ sys.parents: int [1:4] 0 0 2 3     # gg()の評価環境の親フレーム番号
 $ sys.frames :Dotted pair list of 4
  ..$ :<environment: 0xda4bf0> 
  ..$ :<environment: 0xda4a4c> 
  ..$ :<environment: 0xda49c0> 
  ..$ :<environment: 0xda4988> 


以下のようなfoo1()を実行すると,その評価環境がスタック位置1にできる.
次に,foo1()内部で実行されるfoo2()はスタック位置2にその評価環境をつくる.
最後に,foo2()内部で実行されるfoo3()はスタック位置3にその評価環境をつくる.
foo3()の実行が終了するとその評価環境と内部変数x3は消滅する.
foo2()の実行が終了するとその評価環境と内部変数x2は消滅する.
foo1()の実行が終了するとその評価環境と内部変数x1は消滅し,大局的環境にもどる.大局的環境には変数x0がのこる.

foo1 <- function(){
  x1 <- foo2()
  cat(str(foo1))
  return(x1)
}

foo2 <- function(){
  x2 <- foo3()
  cat(str(foo3))
  return(x2)
}

foo3 <- function(){
  x3 <- pi
  cat(str(x3))
  return(x3)
}


ソースコードの読み込み

関数sys.source()はファイルからソースコードを読み込み,指定された環境に登録する.

foo <- function(x){
  # This is a function.
  return(x)
}
y <- 1:3
save("foo", "y", file = "foo.R", ascii = TRUE)
tempenv <- new.env()
sys.source("foo.R", envir = tempenv, keep.source = TRUE)
ls(env = tempenv)
tempenv$foo
tempenv$y

> ls(env = tempenv)
[1] "foo" "y"  
> tempenv$foo
function(x){
  # This is a function.
  return(x)
}
<environment: 0xd045ac>
> tempenv$y
[1] 1 2 3


永続的付値

永続的付値演算子 <<- は,対象の変数をスコープ規則にしたがって探索し,みつかった時点でそれに付値する.もしなければ,大局環境の.GlobalEnv中に変数を作成し,それに付値する.

> X
 エラー:  オブジェクト 'X' がありません 
> foo1 <- function(x){
+   bar <- function(y){
+     X <<- y
+     return(X)
+   }
+   X <- 0
+   bar(x)
+   cat(X)
+ }
> foo1(pi)
3.141593
> X
 エラー:  オブジェクト 'X' がありません 
# foo1()の評価環境中のXに値が付値されている.

> bar <- function(y){
+   X <<- y
+   return(X)
+ }
> 
> foo2 <- function(x){
+   X <- 0
+   bar(x)
+   cat(X)
+ }
> 
> foo2(1:3)
0
> X
[1] 1 2 3
# こちらは大局環境にXがつくられ,それに付値されている.

特定環境への値の付値

assign()関数はある環境中の与えられた名前のオブジェクトに値を付値する.
envir引数が与えられなかった場合,付値は現在アクティブな環境で行われる (よって,関数中でそのままassignしても,関数評価の終了とともにオブジェクトは消滅する).
もしinherit = TRUEが引数で与えられると,対象の変数がみつかるまで囲み環境が探索され,みつかった環境中で付値される.

> for(i in 1:6){
+ 	nam <- paste("n", i, sep = "")
+ 	assign(nam, 1:i)
+ }
> ls(pat = "n.$")
[1] "n1" "n2" "n3" "n4" "n5" "n6"
> myf <- function(x){
+   inner.f <- function(x){
+     assign("x.global", x^2, env = .GlobalEnv)
+   }
+   inner.f(x + 1)
+ }
> myf(3)
> x.global
[1] 16


特定環境からの値の取得

get()関数は与えられた名前のオブジェクトを指定された環境中で探し,値を返す.
pos引数には探索対象の環境を指定する.
mode引数には探したい型のオブジェクトが指定できる.

> get("get")
function (x, pos = -1, envir = as.environment(pos), mode = "any", 
    inherits = TRUE) 
.Internal(get(x, envir, mode, inherits))
<bytecode: 0xc784fc>
<environment: namespace:base>


特定環境へのアクセス

確認

exists()関数は,与えられた名前の変数が指定した環境に存在するかどうかを確認する.
where引数には探索対象の環境を指定する.
mode引数には探したい型のオブジェクトが指定できる.

> search()
 [1] ".GlobalEnv"        "tools:RGUI"        "package:stats"    
 [4] "package:graphics"  "package:grDevices" "package:utils"    
 [7] "package:datasets"  "MacJapanEnv"       "package:methods"  
[10] "Autoloads"         "package:base"     
> aa <- 1
> exists("aa", 1) # aaは位置1 (.GlobalEnv) にある
[1] TRUE
> exists("aa", 2) # aaは位置2 (tools:RGUI) にはない
[1] FALSE
> exists("lm", 3) # lmは位置3 (package:stats) にある
[1] TRUE


アクセス

特定環境中の変数には,2重コロン演算子でアクセスできる.

> pi <- 1:3
> pi
[1] 1 2 3
> base::pi # 本来の値が返される
[1] 3.141593

パッケージの場合は,読み込んだうえで値を返してくれる.

> head(SpaceShuttle)
 以下にエラー head(SpaceShuttle) :  オブジェクト 'SpaceShuttle' がありません 
> head(vcd::SpaceShuttle)
  FlightNumber Temperature Pressure Fail nFailures Damage
1            1          66       50   no         0      0
2            2          70       50  yes         1      4
3            3          69       50   no         0      0
4            4          80       50 <NA>        NA     NA
5            5          68       50   no         0      0
6            6          67       50   no         0      0


名前空間

attachNamespaces(),loadNamespaces(),loadedNamespaces()などの関数で名前空間を操作できる.これらはlibrary()関数の内部で使われている.

> loadedNamespaces()
[1] "base"      "graphics"  "grDevices" "MASS"      "methods"  
[6] "stats"     "utils"     "vcd"  


パッケージ中の変数にアクセス

pkg::name は,パッケージpkgのエクスポートされた変数nameの値を返す.
pkg:::name は,パッケージpkgの内部変数nameの値を返す.

> base::"+"
function (e1, e2)  .Primitive("+")


そのほか

topenv()関数でトップレベルの環境をみつける.


同じ環境を複数登録できる.しかし実行速度が遅くなるので避けるべき.

> attach(iris); attach(iris); attach(iris)
> search()
 [1] ".GlobalEnv"        "iris"              "iris"             
 [4] "iris"              "tools:RGUI"        "package:stats"    
 [7] "package:graphics"  "package:grDevices" "package:utils"    
[10] "package:datasets"  "MacJapanEnv"       "package:methods"  
[13] "Autoloads"         "package:base"     


多数の関数が共有するオブジェクトは,一次作業環境に置いておくと便利.
コピーをつくらずに直接参照でき,メモリの節約にもなる.

> tempenv <- new.env()
> assign("aa", 1, env = tempenv)
> foo <- function(x){x^2}
> assign("foo", foo, env = tempenv)
> rm(foo)
> aa
 エラー:  オブジェクト 'aa' がありません 
> get("aa", env = tempenv)
[1] 1
> get("foo", env = tempenv)
function(x){x^2}
> tempenv$foo(4) # リスト風にアクセスできる
[1] 16

集団遺伝学 (遺伝的浮動2)

集団遺伝学の基礎について,授業の復習をメインにまとめてみます.
「適応進化遺伝学」テキストを参考にしています.
自分の学習用ですので,間違いなどあるかもしれません.

遺伝的浮動1からのつづき

有効集団サイズ

遺伝的浮動に関係するのは,繁殖まできちんとたどりつけた個体の数で,たいていこれは見かけ上の集団個体数より小さい.この問題を一般的に扱う概念が有効集団サイズ (effective population number) であり,記号N_eで表記される.

たとえば,
式8より,(1 - \frac{1}{2N_e})^t = (1 - \frac{1}{2N_{t-1}})(1 - \frac{1}{2N_{t-2}}) \cdots (1 - \frac{1}{2N_0})
であり,右辺は(1 - \frac{t}{2N_e}),左辺は1 - \sum^{t-1}_{i=0}\frac{1}{2N_i}近似できる.
よって
\frac{1}{N_e} = \frac{1}{t}(\frac{1}{N_1} + \frac{1}{N_1} + \cdots + \frac{1}{N_t})
となり,どんなに普段Nが大きくても,1回でも急激な集団縮小がおこれば,有効集団サイズは急に小さくなる (それゆえ変異性が大きく影響される) ことになる.これは,ボトルネック効果を表している.


また,雄と雌で繁殖個体数に差がある場合,雄と雌の個体数をそれぞれN_mN_fとすると,
N_e = \frac{4N_m N_f}{N_m + N_f}
が得られる.
つまり,雄が無限にいても,繁殖に関わる雌が1個体しかいなければ,N_e = 4にしかならない.


coalescence time

集団中の遺伝子がcoalescenceするまでにさかのぼらなければならない時間の期待値.逆に考えれば,集団に1個生じた遺伝子が,遺伝的浮動によって固定するまでにかかる時間の期待値も求められる.

集団中のn個の遺伝子のどれかに,t世代前にはじめてcoalescenceがおこる待ち時間をtとし,その確率分布をP(t)とする.t世代前までのどこかでcoalescenceする確率であるC_tF_t (突然変異がおこらないと仮定した場合) とは異なるので注意.


最後の1世代以外はずっとcoalescenceしない確率なので,1世代のあいだにcoalescenceがおこる確率をC_1とおくと,
coalescence確率は P(t) = C_1 (1 - C_1)^{t-1}
待ち時間の期待値は \overline{t} = \sum^{\infty}_{t=1} {}_tC_1(1 - C_1)^{t-1}
で与えられる.
1世代前でcoalescenceしない確率1 - C_1Dとおくと,
\overline{t} = C_1 (1 + 2D + 3D^2 + \cdots)
ここで S_n = 1 + 2D + 3D^2 + \cdots+ nD^{n-1} とおいて,S_n - DS_n を考えると,
 (1 - D) S_n = 1 + D + D^2 + \cdots + D^{n-1} - nD^n
さらに X = 1 + D + D^2 + \cdots+ D^{n-1} とおいて,X - DX を考えると,
X = \frac{1 - D^n}{1 - D}
したがって,
 (1 - D)S_n = \frac{1 - D^n}{1 - D} - nD^n すなわち
 S_n = \frac{1 - D^n}{C_1^2} - \frac{nD^n}{C_1}
ここで,\lim_{n\rightarrow\infty}nDn = 0\lim_{n\rightarrow\infty}D^n = 0なので,
 \lim_{n\rightarrow\infty}S_n = \frac{1}{C_1^2}
よって,
\overline{t} = \frac{1}{C_1} [式10]
つまり,coalescenceがおこるまでに待ち時間の期待値は,1世代でcoalescenceする確率の逆数となる.


2点のサンプルが1世代前にcoalescenceする確率はC_1[n = 2] = \frac{1}{2N_e}なので,集団中から無作為に選んだ2つの遺伝子がcoalescenceする待ち時間の期待値は,
\overline{t_2} = 2N_e 世代となる.


次に,n点のサンプルが1個にcoalescenceする待ち時間を考える.1世代前にn個のサンプルのどれかにcoalescenceがおこる確率C_1(n)は,どれもcoalescenceしない場合の補集合であるから,
C_1(n) = 1 - \frac{2N_e - 1}{2N_e} \frac{2N_e - 2}{2N_e} \frac{2N_e - 3}{2N_e} \cdots \frac{2N_e - (n - 1)}{2N_e}
N_eが十分大きく,nN_eより十分小さければ,近似を用いて
C_1(n) \approx 1 - (\frac{1}{2N_e}(1 + 2 + 3 + \cdots + n-1)) = \frac{1}{2N_e} \frac{n (n - 1}{2} = \frac{n (n - 1)}{4N_e}
となる.3つ以上が同時にcoalescenceする確率は無視できるほど小さいと考えている.


よって式10より,n個のサンプルがn-1個にcoalescenceする待ち時間の期待値は
\overline{t_n} = \frac{4N_e}{n(n-1)}となる.
また,n個のサンプルがすべて1個にcoalescenceする待ち時間 (Most Recent Common Ancesterまでの時間) の期待値は,
T_{MRCA} = \sum_{i=2}^n \overline{t_1} = \sum_{i=2}^n \frac{4N_e}{i(i-1)} = 4N_e \sum_{i=2}^n (\frac{1}{i - 1} - \frac{1}{i}) = 4N_e(1 - \frac{1}{n}) 世代となる. [式11]
サンプル数が十分大きければ,集団の全遺伝子がcoalescenceする待ち時間の期待値は4_N_e世代となる.


突然変異と遺伝的浮動

突然変異によって新しい変異が供給されるときに,遺伝的浮動がどうなるか考える.

固定指数

世代あたりの突然変異率をu,有効集団サイズをN_eとする.突然変異がおこるという仮定のもとでも,集団から無作為抽出した2つの遺伝子がt世代前までにcoalescenceする確率は,式8と同様に
C_t = 1 - (1 - \frac{1}{2N_e})^t [式12]
となる.
その一方,固定指数F_tはこの式では表せない (対象とする配列が同一であることを前提とするため.突然変異があれば,共通祖先の複製に由来するものでも配列が異なり得る).この場合F_tは,集団から無作為抽出した2つの遺伝子がt世代前までのどこかの複製に由来し,かつ突然変異を免れ続けている確率になるので,
F_t = (\frac{1}{2N_e} + (1 - \frac{1}{2N_e})F_{t-1})(1 - u)^2 [式13]
となる.これは漸化式なので,t-1世代までのすべてで突然変異を免れる確率が前提に含まれている.


突然変異によって,同一コピー数はどこかで平衡状態になる.平衡点では,複製コピーが増えた分だけ突然変異が生じるようなイメージ.
平衡状態ではF_t = F_{t-1}なので,これを\hat{F}と書いて式13を変形すると,
\hat{F} = (\frac{1}{2N_e} + (1 - \frac{1}{2N_e})\hat{F})(1 - u)^2 \\ = \frac{(1 - u)^2}{2N_e - (2N_e - 1)(1 - u)^2} \\ \approx \frac{1 -2u}{4N_e u + 1 - \\ \approx \frac{1}{4N_e u + 1} \\ = \frac{1}{\Theta + 1} [式14]
となる (u << 1のため2u = 0u^2 = 0と近似).
4N_em = \Thetaは集団の多様性を決定する要となる重要なパラメータで,集団の1世代あたりにあらわれる変異の数なので,集団突然変異率 (population mutation rate) とも呼ばれる.


平衡状態でのヘテロ接合度

有限時間内で,突然変異率が十分小さければ,まったく同一の変異はおこらないと考えて良く,常に新規のalleleがつくりだされる (無限アレルモデル infinite allele model).突然変異と遺伝的浮動のもとでの集団内多様性は,1-C_tの割合のもともと異なる遺伝子と,複製に由来するが突然変異を受けたC_t - F_tの割合の遺伝子から生じる.平衡状態での集団のヘテロ接合度を,式9のように平衡状態に達した多くの分集団のヘテロ接合度の平均値\overline{H_S}と考えると,
\hat{H} = \overline{H_S} = H_T(1-C_t) + C_t - F_t [式15]
が与えられる.突然変異がなければC_t = F_tで,式9と式15は同じになる.
平衡状態であればC_t = 1なので,式15より,
\hat{H} = 1 - F_t = \frac{4N_em}{4N_em + 1} = \frac{\Theta}{\Theta + 1} [式16]
となる.


N_em << 1つまりu << \frac{1}{N_e}の場合,突然変異の効果が遺伝的浮動の効果よりずっと小さい.このとき\hat{H} \leftarrow 0となり,あるalleleへの固定がおこる.
N_em >> 1つまりu >> \frac{1}{N_e}の場合,突然変異の効果が遺伝的浮動の効果よりずっと大きい.このとき\hat{H} \leftarrow 1となり,集団には常に多様性があることになる.
4N_em >> 1の状態なら自,然選択のはたらきなしでも,集団内に非常にたくさんの変異が保持されることになる.ゲノムレベルで見るとヘテロ接合度は実質的に1 (非常に大きい値) となり,中立論の論拠のひとつを与えている.


平衡状態でのallele種類数

突然変異と遺伝的浮動の平衡状態にある集団では,無作為抽出したn個の遺伝子に何種類のalleleがあるかを予測できる.alleleの数をkとして,その期待値は,
E(k) = 1 + \frac{\Theta}{\Theta + 1} + \frac{\Theta}{\Theta + 2} + \cdots + \frac{\Theta}{\Theta + n - 1} [式17]
となる (Ewensの式).


個体群動態の検証

遺伝的浮動と突然変異の平衡状態では,allele種類数の頻度分布が予測できる.この分布は,alleleがk種類あるときn個のサンプル中にそれぞれが何個ずつあるかの確率
Pr\{n1, n2, \cdots ,n_k\} = \frac{n! \Theta^k}{k! n_1 n_2 \cdots n_k S_n(\Theta)} [式18]
ただし S_n(\Theta) = \Theta (\Theta + 1)(\Theta + 2) \cdots (\Theta + n -1)
から求められる (Ewens-Wattersonの式).
このallele頻度スペクトルは,集団サイズが世代間で一定かつ,中立選択を前提としている.したがって,観察データがこの式の予測とずれた場合には,これらの前提が成り立っていないことが示唆される.
すなわち,
低頻度alleleの種類数が多く,高頻度alleleの種類数が少ない場合:
  • ボトルネック創始者効果がおこって変異性が低下してから十分に時間が経っていない
  • 集団サイズが拡大中で平衡に達していない
  • 遺伝子頻度構成の大きく異なる集団が少数だけ流入して間もない
  • selettive sweepや弱有害による浄化淘汰がはたらいている
低頻度alleleの種類数が少なく,高頻度alleleの種類数が多い場合:
  • 集団サイズが縮小しており,低頻度alleleの消失が優先的におきている
  • 遺伝子頻度構成の大きく異なる集団が多数混合して間もない
  • 平衡選択がおきている
となる.


移住と遺伝的浮動

変異性は移住によっても供給され,一般に移住率は突然変異率よりずっと大きい.また,移住は突然変異と違って,集団分化を妨げる方向にはたらく.


固定指数と移住

単純なisland model (分岐した多くの分集団があり,各集団から割合mの遺伝子を無作為抽出→混ぜあわせ→再分配) のもとでの移住パターンを考える.
mは世代あたりの移住率,N_eは有効集団サイズとする.C_t (集団から無作為抽出した2つの遺伝子がt世代前までにcoalescenceする確率) は,式13のときと同様に,
C_t = (\frac{1}{2N_e} + (1 - \frac{1}{2N_e})C_{t-1})(1 - m)^2 [式19]
で与えられる.突然変異を無視すれば,固定指数F_tは[texC_t]と同じになり,
F_t = (\frac{1}{2N_e} + (1 - \frac{1}{2N_e})F_{t-1})(1 - m)^2 [式20]
これは漸化式なので,t-1世代までのすべてで移住を免れる確率が前提に含まれている.


移住と遺伝的浮動が平衡状態になると,遺伝的浮動によって複製コピーが増えた分だけ,それらが移出しかわりに移入がおこる,というような状態になる.
平衡状態ではC_t = C_{t-1}なので,これを\hat{C}として式13のときのように変形すると,
\hat{C} = \hat{F} \approx \frac{1}{4N_em + 1} [式21]
が得られる.

平衡状態でのヘテロ接合度

island modelのもとでは,移入によって供給される遺伝子は総合集団のallele頻度と同一になる.そのため,平衡状態でのヘテロ接合度\hat{H}は式15・16と同様に考えて,
\hat{H} \approx \frac{4N_em}{4N_em + 1} H_T [式22]
となる.


N_em << 1つまりu << \frac{1}{N_e}の場合,移住の効果が遺伝的浮動の効果よりずっと小さい.このとき\hat{H} \leftarrow 0となり,移住を考慮しない場合と同じになる (あるalleleへの固定がおこる).
N_em >> 1つまりu >> \frac{1}{N_e}の場合,突然変異の効果が遺伝的浮動の効果よりずっと大きい.このとき\hat{H} \leftarrow 1となり,集団分化していないのと同じになる (集団には常に多様性があることになる).
4N_em >> 1の状態なら自,然選択のはたらきなしでも,集団内に非常にたくさんの変異が保持されることになる.ゲノムレベルで見るとヘテロ接合度は実質的に1 (非常に大きい値) となり,中立論の論拠のひとつを与えている.


移住だけでなく突然変異を考慮しても,C_tの漸化式は式19と,平衡式は式21と変わらない.
一方,固定指数は変化する.ある遺伝子が移入でも突然変異でもない確率は
(1 - m)(1 - u) \approx 1 - (m + u)
となり (mu << m + u << 1という近似を利用),
漸化式は
F_t = (\frac{1}{2N_e} + (1 - \frac{1}{2N_e})F_{t-1})(1 - (m + u))^2 [式23]
平衡式は
\hat{F} \approx \frac{1}{4N_e (m + u) + 1} [式24]
となる.
通常はm >> uであるため,
\hat{F} \approx \hat{F} \approx \frac{1}{4N_em + 1}
\hat{H} = \frac{4N_em}{4N_em + 1} H_T
と考えて差し支えない.


移住率

移住と遺伝的浮動が平衡に達しているとき,F_{ST} = \hat{F}とみなすことができる.式21より,
N_em = \frac{1}{4}(\frac{1}{F_{ST}} - 1) [式25]
と推定される.
複数の集団の遺伝子頻度を測定し,ヘテロ接合度をHWEから推定することで,island modelのもとでの世代あたり移入個体数を推定できる.多くの仮定が含まれるが,移住率の有効な目安になる.

環境 (その1)

Rの環境について,2回に分けてまとめます.


第1回目は,環境全般とスコープ規則について.
第2回目は,環境を操作する関数についての予定.


なお,『Rプログラミングマニュアル』(間瀬茂, 2007年, 数理工学社) を参考にしました.
言語としてのRの基礎や,豊富な役立つtipsが載っており,とても参考になる本です.


環境

環境は,frame (オブジェクトのシンボル名と対になる値の集合) と,その囲み環境 (enclosing environment) へのポインタであるenclosureからなる.
囲み環境には,関数 parent.env() でアクセスできる.
環境はツリー構造をしている.ルートノードは emptyenv() という空の環境で,R base package の囲み環境になっている.

> globalenv()
<environment: R_GlobalEnv>
> parent.env(.GlobalEnv)
<environment: 0x1f88678>
attr(,"name")
[1] "tools:RGUI"
> baseenv()
<environment: base>
> parent.env(baseenv())
<environment: R_EmptyEnv>


関数 new.env() は,指定した環境を囲み環境とする新たな空の環境をつくる.

> (env <- new.env(parent = .GlobalEnv))
<environment: 0x1601c9f8>
> parent.env(env)
<environment: R_GlobalEnv>


環境中のオブジェクトには,関数 ls(), get(), assign(), eval() などでアクセスできる.

> ls(.GlobalEnv)
[1] "env"    "R_LIBS"
> head(ls(baseenv()))
[1] "^"        "-"        "-.Date"   "-.POSIXt" ":"        "::"      


環境は他のオブジェクトと異なり,参照渡しで,コピーされない.よって,複数の変数が同一環境を指していた場合,一方を変更すると他方も変更される.

> # Environment.
> env1 <- new.env()
> env2 <- env1
> assign("x", rnorm(3), envir = env1)
> ls(envir = env2)
[1] "x"
> get("x", envir = env1)
[1] -0.29572096 -0.06950039  0.80579259
> get("x", envir = env2)
[1] -0.29572096 -0.06950039  0.80579259
> # Object.
> obj1 <- numeric(3)
> obj2 <- obj1
> assign("obj1", rnorm(3))
> get("obj1", envir = .GlobalEnv)
[1] -0.8161385 -0.5597408 -0.6494503
> get("obj2", envir = .GlobalEnv)
[1] 0 0 0


サーチパス

現在のパスは関数 search() で得られる.

> search()
 [1] ".GlobalEnv"        "tools:RGUI"        "package:stats"    
 [4] "package:graphics"  "package:grDevices" "package:utils"    
 [7] "package:datasets"  "MacJapanEnv"       "package:methods"  
[10] "Autoloads"         "package:base"     


ある環境中のオブジェクト一覧は関数 ls() で見られる.

> head(ls(pos = "package:stats"))
[1] "acf"        "acf2AR"     "add.scope"  "add1"       "addmargins"
[6] "aggregate" 
> head(ls(pos = 3))
[1] "acf"        "acf2AR"     "add.scope"  "add1"       "addmargins"
[6] "aggregate" 


attach() で登録されたオブジェクトや読み込まれたパッケージは,.GlobalEnv以下に,読み込まれた順に付加されていく.

> library(glmmML)
> x <- data.frame(x = rnorm(3))
> attach(x)
The following object(s) are masked _by_ '.GlobalEnv':

    x
> search()
 [1] ".GlobalEnv"        "x"                 "package:glmmML"   
 [4] "tools:RGUI"        "package:stats"     "package:graphics" 
 [7] "package:grDevices" "package:utils"     "package:datasets" 
[10] "MacJapanEnv"       "package:methods"   "Autoloads"        
[13] "package:base"     
> parent.env(as.environment("package:glmmML"))
<environment: 0x1972078>
attr(,"name")
[1] "tools:RGUI"
> parent.env(as.environment(2))
<environment: package:glmmML>
attr(,"name")
[1] "package:glmmML"
attr(,"path")
[1] "/Library/Frameworks/R.framework/Versions/2.14/Resources/library/glmmML"


名前空間

他の環境からのマスクに影響されることなく,パッケージ内でオブジェクトの参照を解決する仕組み.

> IQR # package:stats の中にあり,公開されている.
function (x, na.rm = FALSE, type = 7) 
diff(quantile(as.numeric(x), c(0.25, 0.75), na.rm = na.rm, names = FALSE, 
    type = type))
<bytecode: 0x1e7aa78>
<environment: namespace:stats>
> predict.nls # package:stats の中にあるが,公開されていない.
 エラー:  オブジェクト 'predict.nls' がありません 
> getAnywhere(predict.nls)
A single object matching ‘predict.nls’ was found
It was found in the following places
  registered S3 method for predict from namespace stats
  namespace:stats
with value

function (object, newdata, se.fit = FALSE, scale = NULL, df = Inf, 
    interval = c("none", "confidence", "prediction"), level = 0.95, 
    ...) 
{
    if (missing(newdata)) 
        return(as.vector(fitted(object)))
    if (!is.null(cl <- object$dataClasses)) 
        .checkMFClasses(cl, newdata)
    object$m$predict(newdata)
}
<bytecode: 0xdcc014>
<environment: namespace:stats>


closure環境,評価環境

関数が定義されると,それが定義された環境を含むclosure環境がつくられる.
関数が呼び出されると,その評価 (evaluation) 環境が,closure環境を囲み環境としてつくられる (必ずしもその環境を呼び出した環境ではない).

> a <- rnorm(3)
> foo <- function(x){round(x)}
> environment(foo)
<environment: R_GlobalEnv>
> environment(foo) <- baseenv() # 環境を変更
> a
[1]  0.7547118 -0.6926263 -0.2032858
> foo(a)
[1]  1 -1  0


関数内部である変数が必要になると,評価環境,そのclosure環境,さらにその囲み環境…と順に検索され,大局的環境またはその関数を含むパッケージ環境に至ると,あとはサーチパスに従って規定環境まで検索され,結局みつからなければ最後は空環境にもどり,エラーがでる.

> foo <- function(x){
+   y <- 10
+   x + y
+ }
> foo(1:10)
 [1] 11 12 13 14 15 16 17 18 19 20
> y <- -1
> foo(1:10) # foo の環境内に定義された y が使われる.
 [1] 11 12 13 14 15 16 17 18 19 20
> bar <- function(x){x + y}
> y <- 10
> bar(1:10)
 [1] 11 12 13 14 15 16 17 18 19 20
> y <- -1
> bar(1:10) # グローバル環境に定義された y が使われる.
 [1] 0 1 2 3 4 5 6 7 8 9


スコープ

アクティブな環境へは,呼び出しスタックを通じてアクセスできる.関数が呼び出されるたびにcontextという内部構造がつくられ,スタックされる.関数評価が終わると,contextは呼び出しスタックから取り除かれる.利用できる変数を呼び出しスタックの上部に,時間軸に沿って定義することをdynamic scopeと呼び,関数が定義された環境中の値が使われるのはstatic scope (lexical scope) である.呼び出しスタックへのアクセスには,一連のsys. 関数が使われる.
スコープ規則とは,未拘束変数の値をみつける際に使われる規則群のこと.関数の仮引数は拘束変数,関数内で定義された変数は局所変数,その他は未拘束変数である.未拘束変数に対しては,関数の付属環境 → その囲み環境 → さらにその囲み環境 と,最終的に大局的環境に至るまで検索がなされる.

> foo <- function(x){
+   y <- 10
+   bar <- function(x){x + y + z}
+   return(bar)
+ }
> (foo2 <- foo())
function(x){x + y + z}
<environment: 0xd43440>
> get("bar", envir = environment(foo2))
function(x){x + y + z}
<environment: 0xd43440>
> foo2(1)
 以下にエラー foo2(1) :  オブジェクト 'z' がありません 
> foo2(1, 2)
 以下にエラー foo2(1, 2) :  使われていない引数 (2) 
> z <- 2
> foo(1)
function(x){x + y + z}
<environment: 0xe66f64>
> foo2(1)
[1] 13


関数引数

関数評価の最初のステップで,仮引数と実引数のマッチングが3段階で行われる.
1. 名前タグをもつ仮引数と実引数の完全マッチング (複数該当した場合はエラー).
2. 部分マッチング (複数該当した場合はエラー).
3. 引数位置によるマッチング.
関数引数は値渡し (call-by-value) が原則で,引数は与えられた値で初期化され,仮引数の名前をもつ局所変数であるかのようにふるまう.
関数引数は遅延評価 (lazy evaluation) であり,必要になるまで評価されない.

> foo <- function(x){
+   y <- 1
+   print(y)
+   x * y
+ }
> foo(y <- 0) # 関数内のyは変更されないが,
[1] 1
[1] 0
> y # 大局環境のyは変更されている.
[1] 0


予約オブジェクト

遅延評価のコアとなる特殊なオブジェクト.ある関数が呼び出されたときに.引数が照合され,それぞれの仮引数が予約 (promise) オブジェクトに結びつけられ,仮引数に与えられた表現式と,その関数が呼び出された環境へのポインタが保存される.引数が実際にアクセスされると,保存された表現式が,ポインタの示す環境中で評価され,その値が予約に保存される.以降のアクセスでは,その保存された値が使われる.
関数substitute() は予約から中身の表現式を取り出す.関数delayedAssign() は表現式から予約オブジェクトをつくりだす.

> delayedAssign("x", {y <- 2; cat("These are", y, "numbers.\n"); runif(y)})
> x # 1回目は表現式全体を評価.
These are 2 numbers.
[1] 0.5436929 0.7225412
> x # 2回目は代入された値のみ評価.
[1] 0.5436929 0.7225412

予約は,その値が必要になった時点で強制評価 (forcing) される.

> foo <- function(x, y = deparse(x)){
+   x <- 0
+   print(y) # ここで初めてyが評価される.
+ }
> foo(1)
[1] "0"
> 
> bar <- function(x, y = deparse(x)){
+   y # ここですでにyが評価される.
+   x <- 0
+   print(y)
+ }
> bar(1)
[1] "1"


第2回へつづく.

集団遺伝学 (遺伝的浮動1)

集団遺伝学の基礎について,授業の復習をメインにまとめてみます.
「適応進化遺伝学」テキストを参考にしています.
自分の学習用ですので,間違いなどあるかもしれません.


遺伝的浮動

集団の大きさが有限であれば,偶然の効果によるゆらぎの影響力が無視できなくなってくる.
ここでは,集団サイズを有限としたときに,「偶然」だけでどこまでのことがおこり得るかを,確率論的なモデルで取り扱う.

概念

遺伝的浮動 (random genetic drift) とは,遺伝子頻度の偶然による変動のこと.自然集団で,実際に次代をつくるのに寄与する配偶子が少数に限られることから生じる.抽出効果 (sampling effect) ともいう.
モデルとしては,配偶子プールは大きさ無限大だが,次代の有限なN個体に寄与する配偶子として抽出されるのは,雄性・雌性配偶子あわせて2N個のみ (2倍体雌雄生物の場合),というようなもの.ポイントは,抽出された配偶子中の遺伝子頻度と配偶子プールの遺伝子頻度が必ずしも一致しないこと.この過程が繰り返せば,偶然だけでallele頻度が変化していく.

遺伝子頻度の変化

ひとつの集団中において,あるallele Aの初期頻度をp_0とする.抽出した2N個の配偶子中にi個のAが含まれる確率は二項分布にしたがうので,
P(\frac{i}{2N}) = {}_{2N}C_i p_0^i (1 - p_0)^{2N-i}
となる.


Nが十分大きくp_0が十分小さい場合,確率はポアソン分布で近似でき,
 P(\frac{i}{2N}) = \frac{m^i}{i!} e^{-m}
と書ける.mは期待値で,2Np_0 (下記参照).


1世代後のAの個数の期待値は2Np_0なので,遺伝子頻度の期待値は
E(p_1) = p_0
つまり,初期頻度から変化しない.


一方,1世代後のAの個数の分散は2Np_0 (1 - p_0)なので※1,遺伝子頻度の分散は
 V(p_1) = \frac{p_0 (1 - p_0)}{2N}
つまり,集団が小さいほど大きい.

このことから,世代を経ても遺伝子頻度の期待値は初期頻度と変わらないが,分散は世代時間に対して直線的に大きくなっていく.ということがわかる.
ちなみに,\frac{t}{2N} << 1のとき,遺伝子頻度の分散の期待値は,
V(p_t) = p_0 (1 - p_0) (1 - (1 - \frac{1}{2N})^t) [式1]
であり,これは
V(p_t) = p_0 (1 - p_0) \frac{t}{2N} [式2]
と近似できるそう※2


拡散モデル

初期頻度から出発して,遺伝的浮動を経て固定に至るまでの過程は,拡散モデルで近似されている (近似しないと,膨大な計算量が必要になる).拡散モデルを理解するには高度な数学の知識が必要になるそうなのでここでは取り上げないが,結果に関しては後に触れる.


集団分化

遺伝子頻度の分散

遺伝的浮動は,集団の遺伝子頻度構成の分化をもたらす.そして,集団分化の程度は,集団間での遺伝子頻度の分散として表現できる (分散が大きいほど分化してると考えることができる).式1より,遺伝子頻度の分散は0 (最小: どの集団も同じ頻度構成) から p_0 (1 - p_0) (最大: 全集団のうちp_0ではAが固定,1-p_0ではAが消失) までの値をとり得る.そこで,最大値に対する分散の値が以下で定義され,
 \frac{V(p\)}{p_0 (1 - p_0)} = F_{ST} = 1 - (1 - \frac{1}{2N})^t [式3]
これが分集団化の指標 (遺伝的浮動の指標でもある) となる.


ヘテロ接合度

新たな変異が供給されなければ,遺伝的浮動の過程は,集団の多様性を減少させていく.これはヘテロ接合度が減少していくことでもあるので,式3はヘテロ接合度を使って書き直せる.


まず,分散の公式は以下のとおり.
V(p\) = E(p^2) - E(p\)^2


ここで,E(p^2)はHWEを仮定したallele Aの,分集団でのホモ接合度の期待値とみなすことができる (!).そして E(p\)は,総合集団 (すべての分集団の混合集団) の遺伝子頻度であり, E(p\)^2は総合集団のホモ接合度の期待値とみなすことができる.(ここ重要)


分散の公式から導かれるのは,遺伝的浮動が進むと,分集団のホモ接合度は総合集団のホモ接合度より大きくなっていく (固定や消失がおこっているからですね),ということ.
すべてのalleleについて上記の式を適用し,足し合わせると,以下のようになる.
\Sigma V(p\) = \Sigma E(p^2) - \Sigma E(p\)^2 [式4]


式4の左辺は,式3を代入して展開すると,
 \Sigma V(p\) = F_{ST} \Sigma p_0 (1 - p_0)\\ = F_{ST}(\Sigma p_0 -  \Sigma p_0^2)\\ = F_{ST}(1 -  \Sigma p_0^2)
となる.
\Sigma p_0^2は,HWE下の総合集団での時間0におけるホモ接合度を意味するので,1 - \Sigma p_0^2は同条件下での総合集団のヘテロ接合度  H_T (total population のT) である.
よって,式4の左辺は
\Sigma V(p\) = H_T F_{ST}
と書ける.


一方,式4の右辺の\Sigma E(p^2)は,HWEを仮定した分集団のホモ接合度の期待値なので,分集団のヘテロ接合度 \overline{H_S} (subpopulationのS)を用いて1 - \overline{H_S}と書ける.右辺の \Sigma E(p\)^2は同じく総合集団のホモ接合度の期待値なので,1 - H_Tと書ける.


これより,式4は
H_T F_{ST} = \overline{H_S}
と書き表すことができ,この式より
F_{ST} = \frac{H_T - \overline{H_S}}{H_T} [式5]
が得られる.F_{ST}は,遺伝的浮動の過程 (分集団化の進行) で,分集団のヘテロ接合度が総合集団のヘテロ接合度より小さくなっていく程度を表す指標である.


さらに,式5より
 \overline{H_S} = H_T (1 - F_{ST})
であり,式1を利用すると
 \overline{H_S} = H_T (1 - \frac{1}{2N})^t
が得られる.分集団のヘテロ接合度は,毎代 \frac{1}{2N}の割合で低下していくということであり,集団サイズが小さいほど速く0に近づくことを示している.


集団の階層構造

分集団内の実際のヘテロ接合度H_I (individualsのI) とHWEから推定したヘテロ接合度H_Sの違いの大きさを
F_{IS} = \frac{H_S - H_I}{H_S}
という指標で評価できる.
生殖隔離などがあり,集団のなかにさらに分集団がある (集団に階層構造がある) とき,各分集団ではalleleが固定・消失して多様性が減少していくのに対し,それぞれの分集団は異なるalleleを保持しているので,集団全体で見るとalleleは多様であり得る.このとき,HWEから計算される集団全体のヘテロ接合度は,実際の値より高くなる.
F_{IS}は集団の階層構造という観点より,HWEからの逸脱を評価する指標である.


分集団間での平均をとり,
F_{IS} = \frac{\overline{H_S} - \overline{H_I}}{\overline{H_S}}
と定義すると,
1 - F_{IS} = \frac{\overline{H_I}}{\overline{H_S}} であり,
式5より
1 - F_{ST} = \frac{\overline{H_S}}{H_T} でもあるから,
(1 - F_{IS})(1 - F_{ST}) = \frac{\overline{H_I}}{H_T} となる.
(1 - F_{IS})(1 - F_{ST}) = 1 - F_{IT} [式6]
と定義すると,
F_{IT} = \frac{H_T - \overline{H_I}}{H_T} [式7]
となる.
式7は,HWEから推定した総合集団のヘテロ接合度H_Tと,実際のヘテロ接合度の分集団平均\overline{H_I}の違いの大きさを評価する指標である.
分集団化は,交配範囲が限定されているという意味で,近親婚構造とみなすこともできる.


固定指数

用語

・集団から抽出した2つの遺伝子が,共通の祖先遺伝子に由来することを,coalescenceするという.
・ある配列セットに対し,t世代前までのあいだのどこかでcolescenceがおこる確率をC_tとする.
・それらの配列が同一である (つまり突然変異がおこっていない) 場合,由来により同一 (identical by descent) という関係となる.
t世代前までの複製以降IBDである確率を,固定指数 (fixation index) といい,F_tと表記する.
・固定指数は,集団内にひとつの遺伝子の由来がどの程度広まっているかを表すと見ることもでき,集団構成員の近親度を表す指標として近親婚係数 (inbreeding coefficient) とも呼ばれる.


coalescence確率

集団の大きさをNとして,F_tを数式で表現する.簡便のため,突然変異はおこらないと考える.
まず,集団から無作為に選んだふたつの遺伝子が,1世代前にcoalescenceする確率は,親世代の2N個の遺伝子のうち特定のひとつに由来する確率なので,\frac{1}{2N}となる.
次に,1世代前ではcoalescenceせず,それ以前のt-1世代の間にcoalescenceする確率は,(1 - \frac{1}{2N}) F_{t-1}である.
よって,t世代前までのあいだにふたつの配列がcoalescenceする確率は,これらふたつの確率を足して,
F_t = \frac{1}{2N} + (1 - \frac{1}{2N}) F_{t-1}
となる.
これは漸化式なので,解くと,
F_t = 1 - (1 - \frac{1}{2N})^t [式8]
となる.
(初期世代ではふたつの遺伝子はcoalescenceしていないので,F_0 = 0)

式8は,少なくともt世代前までにcoalescenceがおこる確率 (1 - t世代前まで毎代coalescenceがおこらない確率) と見ることもできる.
式8からわかるように,F_tの増加は集団が小さいほど急である.


固定指数と集団分化指数の関係

式1と式8を比べると,集団内にひとつの由来の遺伝子が広まっている (F_tが小さい) ほど遺伝子頻度の分散の期待値が小さいことがわかる.
また,式3と式8を比べると,F_tは集団分化の指標であるF_{ST}と一致することがわかる.つまり,ふたつの遺伝子の由来が遠ければ遠いほど,同じだけ集団分化も進んでいることがわかる.

式5のF_{ST}F_Tに書き換えると,
\overline{H_S} = H_T (1 - F_t) [式9]
となる.分集団中の割合F_tの遺伝子は始原集団のある1分子からの複製に由来しているため,新たな変異が供給されなければ,ヘテロ接合度は0である.また,分集団中の割合1- F_tの遺伝子は始原集団中の互いに異なる分子に由来しおり,始原集団での遺伝子構成を今も残している (つまりヘテロ接合度はH_T) と考えることができる.したがって,分集団中のヘテロ接合体の割合は,式9のように,これら両方の積として得られる.


ところで,F_t = F_{ST}が成り立つのは,ひとつの始原集団から同じサイズの集団が同時に分化していくという理想状態のみであり,実際の系統関係はこのようにはおこらない.それでもF_{ST}は集団分化の指標として有用で,集団のヘテロ接合度が遺伝的浮動によって変化してく様子を近似するのに役立つ.
ただ,F_{ST}は無限の分集団間の全体について分化の程度を表すので,これをごく少数の集団間の遺伝的距離として使うのは不適切 (確率的なゆらぎのため,総合集団の遺伝子頻度が始原集団の遺伝子頻度と同じとは限らなくなる).


遺伝的浮動2 につづく

※1 分散の公式は,V(X) = E(X^2) - E(X)^2
※2 『分子進化遺伝学』根井正利著, 五條堀孝・斎藤成也訳. 培風館 1990年. 参照のこと.

集団遺伝学 (遺伝的変異と自然選択)

集団遺伝学の基礎について,授業の復習をメインにまとめてみます.
「適応進化遺伝学」テキストを参考にしています.
自分の学習用ですので,間違いなどあるかもしれません.


遺伝的変異

用語

遺伝的変異 (genetic variation): 集団内で遺伝子に多様性があること
多型 (polymorphism): 遺伝的変異をもつ状態
遺伝子頻度 (allele frequency): 集団中の対立遺伝子 (allele) の頻度
ヘテロ接合 (heterozygote): 2倍体生物の2つのalleleが異なる場合
遺伝子型 (genotype): 2つのalleleの組み合わせ

遺伝的多様性を記述する尺度

alleleの種類数 (number of alleles)
ヘテロ接合度 (heterozygosity): 狭義では,ヘテロ接合個体の割合.広義では,集団から無作為に選んだ2つの遺伝子が異なる種類のalleleである確率 (これを遺伝子多様度 (gene dyversity) ともいう).

Hardy-Weinberg principle

要点は,遺伝子頻度 (p, q) を与えれば,遺伝子型頻度の期待値 (p^2, 2pq, q^2) を計算できる,というところ.期待値と観察値のずれを,χ二乗検定やフィッシャーの直接確率検定で評価できる.
HWEは事実上1世代の遺伝子型頻度の (短い期間の) 変動を問題にしているため,実際には多くの集団でHWEが観察される.しかし,HWE成立条件のうち分集団構造だけは例外.集団構造の遺伝的分化が何世代にもわたって蓄積した複数分集団をひとつの交配集団とみなすと,任意交配からの大きな逸脱となる.このため,HWEは対象集団が均一な任意交配の集団化どうかの指標として使われる意味あいが強い.
分集団化のうち,生殖隔離がおこっていると,全ゲノムで遺伝子頻度が分化する.その一方,特定形質に対して異なる選択がはたらいていた場合は,特定の遺伝子型頻度のみが分化する.


自然選択による遺伝子頻度の変動

遺伝子淘汰 (genic selection)

野生型allele Aに対し,突然変異により生じたalleleをaとおく.遺伝子型AA, Aa, aaの淘汰値を,1, 1+s, 1+2sとして,世代tでの遺伝子aおよびAの頻度をそれぞれq_t, p_tとして,q_tの時間変化を計算してみると,
\Delta q_t = q_{t+1} - q_t = s p_t q_t
となる (sは非常に小さいので,s^2 = 0と近似できる).
これを微分方程式\frac{d q_t}{dt} = s p_t q_tにして解くと,
ln(\frac{q_t}{p_t}) = ln(\frac{q_0}{p_0}) + st
となり,
q_t = \frac{q_0}{q_0 + p_0 e^{-st}}
が得られる.
集団サイズが無限大であれば,変異型の比が時間とともに,指数的に増加してくことがわかる (allele頻度の対数が,時間に比例して増えていく).
これはつまり,有利な突然変異が現れれば,比較的短時間で既存遺伝子を凌駕し得ることを意味している.

超優勢選択 (overdominant selection)

ヘテロの状態が最高の適応度をもつ場合,遺伝子頻度は平衡に達する.遺伝子型AA, Aa, aaの淘汰値を1, 1+s, 1+t ( s > t \ge 0) として,Aの頻度をpaの頻度をqとおくと,
\Delta q = \frac{pq(tq + s - 2sq)}{1 + 2spq + tq^2}
となる.この式を微分方程式に直して解くと,qの時間変化が得られる.
遺伝的浮動の影響を無視すれば,初期頻度に関わらず一定時間の後に,平衡に達する.平衡に達したときのaの頻度は\Delta q = 0で与えられ,q = \frac{s}{2s - t}となる.