이 주제는 나의 블로그에서도 여러 차례 다루어진 바가 있고, 구글을 검색해 보면 꽤 많은 질문과 답이 올라와 있다. 파이썬을 능숙하게 다루는 사람이라면 별로 어렵지 않게 스크립트를 짜서 원하는 바를 해결할 수 있을 것이다. 나 역시 Perl을 이용하여 이따금씩 필요한 염기서열을 multi-FASTA file에서 추출하고는 한다.
KentUtils라 불리는 Jim Kent의 유틸리티 중에서 faSomeRecords(다운로드 링크)가 이러한 목적에 딱 맞는 '이미 잘 알려진' 프로그램이다.'-exclude' 옵션을 이용하면 목록 파일(listFile)에 있는 서열을 제거하는 정반대의 동작을 한다.
$ faSomeRecords
faSomeRecords - Extract multiple fa records
usage:
faSomeRecords in.fa listFile out.fa
options:
-exclude - output sequences not in the list file.
BBMap에 포함된 유틸리티 중에서 Filterbyname.sh 또한 faSomeRecords와 거의 같은 동작을 한다. 단, fastq 파일도 다룰 수 있다는 것은 장점이다.
# By default, "filterbyname" discards reads with names in your name list, and keeps the rest. To include them and discard the others, do this: $ filterbyname.sh in=003.fastq out=filter003.fq names=names003.txt include=t
잘 알려진 프로그램을 사용하는 방법을 알아 보았으니 이제부터는 리눅스에 기본적으로 포함된 유틸리티(awk 등)을 사용하여 다소 원시적인 방법으로 똑같은 일을 해 보자. 단, 전제 조건을 살짝 비틀어 보았다. 추출할 서열의 ID가 아니라 일련번호(1부터 시작)를 별도의 목록 파일에 갖고 있다고 가정한다.
FASTA file의 서열 ID, description 및 염기서열(줄바꿈을 제거하여 한 줄로 펼쳐야 함 - 흔히 unwrapping이라 부른다)을 탭으로 구분된 TSV 파일을 만들고, 이를 목록 파일과 비교하여 공통된 컬럼을 갖는 라인을 join으로 병합한 뒤 해당 라인만 추출하여 다시 FASTA file로 전환하는 것이 매우 일반적인 방법이다. 이러한 방법을 쓰는 경우에는 꼭 서열 일련번호가 아니라 서열 ID를 써서 원하는 서열을 뽑는 것도 무난하게 할 수 있다.
최근 알게 된 long read용 binning program인 LRBinner의 출력 파일('Bins.txt') 정보를 이용하여 입력 FASTA file에서 각 bin에 해당하는 염기서열을 별도의 파일로 뽑아내는 방법을 알아보도록 한다. Bins.txt는 입력 파일에 실린 각 염기서열이 속하는 bin 번호를 다음과 같이 보여준다. 참 건조하고 재미 없는 포맷이다.
2 1 3 2 0 ...
2번 bin에 속하는 염기서열은 1번과 4번에 해당한다. 2라는 값을 갖는 라인의 번호를 뽑아서 LineNum.txt 파일에 저장하려면 다음과 같이 awk 명령어를 실행한다.
awk '$1==2{print NR}' Bins.txt > LineNum.txt
nr이라는 유틸리티는 인수로 주어진 텍스트 파일을 라인 번호와 함께 출력한다. 당장 아무 텍스트 파일이나 가져다가 'nr file.txt'를 실행해 보라. 행 시작 부분의 공백을 제거하고 구분자를 탭으로 전환하려면 다음과 같이 약간의 수정을 거쳐야 한다. 맨 왼쪽에 라인 번호를 삽입한 텍스트 파일은 앞으로도 종종 쓸모가 있을 것이다.
nl -n ln -s $'\t' Bins.txt > BinWithLineNum.txt
이상에서 설명한 방법으로 추출해야 할 염기서열의 번호를 별도의 파일('LineNum.txt')에 저장하였다고 하자. 그러면 우선 입력물인 input.fa 파일을 TSV 파일('SeqWithNum.tab')로 펼치되, 첫 번째 컬럼에는 1에서부터 시작하는 서열 일련번호를 넣기로 한다.
awk '/^>/ { id=$1; sub(/^>/,"",id); $1=""; printf("%s%d\t%s\t%s\t",(N>0?"\n":""),N+1,id,$0); N++; next; } { printf("%s", $0); } END {printf("\n");}' input.fa > SeqWithNum.tab
awk 명령이 좋은 점은 따옴표로 둘러싼 긴 명령어 내부에서 가독성을 위해 줄바꿈을 해도 된다는 것. 단, END와 {..} 사이에 줄바꿈을 넣으면 에러가 남을 확인하였다. BEGIN과 {..} 사이도 마찬가지일 것이라고 생각한다. 이 awk 명령어를 여기에서 설명하지는 않겠다. 첫 번째 컬럼에 일련번호를 삽입하고, descripton에 해당하는 여러 단어를 하나의 컬럼에 들어가게 하는 등 꽤 많은 공을 들였다. 하지만 이 awk 명령어는 나의 창작은 아니다. 또한 이것이 awk를 이용하여 FASTA file을 펼치는 유일하거나 가장 능률적인(짧은?) 방법도 아니다. 가끔씩 이 awk 명령어를 바라보면서 왜 이렇게 짜여졌는지 생각을 하는 것도 좋은 공부가 될 것이다.
LineNum.txt와 SeqWithNum.tab이라는 두 개의 파일이 만들어졌다. LineNum.txt 파일에 수록된 라인 번호에 해당하는 것을 SeqWithNum.tab에서 골라낸 뒤, 이를 FASTA file로 전환하는 것이 오늘의 목표이다. SeqWithNum.tab 파일의 첫 번째 컬럼에 1로부터 시작하는 라인 번호가 들어 있어서 점검하기는 편하다.
오늘의 최종 목표를 달성하는 방법도 몇 가지가 있을 것이다. 우선 첫 번째 방법을 소개한다. 사실 이 방법을 쓰는 경우라면 SeqWithNum.tab 파일에 라인 번호를 별도의 컬럼으로 넣지 않아도 된다. 자세한 설명은 생략한다.
awk 'FNR==NR{a[$1];next}(FNR in a){print}' LineNum.txt SeqWithNum.tab
표준 출력으로 SeqWithNum.tab에서 조건에 맞는(즉 LineNum.txt에서 라인 번호를 지정한) 라인이 나올 것이다. 따라서 이를 FASTA file로 전환하여 저장하려면 파이프('|') 기호 뒤에 다음 명령어를 넣어야 한다(A).
awk -F"\t" '{printf(">%s %s\n%s\n",$2,$3,$4)}' | seqret fasta::stdin fasta:outfile.fa
왼쪽의 awk 명령어만을 사용해도 FASTA 기본 골격으로 출력은 된다. 그러나 서열 영역은 한 줄로 펼쳐진 형태가 되므로 이를 EMBOSS의 seqret로 처리하여 wrapping을 실시하였다.
최종 목표를 달성하는 두 번째 방법을 소개한다. 그것은 바로 sort와 join 명령어를 이용하는 것인데, 검색을 해 보면 두 텍스트 파일에서 특정 컬럼을 서로 비교하여 같은 값을 갖는 라인(혹은 선별된 컬럼)을 추출하는 일반적인 솔루션으로 소개하고 있다. 예를 들어서 두 개의 파일 file_A와 file_B가 있고, 각각의 첫 번째 컬럼이 같은 값을 갖는 경우 file_A의 1, 2번째 컬럼과 file_B의 4, 5번째 컬럼을 추출하 싶다면 다음과 같이 해야 한다. 컬럼은 탭 문자로 구분되어 있다고 가정한다.
join -t $'\t' -1 1 -2 1 -o 1.1,1.2,2.3,2.5 <(sort -t $'\t' -k1,1 file_A) <(sort -t $'\t' -k1,1 file_B)
생각보다 꽤 복잡하다. 리눅스 shell에서 컬럼 구분자를 탭문자로 설정하는 것(-t $'\t')도 그렇고, file_A와 file_B가 join field(같은 값을 갖는지 비교할 컬럼)에 대해서 이미 sort가 되어 있다 해도 여기에서 보인 process substitution에 의해서 join 명령어 내에서 sort를 다시 하지 않으면 경고문을 발생한다. '--nocheck-order' 옵션을 주면 될 것 같으나 제대로 작동을 하지 않아서 몹시 마음에 들지 않는다. sort나 join 모두 사용법을 정확히 알지 않으면 낭패를 겪기 쉬운 명령어이다.
따라서 LineNum.txt와 SeqWithNum.tab의 두 파일이 준비된 경우에 각 파일의 첫 번째 컬럼을 비교하여 같은 값을 갖는 라인을 SeqWithNum.tab에서 추출하여 FASTA file로 되돌리는 '두 번째의 복잡한' 방법은 다음과 같다.
join -t $'\t' -1 1 -2 1 <(sort -t $'\t' -k1,1 LineNum.txt) <(sort -t $'\t' -k1,1 SeqWithNum.tab)
이 명령어 뒤에 파이프 기호를 쓴 뒤 위에서 소개한 (A) 명령어를 써 넣어야 FASTA file로 결과가 저장된다. 이 방법은 아직 완벽하지 않다. process substitution 내에서 일반 문자의 순서로 소팅했기 때문에 라인 넘버로는 2, 3, 15 순서가 아니라 15, 2, 3...의 순서로 출력된다. sort 명령어의 단독 실행에서는 -k1n,1 또는 -nk1,1 옵션을 제공하여 숫자 기준의 정렬이 되지만, 이상하게도 join 명령어 안에서 <(sort ...) 형태로 실행하면 sort가 되지 않았다는 오류 메시지가 나온다. 따라서 단순하게 -k1,1 또는 -k1 옵션을 주어서 문자 순서로 sort를 해야만 한다.
오늘 학습의 교훈:
- Don't reinvent the wheel. 창의력을 시험하고 싶다면 바퀴를 다시 발명해도 좋다.
- awk에서 가끔은 printf() 함수를 사용할 필요가 있다.
- sort 명령의 -k(key)는 시작과 끝을 지정하는 키 2개를 지정하는 것이 기본이 되어야 한다.
- 예전에 쓴 글 Awk를 이용한 간단한 텍스트 파일의 조작(join)을 가끔 들추어 보면서 awk 명령어의 기본을 잊지 말도록 하자.
Numeric sort 문제가 없는 join 방법
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | # 3개 컬럼으로 구성된 파일을 만든다. $ awk '/^>/ { id=$1; sub(/^>/,"",id); $1=""; # 아래 코드는 seq_id, description, sequence의 3개 컬럼을 만든다, printf("%s%s\t%s\t",(N>0?"\n":""),id,$0); # 아래 코드는 1로부터 시작하는 번호를 첫 번째 컬럼에 넣는다(총 4개 컬럼). # printf("%s%d\t%s\t%s\t",(N>0?"\n":""),N+1,id,$0); N++; next; } { printf("%s", $0); } END {printf("\n");}' input.fa > Seq.tab # 라인 번호를 첫 컬럼으로 삽입한다. $ nl -n rz -s $'\t' Seq.tab > SeqWithNumLeadingZeros.tab # LineNum.txt의 숫자에 leading zero를 채운다. # 그러려면 SeqWithNumLeadingZeros.tab에서 자릿수가 얼마나 필요한지를 알아내야 한다. $ for n in $(cat LineNum.txt) > do > printf "%06d\n" $n >> LineNumLeadingZeros.txt > done # join 실행(sort가 필요하지 않음) $ join -t $'\t' LineNumLeadingZeros.txt SeqWithNumLeadingZeros.tab |\ awk -F"\t" '{printf(">%s %s\n%s\n",$2,$3,$4)}' |\ seqret -auto fasta::stdin fasta:outfile.fa |
SeqKit의 fx2tab이라는 명령어를 쓰면 FASTA file을 TSV file로 전환할 수 있다. 물론 내가 원하는 형태와는 약간 다르다.
어쩌면 나는 바퀴를 새로 발명하려고 애를 쓰고 있었는지도 모른다...
How to convert fasta file to tab delimited file
1 2 3 4 | from Bio import SeqIO for record in SeqIO.parse('input.fa', 'fasta'): print('{}\t{}\t{}'.format(record.id, record.description, record.seq)) |
댓글 없음:
댓글 쓰기