寻找算法来做长的成对的核苷酸比对

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)。似乎Needleman-Wunsch 相似矩阵也应该占用 3.1 gigs 的内存。不确定我是否没有正确记住 big-o 表示法，或者这实际上是在给定 R scripting 环境的开销的情况下构建对齐所需的内存开销。

是否有人建议使用更好的比对算法来比对更长的序列？已经使用 BLAST 完成了初始比对，以找到要比对的参考基因组区域。我对BLAST 正确放置插入缺失的可靠性并不完全有信心，而且我还没有找到与 biostrings 提供的用于解析原始 BLAST 对齐的 api 一样好的 api。

顺便说一下，这里有一个复制问题的代码sn-p：

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