尋找算法做長對核苷酸比對

Question

我試圖掃描可能的SNPs和indels通過將支架與參考基因組的子序列進行比對。 （原始讀取不可用）。我正在使用R/bioconductor和Biostrings軟件包中的`pairwiseAlignment函數。 這是工作的罰款更小的支架，但未能當我試圖與錯誤消息作爲對準支架56kbp：

錯誤QualityScaledXStringSet.pairwiseAlignment（圖案=圖案， ：無法分配大小17179869183.7 GB的內存塊

我不知道這是否是一個錯誤或沒有;我的印象是，通過pairwiseAlignment使用的Needleman-Wunsch algorithm是O(n*m)我本以爲這意味着計算需求是3.1E9操作(56K * 56k ~= 3.1E9)的順序它本身相似性矩陣也應該佔用3.1G內存的順序。不確定是否我沒有正確記住big-o符號，或者實際上是在給定R scripting環境的開銷時構建對齊所需的內存開銷。

有沒有人有更好的比對算法用於比對長序列的建議？使用BLAST已經完成了初始比對以找到參比基因組的區域以進​​行比對。我並不完全有信心BLAST的正確放置indels的可靠性，我還沒有找到一個好像biostrings提供的解析原始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 

的錯誤不會用較短的比對（幾百個鹼基）發生。 我還沒有找到錯誤開始發生的長度截止點

Answer 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"; 
}