Rに塩基で表されるデータを読み込む時
普段読み込まないから、知らなかったがread.delimなどでデータの一列丸ごとが「T」一文字のデータを読み込むとlogicalに変換されてしまう(TRUE扱いになる)。
これを防ぐにはオプションでcolClassesを指定しておく(しないとRがcolClassを自動判定)。
data<-read.delim("hoge_typing.tsv", colClasses="factor")
Rで集計しようとして気がついた。
この後、昔作ってもらったHW平衡を計算するスクリプトが使えない状況になったのに気がついた(ToT)
作者さんは、明日お休みらしいので、良い機会だと思い込み、Rのgenetics packageに手を出す。
genetics packageはSNP向けの基本的な統計量(例えば、アレル頻度とかHW平衡の検定)を算出出来る。
例えば
hoge.csv Age,Sex,G1,G2 31,M,A/A,G/T 27,F,A/C,G/G 35,M,C/C,G/T 19,M,A/C,G/T 55,M,,G/G 34,F,A/A,G/G 45,F,A/C,T/T 32,M,A/C,G/T
Rを起動させて
> library(genetics) > gm1 <- makeGenotypes(read.csv("hoge.csv")) > gm1 Age Sex G1 G2 1 31 M A/A G/T 2 27 F A/C G/G 3 35 M C/C G/T 4 19 M A/C G/T 5 55 M <NA> G/G 6 34 F A/A G/G 7 45 F A/C T/T 8 32 M A/C G/T > HWE.test(gm1)
HWE.test(gm1)で計算結果がずらずら出てくる。
欲しいものだけValueを使って取り出すかHWE.exact()やsummary()を使って算出した方がよい感じだ。
library(help=genetics)で内容の一覧。
このgenetics packageの文章はhttp://www.amstat.org/chapters/connecticut/ASAConf/thirdMiniConf/ConferenceSlides/DrWarnes.pptを参考にしています。