3

我正在尝试扫描可能的SNPs并将indels支架与参考基因组的子序列对齐。（原始读数不可用）。我正在使用R/bioconductorBiostrings 包中的 `pairwiseAlignment 函数。这对于较小的脚手架工作得很好，但是当我尝试将 56kbp 脚手架与错误消息对齐时失败了：

QualityScaledXStringSet.pairwiseAlignment 中的错误（模式 = 模式，：无法分配大小为 17179869183.7 Gb 的内存块

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

有没有人建议使用更好的对齐算法来对齐更长的序列？已经使用 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

较短的对齐（数百个碱基）不会发生该错误。我还没有找到错误开始发生的长度截止

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 回答