2020年1月14日火曜日

Rでハプロタイプネットワーク

午前
・データ整理
・個人ゼミ

午後
・来客
・データ整理
・特別業務
・データ整理

共同研究者をかなり待たせてしまっているサソリモドキ論文を進める。

Rでハプロタイプネットワークを書く。ここを参照しました。

libraryがあるので作業は非常に簡単。1点悩んだのが、デフォルトでは円の中に「I、II、III」というローマ数字が入る点。

私の場合、円の中には独自のハプロタイプ番号を付けたい。plothaplotype(**, labels=)で指定できるが、デフォルトで算出される「I、II、III」がどのように決まるのか良く分からず悩んだ。

結局、ここの「ind.hap」を使って確認して地道に書いた。

install.packages("ape", repos="http://cran.ism.ac.jp/")
install.packages("pegas", repos="http://cran.ism.ac.jp/")
setwd("***")
library(pegas)
library(ape)
COI<-read.dna("***.fasta", format="fasta") #FASTAデータの読み込み
COIHaps <- haplotype(COI) #シーケンスデータからハプロタイプの検出と整理
COINet <- haploNet(COIHaps) #シーケンスデータからハプロタイプの検出と整理
ind.hap<-with(
stack(setNames(attr(COIHaps, "index"), rownames(COIHaps))),
table(hap=ind, individuals=rownames(COI)[values])
)
ind.hap #デフォルトで表示されるハプロタイプ番号を確認する。下記の「I、II、III」が円の中に記される


hap_ID <- c("St57","St58","St59") #I=St57、II=St58、III=St59にする
COIHaps <- haplotype(COI, labels=hap_ID) #ラベルを追加して、シーケンスデータからハプロタイプの検出と整理のやり直し
COINet <- haploNet(COIHaps) #シーケンスデータからハプロタイプの検出と整理
plot(COINet, size = attr(COINet, "freq"), fast = T, col=2, bg=2, scale.ratio = 1, cex = 1)

ハプロタイプ番号を挿入する点のみ少し面倒だけど、とても簡単にネットワーク樹が作成できます。アルゴリズムによって、いくつか関数があるみたい。


円の大きさの制御など、細かな点は勉強が必要だけど、plot()で出力しているので、そこまで難しくはないハズ。

12月にはボロボロに負けた早稲田。新国立で令和初のチャンピオンに!