일반적으로 FASTA 형태의 파일이라 하면 서열 데이터를 60 문자 단위로 끊어서 줄바꿈을 한다. 이를 wrapping이라 한다. Pacific Biosciences의 웹사이트에 나온 FASTA file format의 best practice에는 다음과 같은 항목이 있다.
Sequences in FASTA files should be wrapped at a uniform line length, to enable indexing. A common convention is to wrap lines at 60 characters.Wrapping이 인덱싱을 용이하게 하는 측면도 있었구나... 어쩐지 어떤 프로그램은 60 글자 단위를 지키지 않으면 에러를 토해 내더라니.
만약 FASTA 파일이 서열 ID 한 줄, 서열 자료는 줄바꿈 없이 한 줄로 구성되어 있다면 이를 다루는 후속 프로그램을 짜는 일도 무척 간편해진다. EZBioCloud에서 세균 유전체 파일을 받아서 열어보라. 바로 이런 unwrapped FASTA 형식으로 되어 있음을 알 수 있다. 아마도 데이터베이스 서버 내에 문자형 데이터로 저장되어 있어서 그런 것은 아닌지 모르겠다.
FASTA file의 wrapping 및 unwrapping을 할 수 있는 '생명정보학용' 유틸리티가 당연히 지구상 어딘가에 존재하겠지만(예: BBMap의 reformat.sh), Bash 환경과 이를 바탕으로 기본적으로 설치된 일반 유틸리티를 이용할 수 있다면 더욱 좋지 않겠는가? 오늘은 그 방법에 대해서 정리해 보려고 한다. 여기에서 소개하는 방법은 유일한 것이 아니다. 더 간단하고 유려한 방법이 있다면 그것을 쓰면 된다.
FASTA file unwrapping
줄바꿈이 적용된 FASTA file의 서열 데이터를 한 줄로 펼쳐보자. BioStars를 빛내는 Pierre Lindenbaum의 글에서 많은 영감을 받았다. N이라는 변수를 활용하는 트릭이 멋지다.
$ awk '/^>/ {printf("%s%s\n",(N>0?"\n":""),$0);N++;next;} {printf("%s",$0);} END {printf("\n");}' infile.fa
맨 마지막의 END{..} block을 빼먹으면 출력물을 파일로 저장하여 나중에 vi로 열었을 때 [noeol], 즉 이 파일이 개행문제로 끝나지 않았음을 표시하는 경고 문자가 뜬다. Unwrapped FASTA file은 하나의 서열이 정확하게 2 줄로 이루어지므로 파일을 쪼개는 등 다른 작업을 하기에 매우 편리하다.
FASTA to tab-delimited file
FASTA file을 unwrap하는 것에서 더욱 나아가서 아예 서열 ID를 비롯한 모른 자료를 한 줄로 표현해 보자. 이때 sequence description에 해당하는 것을 처리하는 것이 약간 성가신 일이 된다. 편의상 description 필드에 있는 첫번째 단어만을 취하여 하나의 컬럼으로 삼기로 하자. 즉 FASTA file로부터 서열 ID, description, 그리고 줄바꿈 없는 서열이라는 세 개의 컬럼이 탭으로 구분된 텍스트 파일을 만드는 것이다.
$ awk '/^>/ {printf("%s%s\t%s\t",(N>0?"\n":""),$1,$2);N++;next;} {printf("%s",$0);} END {printf("\n");}' infile.fa
그러나 이 방법은 아직 완벽하지 않다. 실행을 해 보면 알겠지만, 서열 ID 앞의 '>'가 그대로 출력되기 때문이다. 이를 해결하려면 sub() 함수를 써야 한다. sub() 함수 자체를 printf() 함수의 인자로 공급하면 제대로 작동이 안되므로, 해당되는 블록의 맨 앞에서 실행하여 패턴에 따르는 치환을 먼저 실시해야 한다. sub() 함수 내에 포함된 필드 번호는 $0과 $1 무엇을 쓰더라도 결과는 같다.
$ awk '/^>/ {sub(/^>/,"",$1); printf("%s%s\t%s\t",(N>0?"\n":""),$1,$2);N++;next;} {printf("%s",$0);} END {printf("\n");}' infile.fa
이번에는 서열의 일련번호를 첫번째 컬럼에 삽입해 보자. 즉 일련번호, 서열 ID, description 및 서열 데이터라는 네 개의 컬럼을 갖는 출력물을 만들어 내는 것이다. N+1을 인쇄하게 만든 이유는 일련번호를 1부터 시작하게 만들기 위함이다. 전산 전공자라면 0부터 시작하는 것을 더 좋아할 수도 있겠다.
$ awk '/^>/ {sub(/^>/,"",$1); printf("%s%d\t%s\t%s\t",(N>0?"\n":""),N+1,$1,$2);N++;next;} {printf("%s",$0);} END {printf("\n");}' infile.fa
FASTA file wrapping
서열 데이터가 한 줄로 펼쳐진 FASTA file을 다시 wrapping하려면 EMBOSS 패키지의 seqret을 쓰는 것이 편하다. awk로 처리하려면 상당히 귀찮다. 읽어서 그냥 출력하게 만들면 wrapping이 되어 나온다. 그러면 tab-delimted file로 풀어헤쳐진 것을 wrapped FASTA file로 만들어 보자. 다음과 같이 먼저 awk를 써서 unwrapped FASTA로 만든 다음 seqret에게 보내면 된다. EMBOSS의 매뉴얼은 이러한 용법을 쉽게 알 수 있도록 작성되지를 않았다는 것이 문제이다.
$ awk '{printf(">%s %s\n%s\n",$1,$2,$3)}' tab-delimited-file | seqret fasta::stdin fasta:outfile.fa
응용: Multi-FASTA description을 기준으로 특정 서열을 추출하기
FASTA file(infile.fa)의 description field에 gene 이름이 담겨 있다고 가정하자. gene list가 별도의 파일(list.txt)에 담겨 있으며, 이에 해당하는 서열을 FASTA file에서 추출하고 싶다. 어떻게 하면 좋을까? 우선 FASTA file을 tab-delimited file로 분해한 다음, gene 목록 파일과 join을 시킨 뒤 다시 FASTA로 전환하면 된다. join 전에는 공통 키를 기준으로 먼저 sort를 해 놓아야 한다. '>'를 떼었다 붙이는 번거로운 작업은 하지 않았다. list.txt 파일에 추출하고 싶은 서열의 ID가 들어있는 상태라면 '>'를 떼는 작업이 필요할 것이다. Process substitution을 사용하면 전체를 one-liner로 만들어 버릴 수 있지만 그렇게 되면 코드의 가독성이 떨어지니 오늘은 중간 단계 파일을 만들어서 다음 명령어에 인수로 공급하는 방법을 쓰기로 하겠다.
$ awk '/^>/ {printf("%s%s\t%s\t",(N>0?"\n":""),$1,$2);N++;next;} {printf("%s",$0);} END {printf("\n");}' infile.fa | sort -t $'\t' -k2,2 > infile.tsv.sorted $ sort -t $'\t' -k1,1 list.txt > list.txt.sorted $ join -t $'\t' -1 2 -2 1 -o 1.1,1.2,1.3 infile.tsv.sorted list.txt.sorted > extracted_sequences.tsv $ awk '{printf("%s %s\n%s\n", $1, $2, $3)}' extracted_sequences.tsv | seqret fasta::stdin fasta:outfile.fa
마지막에 소개한 응용편에서 약간 어렵지만 꽤 유용한 것은 join 명령어였다. 물론 list.txt에는 gene 정보가 있지만 infile.fa에는 없다면 어떻게 할 것인가? 이러한 예외적인 조건까지 다루도록 설계하지는 않았다. 하지만 옵션에 따라서 outer join도 가능하니 관심이 있다면 join 유틸리티의 사용법을 공부하면 된다.
댓글 없음:
댓글 쓰기