3

私は可能性をスキャンしようとしてSNPsおりindels、足場を参照ゲノムのサブシーケンスに整列させています。(生の読み取りは利用できません)。R/bioconductorBiostrings パッケージの「pairwiseAlignment」関数を使用しています。これは小さなスキャフォールドでは問題なく機能していましたが、56kbp スキャフォールドとして整列しようとすると失敗し、次のエラー メッセージが表示されました。

QualityScaledXStringSet.pairwiseAlignment のエラー (パターン = パターン、: サイズ 17179869183.7 Gb のメモリ ブロックを割り当てることができません

これがバグかどうかわかりませんか？; Needleman-Wunsch algorithm私は、 used bypairwiseAlignmentは、計算の要求がoperationsO(n*m)の順序にある​​ことを意味すると私が考えたであるという印象を受けました。同様に、類似度マトリックスも 3.1 ギガのメモリを占めるようです。big-o表記を正しく覚えていないのか、それとも環境のオーバーヘッドを考慮してアライメントを構築するために実際に必要なメモリオーバーヘッドなのかはわかりません。3.1E9(56K * 56k ~= 3.1E9)Needleman-WunschR scripting

より長いシーケンスのアラインメントに使用するより良いアラインメント アルゴリズムの提案はありますか? 最初のアラインメントは、BLAST を使用してすでに行われており、アラインメントする参照ゲノムの領域を見つけています。インデルを正しく配置するための の信頼性に完全に自信があるわけでBLASTはなく、生の BLAST アラインメントを解析するためのバイオストリングによって提供されるものと同じくらい良い API をまだ見つけることができませんでした。

ところで、問題を再現するコード スニペットを次に示します。

library("Biostrings")
scaffold_set = read.DNAStringSet(scaffold_file_name) #scaffold_set is a DNAStringSet instance
scafseq = scaffold_set[[scaffold_name]] #scaf_seq is a "DNAString" instance

genome = read.DNAStringSet(genome_file_name)[[1]] #genome is a "DNAString" instance

#qstart, qend, substart, subend are all from intial BLAST alignment step
scaf_sub = subseq(scafseq, start=qstart, end=qend) #56170-letter "DNAString" instance
genomic_sub = subseq(genome, start=substart, end=subend) #56168-letter "DNAString" instance

curalign = pairwiseAlignment(pattern = scaf_sub, subject = genomic_sub)
#that last line gives the error: 
#Error in .Call2("XStringSet_align_pairwiseAlignment", pattern, subject,  : 
#cannot allocate memory block of size 17179869182.9 Gb

このエラーは、短いアラインメント (数百の塩基) では発生しません。エラーが発生し始める長さのカットオフをまだ見つけていません

4

1 に答える 1

-1

なのでClustalをアライメントツールとして使っています。特定のパフォーマンスについてはわかりませんが、大量の複数の配列アラインメントを行うときに問題が発生したことはありません. 以下は、.fasta ファイルのディレクトリ全体を実行して整列させるスクリプトです。入出力のニーズに合わせて、システム コールのフラグを変更できます。clustalのドキュメントを見てください。これは Perl で作成されたもので、アラインメントに R をあまり使用しません。スクリプト内の実行可能パスを編集して、コンピューター上の clustal の場所と一致させる必要があります。

#!/usr/bin/perl 


use warnings;

print "Please type the list file name of protein fasta files to align (end the directory    path with a / or this will fail!): ";
$directory = <STDIN>;
chomp $directory;

opendir (DIR,$directory) or die $!;

my @file = readdir DIR;
closedir DIR;

my $add="_align.fasta";

foreach $file (@file) {
 my $infile = "$directory$file";
 (my $fileprefix = $infile) =~ s/\.[^.]+$//;
 my $outfile="$fileprefix$add";
 system "/Users/Wes/Desktop/eggNOG_files/clustalw-2.1-macosx/clustalw2 -INFILE=$infile -OUTFILE=$outfile -OUTPUT=FASTA -tree";
}
于 2012-10-09T05:43:11.927 に答える