【R】Hardy-Weinberg平衡（Fisherの正確確率検定）

はじめに

仕事で統計遺伝学が必要になりそうなので、実践でわかる!Rによる統計遺伝学を始めた。

備忘録的に記していく。

使用するデータ

使用するデータは本で指定されたURLから。

FAMussデータ
http://pub.maruzen.co.jp/book_magazine/support/r_toukeiidengaku/FMS_data.txt

必要なパッケージ

geneticsパッケージをインストールしておく。

install.packages("genetics")

そして読み込む。

library(genetics)

データの読み込み

FAMuSSデータを読み込む。
**read.delim()**関数は、データ区切りがタブの場合。

参考
ファイルからデータを読み込むとき

fms <- read.delim("http://pub.maruzen.co.jp/book_magazine/support/r_toukeiidengaku/FMS_data.txt",fileEncoding="cp932")

**attach()**関数を使うことで、データフレーム名を使わず変数名だけを用いて変数を指定する。

attach(fms)

Fisherの正確確率検定を用いたHardy-Weinberg平衡の検定

geneticsパッケージの**HWE.exact()**関数を使って求める。

5未満の期待度数が全セルの20%以上ある場合はピアソンのχ二乗検定を使わずにFisherの正確確率検定を使う。

民族がAfrican Amだけを用いる。

Akt1SNPAfricanAm <- genotype(akt1_t10726c_t12868c[Race == "African Am"], sep = "")
HWE.exact(Akt1SNPAfricanAm)

実行結果。

> Akt1SNPAfricanAm <- genotype(akt1_t10726c_t12868c[Race == "African Am"], sep = "")
> HWE.exact(Akt1SNPAfricanAm)

	Exact Test for Hardy-Weinberg Equilibrium

data:  Akt1SNPAfricanAm
N11 = 18, N12 = 18, N22 = 1, N1 = 54, N2 = 20, p-value = 0.2353

この結果からp値は0.2353で、この集団では帰無仮説のHardy-Weinberg平衡を棄却できない。

最終的なスクリプト

# Rをきれいにする
rm(list = ls())

# データの読み込み
fms <- read.delim("http://pub.maruzen.co.jp/book_magazine/support/r_toukeiidengaku/FMS_data.txt",fileEncoding="cp932")

# ライブラリの読み込み
library(genetics)

# attach()関数でデータフレームを指定
attach(fms)

# Fisherの正確確率検定を用いたHardy-Weinberg平衡の検定
# 変数Raceの水準の確認
table(Race)

Akt1SNPAfricanAm <- genotype(akt1_t10726c_t12868c[Race == "African Am"], sep = "")
HWE.exact(Akt1SNPAfricanAm)