Barrnap 결과물(GFF3)에서 16S rRNA 서열만 뽑아내기

Barrnap(BAsic Rapid Ribosomal RNA Predictor, 링크)를 사용하여 박테리아의 유전체 서열로부터 ribosomal RNA를 예측하였다. 결과물은 GFF3 파일이다. 여기에서 16S rRNA sequence를 뽑아내려면 어떻게 하면 좋을까? 막상 Perl을 이용하여 GFF와 FASTA file을 파싱하려니 여간 성가신 일이 아니다. 혹시 인터넷을 뒤지면 방법이 나오지 않을까?

역시... 친절한 shell script가 하나 나타났다. bedtools와 samtools를 기반으로 돌아가는 스크립트인 16S_extraction_Barrnap.sh라는 것이 존재하는 것이다.

https://github.com/raymondkiu/16S_extraction_Barrnap

핵심은 다음과 같이 달랑 네 줄이 전부이다. 중요한 명령어는 굵은 글씨로 강조하였다.

grep '16S' $gfffile > 16S-gff.gff;
bedtools getfasta -fi $fastafile -bed 16S-gff.gff -fo 16S-fasta.fna;
grep -m 1 ">" 16S-fasta.fna|sed 's/>//g' > 16S-id.txt;
xargs samtools faidx 16S-fasta.fna < 16S-id.txt > $fastafile-16S.fna

첫단계에서는 Barrnap 결과물인 GFF 파일에서 16S ribosomal RNA feature만 추출한다. 유전체에는 여러개의 16S rRNA가 존재하므로 두번째 줄에서 bedtools getfasta를 이용하여 서열을 추출하면 multiple fasta file(16S-fasta.fna)이 생긴다. 세번째 줄에서는 가장 먼저 출현하는 16S rRNA의 ID 줄(예: >CP013254.1:1163436-1164962)을 뽑아내어 '>'를 제거한 뒤 16S-id.txt에 저장한다. 마지막 명령어(samtools faidx)에서는 16S-id.txt에 저장된 서열 ID에 해당하는 것을 16S-fasta.fna에서 뽑아내는 것이다.

16S-id.txt 파일에는 서열 ID가 하나만 들어있으므로 사실은 마지막 행에서 xargs를 쓸 필요는 없었다.

bedtools와 samtools를 이런 용도로 사용해 본 적은 없었다. 기껏해야 feature와 관련된 수치를 뽑거나, read mapping을 하는 파이프라인 실습 자료가 지시하는 대로 기계적으로 활용했을 뿐이었다. 특히 samtools faidx는 FASTA sequence file의 인덱스 생성에만 쓰는 것으로 생각하였었다. 매뉴얼을 찾아보니 samtools faidx는 'index/extract FASTA'라고 설명이 되어있다.

Multiple fasta file로부터 특정 ID의 서열을 추출하는 용도로는 BioPerl을 사용하여 적당히 코딩한 것을 사용해 왔었는데, 이렇게 samtools를 쓸 수도 있다는 것을 오늘 처음 알았다. 단, 뽑아낼 서열이 별도의 텍스트 파일에 존재하는데 그 수가 여럿이라면(id1, id2, id3...idN) 다음과 같이 쓰거나,

samtools faidx 16S-fasta.fna id1 id2 id3 ... idN

혹은 위에서 보인 것처럼  xargs를 써야 한다. 내가 만든 스크립트는 서열 ID를 수록한 파일이 여러 라인으로 되어 있을 것을 감안한 것이었다.