염기서열 정보가 포함된 GFF에서 이를 떼어내려면 어떻게 하면 좋을까? 파일을 라인 단위로 훑다가 '##FASTA'라는 라인을 만나게 되면 여기서부터 파일의 끝까지는 무시하게 만들려면 어떻게 할까?
여기서 잠깐! GenBank 파일이 복수의 서열을 포함하고 있으면 어떻게 하나? Draft assembly이거나, 혹은 플라스미드를 포함하는 박테리아의 RefSeq assembly(GCF_어쩌고저쩌고..)가 이러한 경우에 해당한다. 전혀 걱정할 것이 없다. GenBank 파일에서는 '//'을 구분자로 분리된 각 레코드가 서열 정보를 개별적으로 갖고 있지만, 이를 bp_genbank2gff3.pl로 처리하면 서열이 몇 개나 존재하든지 관계없이 모든 염기서열은 GFF 파일의 뒷부분에 모이기 때문이다.
그러면 다시 본론으로 돌아가서 '##FASTA'로 시작하는 라인부터 끝부분까지를 제거하는 방법을 알아보자. 웹 검색을 전혀 하지 않고 머릿속에 들어있는 지식만 가지고 해결하려면 Perl로 몇 줄 쓰면 된다. 하지만 쉘에서 돌아가는 일반 유틸리티를 이용하여 한 줄의 명령으로 가능할 것 같다.
이런 해결 방안이 있다. bp_genbank2gff3.pl에서 -r[--noinfer] 옵션은 박테리아 유전체 정보에서는 별로 중요하지 않은 exon/mRNA subfeature를 추론하지 말라는 뜻이다. 결과를 표준 출력으로 내보낸 뒤(-o -), sed stream editor를 돌리면 된다. 그런데 뭐가 이렇게 간단한가?
$ bp_genbank2gff3.pl -r input.gbk -o - | sed -n '/^##FASTA/q;p' > output.gff
원리를 설명하고 싶은데, 그게 쉽지가 않다. 원리를 완벽하게 이해하지 않고 쓰는 것은 약간 위험하지만, 내가 Perl을 처음 공부할 때도 그렇지 않았던가? 비록 지저분한 코드라 하더라도 일단 작동하는 코드를 만들어라. 주석을 달고 개선하는 것은 그 다음이다. 일단 -n 옵션의 의미부터 알아보자. sed의 기본 동작은 입력 파일을 줄 단위로 일단 출력해 놓은 다음, 명령에 따라 행동을 하는 것이다. -n은 명시적으로 인쇄를 하라는 명령을 제공해야한 인쇄를 하라는 뜻이다. sed 'p' file을 입력해 보라. 모든 줄이 두번씩 인쇄될 것이다. sed를 cat처럼 이용하려면 sed -n 'p' file이라고 실행해야 한다. 마찬가지로 sed를 grep처럼 사용하려면 sed -n '/PATTERN/ p' file이라고 입력해야 한다.
만약 어떤 패턴이 매치된 라인의 다음 줄부터 마지막까지를 인쇄하고 싶다면 이것도 sed로 해결할 수 있다. Eric's sed one-liner의 제4부 Selective Printing of Certain Lines에 그 해답이 있다. 나머지 파트도 숙독해 보면 놀라운 sed의 능력을 알게 될 것이다. sed와 awk의 기능을 살펴보다가 늘 '왜 내가 이런 기능을 Perl로 짜려고 했던 것이지?'하는 결론을 내리고는 한다. Eric Pement의 웹사이트에 마련된 sed, awk 및 Perl 정보가 많은 영감을 줄 것 같다. 그가 정리한 sed one-liner 사례 모음인 sed1line.txt를 보고 있노라면 경이롭기까지 하다.
오늘의 글을 끝맺는 한 마디.
Don't reinvent the wheel.
댓글 없음:
댓글 쓰기