r - 寻找算法进行长对明智的核苷酸比对

我正在尝试通过将支架与参考基因组的子序列对齐来扫描可能的SNPs和indels。 (原始读取不可用)。我正在使用R/bioconductor和Biostrings包中的`pairwiseAlignment函数。
这对于较小的脚手架工作正常，但是当我尝试将56kbp脚手架与错误消息对齐时失败:

QualityScaledXStringSet.pairwiseAlignment(pattern = pattern，
:无法分配大小为17179869183.7 Gb的内存块

我不确定这是否是错误？ ;我的印象是Needleman-Wunsch algorithm使用的pairwiseAlignment是O(n*m)，我认为这将暗示计算需求大约为3.1E9操作(56K * 56k ~= 3.1E9)。看来Needleman-Wunsch相似性矩阵也应该占用3.1 gigs的内存。不知道我是否没有正确记住big-o表示法，或者在给定R scripting环境开销的情况下，实际上这是构建对齐所需的内存开销。

是否有人建议使用更好的比对算法来比对更长的序列？已经使用BLAST完成了初始比对，以找到要比对的参考基因组区域。我不完全相信BLAST正确放置indel的可靠性，而且我还无法找到与解析原始BLAST序列的生物字符串所提供的API一样好的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

较短的对齐方式(数百个碱基)不会发生该错误。
我尚未发现错误开始发生的长度截止值

最佳答案

因此，我将Clustal用作对齐工具。不确定具体的性能，但是在进行大量的大序列比对时却从未给我带来任何问题。这是一个脚本，它将运行.fasta文件的整个目录并对其进行对齐。您可以修改系统调用上的标志以适合您的输入/输出需求。只需查看相关的文档。这是在Perl中，我不会在对齐中使用太多R。您需要在脚本中编辑可执行路径，以匹配计算机上的提示位置。

#!/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";
}