이 주제는 나의 블로그에서도 여러 차례 다루어진 바가 있고, 구글을 검색해 보면 꽤 많은 질문과 답이 올라와 있다. 파이썬을 능숙하게 다루는 사람이라면 별로 어렵지 않게 스크립트를 짜서 원하는 바를 해결할 수 있을 것이다. 나 역시 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이 아니라 00001, 00002, 00003과 같이 자릿수를 꽉 채운 0('leading zero')이 있는 것으로 바꾸면 해결이 된다. 단, seq_id, description 및 서열의 3개 컬럼으로 이루어진 TSV 파일을 먼저 만든 뒤 nl 명령을 사용하여 첫 컬럼에 라인 번호를 넣는다. nl 명령은 라인 번호를 오른쪽으로 정렬한 뒤 leading zero를 넣는 옵션('-n rz')이 있다. 전체 명령어는 다음과 같다.
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))
|
