土と生き物

2020年8月31日月曜日

データが無い期間

午前
・標本整理

午後
・修論手伝い（データ解析）
・学会業務

修論のデータ解析。データが無い期間を0個として表現する方法でつまづいた、、、未だに、つまづいている。

これといった用事はなかったのに、午後は時折届くメールへの対応でバタバタした。

2020年8月28日金曜日

午前
・データ解析

午後
・データ解析

消耗品が切れてしまったので、一日中データ解析。数ヶ月待たせてしまったが、共同研究者にやっとdraft reportができた。

ここ数日苦戦したのが、haplotype networkの描画。結局、解決はしていないけど。

library(pegas) を使った方法だと、プロットのサイズの制御ができなかった。１枚だけなら良いのだが、今回は６つのサブグループがあり、グループ間でサイズの統一ができず断念。

library(haplotypes) では、プロットの配置の際に、シミュレーション？が行われているようで、出力するたびに微妙に配置が異なる。

完全な制御はできていないが、

coo1 <- plot(p, net=1, mode="fruchtermanreingold",,,,)
plot(p, net=1, coord=coo1, label.pos=5,,,,)

とすると、大きくズレることはない気がする？気のせいかも。

とはいえ、どうにかこうにか、haplotipe network図を作成し、

系統樹を作成し、

地図に投影できた。ゴチャゴチャ、、、でも、しばらくは考えたくない。一旦、お返しします。

2020年8月27日木曜日

どハマり

午前
・データ解析
・統計ゼミ

午後
・個人ゼミ
・データ解析

昨日で解決したと思ったHaplotype netwokの描画で、どハマり中、、、。

2020年8月26日水曜日

ハプロタイプ；シーケンス名変更；「.fasta」を出力

午前
・標本整理

午後
・授業準備
・データ解析

かなり苦戦したが、どうにかできるようになった。

library(haplotypes)を用いてハプロタイプごとにシーケンスをまとめ、さらに、シーケンス名を変更したのち、「.fasta」を出力する方法。

#library(haplotypes) #「::」で指定する
#library(ape) #「::」で指定する
x <- haplotypes::read.fas("???.fasta")
h <- haplotypes::haplotype(x) #to class"haplotype"
#h@nhap #ハプロタイプ数の確認
h_dna <- haplotypes::as.dna(h) #to class"DNA"
names(h_dna) <- c(1:34) #シーケンス名の変更；class"DNA"で使う
h_dnabin <- haplotypes::as.DNAbin(h_dna) #to class"DNAbin"
ape::write.dna(h_dnabin, file = 'test.fasta', format = 'fasta')

library()は使わず、「::」で指定した方が良い？詳しく検証していないが、２つのlibraryで同じ関数名が使われている？

関数が「class」に依存するので、順番を間違えるとハマる。「class:haplotype」→「class:DNA」→「class:DNAbin」で「ape::write.dna」が実行できる。

また、「names()」は「class:DNA」で使う。

2020年8月25日火曜日

packages"haplotypes"

午前
・標本整理

午後
・データ解析
・手伝い
・会議

変異が大きいため、pegas::haplotype、pages::haploNetを使った方法では、ハプロタイプネットワークがうまく描けなかった。

replot()を使って、頑張ればどうにかなりそうだけど、あまり現実できではない。

で、昨夜、packages"haplotypes"というのを見つけた。

ハプロタイプネットワークの老舗ソフトTCSでも利用されている最節約法によるサブネットワーク化？を計算することができる。pagesでもできるかも？

install.packages("haplotypes")
library(haplotypes)
#マニュアルでは「.fas」となっているが、うまくいかなかった。
x <- read.fas("***.fasta")
h <- haplotypes::haplotype(x)
p <- parsimnet(x, prob=.95)
rownames(p@d[[1]])<-c(1:8)
plot(p,
   net=1, #network番号の指定
   label.pos=5,
   vertex.cex=c(rep(5, p@nhap[[1]]), #円の大きさの指定
		rep(1, c(nrow(p@d[[1]])-p@nhap[[1]]))),
   vertex.col=c(rep("magenta", p@nhap[[1]]), #の色の指定
   		rep("deepskyblue3", c(nrow(p@d[[1]])-p@nhap[[1]]))))

パッケージのマニュアルでは、「p@nhap」となっているが、「p@nhap[[1]]」と「[[1]]」の指定が必要。この数字はネットワーク番号に合わせて変更する。

例えば、２つ目のネットワークの緑色で出力する場合は下記となる。

rownames(p@d[[2]])<-c(7:11)
plot(p,
   net=2,
   label.pos=5,
   vertex.cex=c(rep(5, p@nhap[[2]]),
		rep(1, c(nrow(p@d[[2]])-p@nhap[[2]]))),
   vertex.col=c(rep("green", p@nhap[[2]]),
   		rep("deepskyblue3", c(nrow(p@d[[2]])-p@nhap[[2]]))))