긴 쌍 현자단 뉴클레오티드 배열을 수행하는 알고리즘을 찾고

Question

나는 참고 게놈의 하위 서열에 스캐 폴드를 정렬하여 가능한 SNPs 및 indels을 스캔하려고합니다. (원시 읽기는 사용할 수 없습니다). R/bioconductor과 Biostrings 패키지의 pairwiseAlignment 함수를 사용하고 있습니다. 나는 오류 메시지와 함께 같은 56kbp 발판을 정렬하려고 할 때 이 작은 발판을 위해 잘 작동하지만, 실패했습니다 :

오류 QualityScaledXStringSet.pairwiseAlignment (패턴 = 패턴, 에서 : 크기 17179869183.7 기가의 메모리 블록을 할당 할 수 없습니다 ?

나는 이것이 버그인지 아닌지 확실하지 않다, 나는 pairwiseAlignment에 의해 사용되는 Needleman-Wunsch algorithm 내가 3.1E9 작업 (56K * 56k ~= 3.1E9)의 순서로 계산 수요를 의미하는 것이라고 생각 O(n*m)이라는 인상이었다 그것은 또한 Needleman-Wunsch 유사도 행렬은 3.1 기가 메모리의 순서를 따라야합니다. 나는 제대로 큰 - 오 표기법을 기억하고 있지 않다거나 실제로 R scripting 환경의 오버 헤드가 주어진 정렬을 구축하는 데 필요한 것입니다 메모리 오버 헤드가 확실하지합니다.

누가 더 긴 시퀀스 정렬에 사용할 더 나은 정렬 알고리즘에 대한 제안 사항이 있습니까? 초기 정렬은 이미 BLAST를 사용하여 참조 게놈의 영역을 찾아서 정렬했습니다. 나는 제대로 삽입이나 삭제를 배치하기위한 BLAST의 신뢰성을 전적으로 확신하지 그리고 난 아직 원시 BLAST 정렬을 구문 분석 biostrings에서 제공하는 것과 같은 좋은 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 파일의 전체 디렉토리를 실행하고 정렬하는 스크립트입니다. 입/출력 요구에 맞게 시스템 호출에서 플래그를 수정할 수 있습니다. 클러스터 문서를 살펴보십시오. 이것은 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"; 
}