100 個の fasta ファイルがあり、遺伝的距離行列の重なり合うヒストグラムをプロットして、DNA データのブートストラップ複製間にどれだけの重なりがあるかを確認したいと思いますか?
次を使用して、各ファイルを ape に読み取らせる方法を見つけました。
files <- list.files("/Volumes/ALEX_R-HD/", pattern="/Volumes/ALEX_R-HD/xii-27-D")
library("ape")
library("pegas")
library("plyr")
library("dostats")
filenames <- dir(path="/Volumes/ALEX_R-HD/xii-27_D_coccus", full.names="TRUE", pattern="xii-27")
listOfiles <- lapply(filenames, function(x) read.dna(x, format="fasta"))
次に、次を使用して、それぞれの遺伝的距離行列を生成します。
distOfiles <- lapply(listOfiles, function(y) dist.dna(y, model="TN93"))
R コンソールから呼び出すと、遺伝的距離ファイルは次のようになります。
[[1]]
M_51_1_new__ M_51_3_new__ M_51_4_new2__ M_51_5_new2__ M_51_6_new__ M_51_7_new__ M_51_8_new__ madera_1_new__ madera_2_new__ madera_3__ madera_4_new__ madera_5_new__
M_51_3_new__ 0.000000000
M_51_4_new2__ 0.000000000 0.000000000
M_51_5_new2__ 0.000000000 0.000000000 0.000000000
M_51_6_new__ 0.000000000 0.000000000 0.000000000 0.000000000
M_51_7_new__ 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
M_51_8_new__ 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
madera_1_new__ 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343
madera_2_new__ 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.000000000
madera_3__ 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.037124343 0.000000000
and goes on to.... [[100]]
私が問題に遭遇するのは、各ブートストラップが同じウィンドウ内で他のブートストラップの上にプロットされるように、それぞれのヒストグラムをプロットすることです。以下のスクリプトは、それぞれを新しいウィンドウにプロットするだけで、重複しません。
bins=seq(0,0.05,by=0.001)
HistOfiles <- lapply(distOfiles, function(z) hist(z, breaks=bins, main="Histogram of D. coccus Mexico-types TN93 distances", ylim=c(0,1500), xlab="TN93 distance", ylab="frequency", col=rgb(0,0,0,0.01), border=rgb(0,0,0,0.01)))
これは、次の方法で難しい方法で実行できることを知っています。
bins=seq(0,0.05,by=0.001)
readfile1 <- read.dna("/Volumes/ALEX_R-HD/xii-27_D_coccus/xii-27-D_coccus1", format="fasta")
distance_TN931 <- dist.dna(readfile1, model="TN93")
bins=seq(0,0.05,by=0.001)
hist(distance_TN931, breaks=bins, main="Histogram of D. coccus Mexico-types TN93 distances", ylim=c(0,1500), xlab="TN93 distance", ylab="frequency", col=rgb(0,0,0,0.01), border=rgb(0,0,0,0.01))
lines(density(distance_TN931), col=rgb(1,0,0,0.01))
par(new=TRUE)
readfile2 <- read.dna("/Volumes/ALEX_R-HD/xii-27_D_coccus/xii-27-D_coccus2", format="fasta")
distance_TN932 <- dist.dna(readfile2, model="TN93")
bins=seq(0,0.05,by=0.001)
hist(distance_TN932, breaks=bins, ylim="", main="", xlab="", ylab="", col=rgb(0,0,0,0.01), border=rgb(0,0,0,0.01))
lines(density(distance_TN932), col=rgb(1,0,0,0.01))
par(new=TRUE)
.......最後のファイルへ
しかし、それは大変な作業になると思います。100 個のファイルには問題ありませんが、1,000 個のファイルを持っている他の人 (たとえば、GenBank データで作業している人など) の場合、これは多すぎるかもしれません。
また、いくつかの Unix を使用して別の方法で別のファイルを \t で区切られた列のリストに貼り付けることも試みました。
paste /Volumes/ALEX_R-HD/xii-27_D_coccus/xii-27-D_coccus* /Volumes/ALEX_R-HD/xii-27_D_coccus/blank > /Volumes/ALEX_R-HD/xii-27_D_coccus/blank
そのファイルは次のように表示されます。ファイルがどのように分離されているかを明確にするために、"" \t
>name1 "\t" >name1 "\t" >name1 ...... 100 times for each row
actgactg "\t" actgaca "\t" actgaca
actgttgc "\t" actgact "\t" actgaca
>name2 "\t" >name2 "\t" >name2
actgactg "\t" actgaca "\t" actgaca
actgttgc "\t" actgact "\t" actgaca
しかし、read.dna を取得して各列を個別のデータ マトリックスとして読み取る方法がわかりません。read.table を取得してファイルを読み取ることができますが、そこでスタックしてしまいます
私は新しいRユーザーであるため、この時点で完全に困惑しています。これに対する解決策をオンラインでたくさん調べましたが、いくつかを含まないことがわかったものはないようです上で説明したようにこれを行うのが難しい方法の変形ですが、おそらく格子は仕事を成し遂げることができますか?