AWK를 이용한 fasta file splitter

여러 염기서열로 이루어진 multi-fasta file을 쪼개에서 각 서열이 하나씩의 파일에 담기게 하는 작업은 매우 흔히 벌어지는 일이다. 나도 10줄 조금 넘는 Perl 스크립트를 만들어서 fasta file을 분리하는데 종종 쓰고는 한다.

awk를 쓰면 달랑 코드 한 줄로 fasta file splitter를 만들 수 있지 않을까? 구글에서 검색을 해 보니 이에 대한 많은 질문과 답이 있었다. 이를 참조하여 다음과 같은 awk 명령어를 만들어 보았다.  awk의 명령에 내에서 shell과 유사하게 파일로 직접 인쇄결과를 리다이렉션할 수 있다는 것을 몰랐다. 파일핸들을 정의하고 여닫는 일을 하지 않아도 되니 Perl보다 더 간단한 코드가 되었다. 아무런 실수가 없을 것으로...믿는다.

$ awk '/^>/{s=$0; sub(/^>/, "", $1); F=$1".fa"; print s > F; next}{print >> F}' input.fasta

awk를 쓰면서 종종 혼동하는 것은 sub() 함수가 치환된 문자열을 반환하지는 않는다는 것이다. 즉 var = sub(regexp, replacement, target) 형태로 사용하면 안된다.

코드를 간단히 분석해 보자. 서열 ID(>seq_id)가 있는 줄을 만나면, 첫번째 컬럼($0)에서 '>'를 제거한 뒤 .fa를 붙여서 출력할 파일의 이름을 만들어 변수 F에 저장한다. 그런데 sub() 함수에서 직접 $1을 건드리게 되므로, 나중에 원본 라인을 그대로 출력하기 위하여 명령어 맨 처음에 s라는 변수에 미리 서열 ID 라인을 저장해 놓는 꼼수(s=$0)를 사용했다. 그 다음에 유의할 것은 '/조건/{조건에 맞는 줄에 대하여; next} {조건에 맞지 않는 나머지 줄에 대해 실행}' 명령을 이해하는 것이다. 여기에서 next가 매우 중요한 역할을 담당한다. 이에 대해서는 예전의 포스팅에서 두번 정도 다룬 것 같다.

이 방법을 이용하면 SRR2162225(Gopalakrishnan et al., 2018)의 fastq 파일을 각 샘플별로 쪼갤 수도 있다. 이 fastq file은 다음과 같은 형식의 sequence ID를 갖고 있다. 내용상으로는 이미 demultiplexing이 된 것에 해당하지만 하나의 fastq 파일에 전부 수록되어 있다.

@ERR2162225.10 Wargo.109865_330283 length=253

두번째 필드를 추출하여 빨간색으로 표시된 것을 파일 이름으로 하는 fastq 파일을 만들어 보자. 현 디렉토리에 43개의 fastq 파일이 생겼다.

$ awk '/^@ERR2162225/{s=$0; sub(/_.+$/, "", $2); F=$2".fastq"; print s > F; next}{print >> F}' ERR2162225.fastq
$ ls
ERR2162225.fastq    Wargo.121296.fastq  Wargo.122637.fastq  Wargo.124418.fastq  Wargo.130511.fastq  Wargo.132266.fastq
Wargo.109865.fastq  Wargo.121409.fastq  Wargo.122639.fastq  Wargo.124419.fastq  Wargo.130801.fastq  Wargo.132268.fastq
Wargo.115690.fastq  Wargo.121585.fastq  Wargo.122724.fastq  Wargo.124421.fastq  Wargo.131051.fastq  Wargo.132538.fastq
Wargo.116774.fastq  Wargo.121891.fastq  Wargo.123009.fastq  Wargo.126747.fastq  Wargo.131300.fastq  Wargo.133270.fastq
Wargo.117349.fastq  Wargo.121892.fastq  Wargo.123197.fastq  Wargo.126802.fastq  Wargo.131302.fastq
Wargo.118369.fastq  Wargo.122114.fastq  Wargo.123305.fastq  Wargo.127261.fastq  Wargo.131512.fastq
Wargo.118370.fastq  Wargo.122410.fastq  Wargo.123704.fastq  Wargo.127262.fastq  Wargo.131541.fastq
Wargo.120662.fastq  Wargo.122412.fastq  Wargo.123986.fastq  Wargo.130510.fastq  Wargo.132265.fastq

다른 사례를 하나 더 들어보자. 다음은 vsearch --derep_fulllength로 dereplication을 거친 염기서열 파일의 일부를 보이고 있다.

>Mock.1;size=492
TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGTAGGCGGTTTTTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGAAAACTTGAGTGCAGAAGAGGA
AAGTGGAATTCCATGTGTAGCGGTGAAATGCGCAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACGCTGATGTGCGAAAGCGTGGGGATCAAACAGG
>Mock.2;size=345
TACAGAGGGTGCGAGCGTTAATCGGATTTACTGGGCGTAAAGCGTGCGTAGGCGGCTTATTAAGTCGGATGTGAAATCCCCGAGCTTAACTTGGGAATTGCATTCGATACTGGTGAGCTAGAGTATGGGAGAGGA
TGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGATGGCGAAGGCAGCCATCTGGCCTAATACTGACGCTGAGGTACGAAAGCATGGGGAGCAAACAGG

size=뒤에 나오는 모든 숫자를 더하면 dereplication 작업 전 FASTA file에 포함된 전체 서열의 수를 알 수 있다. 오늘 익힌 awk의 sub() 함수를 이용하면 간단하다.

$ awk 'BEGIN{s=0}/^>/{sub(/^>.*size=/, "", $0);s+=$0}END{print s}' input.fasta
4268

별로 설명할 필요도 없다. 서열 ID 라인에 description이 있다면 $0을 $1로 바꾸면 그만이다. 아니, 처음부터 code를 $1로 고쳐 놓으면 description이 존재하는지 여부에 상관없이 잘 돌아간다.

Bioinformatics on-liner

One-liner 애호가를 위하여 Stephen Turner의 GitHub 사이트를 소개한다(링크).