2

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ユーザーであるため、この時点で完全に困惑しています。これに対する解決策をオンラインでたくさん調べましたが、いくつかを含まないことがわかったものはないようです上で説明したようにこれを行うのが難しい方法の変形ですが、おそらく格子は仕事を成し遂げることができますか?

4

0 に答える 0