r - 寻找算法来进行长配对核苷酸比对

Question

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

Error in QualityScaledXStringSet.pairwiseAlignment(pattern = pattern, : cannot allocate memory block of size 17179869183.7 Gb

我不确定这是不是一个错误？ ;我的印象是 pairwiseAlignment 使用的 Needleman-Wunsch 算法 是一个 O(n*m)，我认为这意味着计算要求按3.1E9操作顺序(56K * 56k ~= 3.1E9)。看来 Needleman-Wunsch 相似度矩阵也应该占用 3.1 g 的内存量级。不确定我是否没有正确记住 big-o 表示法，或者这实际上是在给定 R 脚本 环境开销的情况下构建对齐所需的内存开销。

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

较短的比对(数百个碱基)不会发生错误。 我还没有找到错误开始发生的长度截止点

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