【R】Hardy-Weinberg平衡（ピアソンのχ二乗検定）

はじめに

仕事で統計遺伝学が必要になりそうなので、実践でわかる!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)

ピアソンのχ二乗検定を用いたHardy-Weinberg平衡の検定

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

FAMussデータのakt1_t10726c_t12868cSNPについて、ピアソンのχ二乗検定を用いてHardy-Weinberg平衡の検定を行う。

Akt1SNP <- genotype(akt1_t10726c_t12868c, sep = "")
HWE.chisq(Akt1SNP)

実行結果。

> Akt1SNP <- genotype(akt1_t10726c_t12868c, sep = "")
> HWE.chisq(Akt1SNP)

	Pearson's Chi-squared test with simulated p-value (based on 10000 replicates)

data:  tab
X-squared = 26.539, df = NA, p-value = 9.999e-05

Hardy-Weinberg平衡に関する検定結果が統計的に有意なら、そのSNPは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)

# ピアソンのχ二乗検定を用いたHardy-Weinberg平衡の検定
Akt1SNP <- genotype(akt1_t10726c_t12868c, sep = "")
HWE.chisq(Akt1SNP)