Perl에서 두 개의 FASTA 파일 (라인 브레이크가있는 파일)을 어떻게 병합합니까?
https://stackoverflow.com/questions/736855
-
09-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian문제
다음 두 파일이 있습니다.
file1.fasta
>0
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
>1
GTTAAGTTATATCAAACTAAATATACATACTATAAA
>2
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC
file2.qual
>0
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 15 40 40
>1
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40
40 40 40 40 40 40 40 40 40 40 40
>2
40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4
40 8 3 29 10 19 18 40 19 15 5
">"로 표시된 각 Fasta 헤더의 "Qual"파일의 라인 브레이크를 기록하십시오. 파일 헤더 수 ( '>')는 두 파일 모두 동일합니다. 숫자 품질의 수 = 시퀀스 길이.
내가하고 싶은 것은이 두 파일을 추가하는 것입니다.
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 15 40 40
GTTAAGTTATATCAAACTAAATATACATACTATAAA 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40 40 40 40 40 40 40 40 40 40 40 40
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC 40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4 40 8 3 29 10 19 18 40 19 15 5
그러나 어떻게 든 아래 코드가 올바르게 수행되지 않습니까? 특히 'Qual'파일의 각 항목의 두 번째 줄은 인쇄되지 않습니다.
use strict;
use Data::Dumper;
use Carp;
use File::Basename;
my $fastafile = $ARGV[0] || "reads/2039F.2.fasta";
my $base = basename( $fastafile, ".fasta" );
my $qualfile = "reads/" . $base . ".qual";
print "$qualfile\n";
open SEQ, '<', $fastafile or die $!; #Seq
open PRB, '<', $qualfile or die $!; #quality
while (my $seq = <SEQ>) {
my $qual = <PRB>;
chomp($seq);
chomp($qual);
if ($seq =~ /^>/ || $qual =~ /^>/) {
next;
}
else {
print "$seq\t$qual\n";
}
}
올바른 방법은 무엇입니까?
해결책
품질 점수의 2 위 (및 모든 후속) 라인이 누락되어 추가 시퀀스 라인도 누락됩니다. 이와 코드 재사용 목적으로 FASTA 시퀀스를 처리하는 방법은 전체 항목/레코드입니다.
local $/ = "\n>";
while (my $seq = <SEQ>) {
my $qual = <PRB>;
chomp($seq); $seq =~ s/^>*.+\n//; $seq =~ s/\n//g;
chomp($qual); $qual =~ s/^>*.+\n//; $qual =~ s/\n/ /g;
print "$seq\t$qual\n";
}
첫 번째 교체에서 Fasta 헤더를 쉽게 캡처 할 수도 있습니다.
다른 팁
문제는 파일을 통해 병렬로 발전하고 있으므로 한 줄이 ">"인 경우 다음 파일에서 ">"가 아닐 수도 있습니다.
데이터를 읽는 방식은 다음과 같이 쌍입니다.
1: >0 2: >0 1: GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT 2: 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 1: >1 2: 40 40 40 40 40 40 40 40 15 40 40 1: GTTAAGTTATATCAAACTAAATATACATACTATAAA 2: >1 1: >2 2: 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40 1: GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC 2: 40 40 40 40 40 40 40 40 40 40 40 1: EOF 2: >2 1: EOF 2: 40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4 1: EOF 2: 40 8 3 29 10 19 18 40 19 15 5
루핑 규칙을 적용한 동일한 데이터 세트가 다음과 같습니다.
1: GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT 2: 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 1: GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC 2: 40 40 40 40 40 40 40 40 40 40 40
따라서 루핑 로직을 분리하거나 파일을 일치시키는 방법을 찾아야합니다.
다음은 찾는 것을 분리하려는 시도가 있지만 테스트하지 않았습니다.
fileIO: {
while( 1 ){
my $seq;
my $qual = q{};
while( 1 ){
$seq = <SEQ>;
last fileIO if not $seq; # stop at end of file
last if $seq !~ /^>/;
}
while( 1 ){
my $qual_in = <PRB>;
last fileIO if not $qual_in; # stop at end of file
last if $qual_in =~ /^>/ and $qual ne q{};
next if $qual_in =~ /^>/ and $qual eq q{};
$qual .= $qual_in;
}
print "$seq \n $qual \n";
}
}
업데이트
위의 코드를 필요에 따라 임의의 파일 핸들에서 청크를 읽는 단일 함수로 다시 얻었으므로 필요에 따라 작동하는 것 같습니다. 물론 나는 실용적인 것을 위해 사용하는 트릭으로 여기에서 약간 실험했습니다.
use strict;
use warnings;
#
# readUntilNext( $fileHandle, \$scalar_ref );
#
# returns 0 when nothing could be read from the fileHandle.
# otherwise returns 1;
#
sub readUntilNext {
my ($fh) = shift;
my ($output) = shift;
my ($output_buffer) = '';
while (1) {
my $line = <$fh>;
if ( !$line ) { # No more data
# No data to flush to user, return false.
return 0 if $output_buffer eq q{};
last; # data to flush to user, loop exit.
}
if ( $line =~ /^>/ ) {
# Didn't get anything, keep looking.
next if $output_buffer eq q{};
# Got something, flush data to user.
last;
}
chomp($line);
$output_buffer .= $line;
}
# Data to flush to user
# Write to the scalar-reference
$$output .= $output_buffer;
return 1;
}
open my $m, '<', 'a.txt';
open my $n , '<', 'b.txt';
# Creates 2 scalar references every loop, and only loops as long
# as both files have data.
while ( readUntilNext( $m, \my $seq ) && readUntilNext( $n, \my $qual ) ) {
print "$seq\t$qual\n";
}
그리고 테스트 된 위의 코드는 정확히 당신이하고 싶은 일을합니다.
그 내 물건에 주목하십시오
while( readUntilNext( $m, \my $seq ) ) {
}
근본적으로 동일합니다
my $seq;
while( readUntilNext( $m, \$seq ) ) {
}
전자가 매번 새로운 스칼라를 생성한다는 사실을 제외하고는 동일한 값이 성공적인 루프로 보이지 않음을 보장합니다.
그래서 그것은 더 나옵니다.
while( 1 ){
my $seq;
last if not readUntilNext($m, \$seq);
do {
# loop body here
}
}
다음은 Perl을 사용하지 않지만 일반 쉘 명령을 사용하는 솔루션입니다.
prompt>grep -v '^>[0-9]' file1.fasta > tmp1
prompt>(tr '\012' ' ' < file2.qual; echo) | sed 's/>[0-9]* /\n/g' | sed 1d > tmp2
prompt>paste tmp1 tmp2
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 15 40 40
GTTAAGTTATATCAAACTAAATATACATACTATAAA 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40 40 40 40 40 40 40 40 40 40 40 40
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC 40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4 40 8 3 29 10 19 18 40 19 15 5
prompt>
나는 페이스트 명령을 수년 동안 검색했습니다 ( "이것은"이것은 매우 기본적인 작동입니다. 누군가 ~ 해야 하다 이미이 문제를 해결하기 위해 무언가를 구현했습니다. ").
두 번째 명령 줄은 먼저 모든 최신 라인을 공백으로 변환하고 ECHO 명령이 추가되어 입력에 최종 신자 라인을 추가하여 (SED가 EOL이없는 선을 무시하기 때문에) 모든 입력 라인을 하나의 단일 라인으로 연결 한 다음 SED 명령에 연결합니다. 다시 나뉘어진다 (이식성 참고 : 모든 SED 프로그램이 임의의 라인 길이와 함께 작동하는 것은 아니지만 GNU SED는 그렇습니다).
