10년도 더 된 옛날 만든 미생물의 일루미나 시퀀싱 데이터를 정리하여 국가바이오데이터스테이션(K-BDS)에 등록하는 '유전체 박물관(고물상?)'에 열심히 매진하는 과정에서 데이터 검수 책임을 맡는 관리자에게 늘 고마운 마음을 갖게 된다. 등록 당시에 포함된 어이없는 오류를 바로잡기 위해 번거로운 부탁을 해도 매번 잘 수용해 주기 때문이다. 아직 GenBank의 서열 정보와 같이 버전 관리가 이루어지고 있지는 않지만, 앞으로는 이 기능이 포함되리라 믿는다.
최근 품질관리 담당자로부터 내 HiSeq 2000 raw sequencing data(FASTQ 파일)를 검수하는 과정에서 동일한 서열 ID가 중복하여 나타난다는 소식을 들었다. 당시에는 한 샘플의 FASTQ 파일이 너무 크면 001, 002...와 같은 이름을 붙여서 분할한 형태로 제공하였었다. 이를 'cat' 명령어로 하나의 FASTQ 파일로 합치는 과정에서 똑같은 원본 파일을 두 차례 반복하여 이어 붙인 것이 원인이었다.
72개의 샘플 중에서 이런 문제가 발생한 것은 6가지였다. 문제를 수정한 뒤 모든 분석 작업, 즉 jellyfish를 이용한 21-mer abundance analysis, phyloFlash, ZGA pipeline 및 GTDB-Tk를 전부 새로 실시하였다. 특히 GTDB-Tk는 공식 문서를 기계적으로 따라하다가 너무 오래된 소프트웨어 및 reference DB를 설치하였음을 발견하고 최근 재설치를 하였기에 어차피 모든 데이터에 대한 재작업이 필요한 상태였다. Sequence read ID가 전부 2회씩 중복하여 존재함에도 불구하고 내가 사용한 그 어떤 소프트웨어도 오류 메시지를 뱉어내지 않았었다. 아마도 프로그램 내부에서는 read ID를 전부 단순한 일련번호와 같은 것으로 치환하여 쓰고 있는 것이 아닐까 한다.
수정된 FASTQ 파일을 이용하여 재분석을 하던 과정 중에 Methanobrevibacter smithii의 어느 균주에 대한 시퀀싱 데이터가 이상하다는 것을 발견하였다. FASTQ 파일에 오류가 있었을 때, 즉 정상에 대하여 2배로 부풀려진 상태에서는 contig의 수가 500개에 육박하였고 CheckM(실제로는 ZGA pipeline 내부에서 실행)에서 계산한 completeness가 40%도 되지 않았었다. 그래서 sequencing coverage가 상당히 부족해서 문제가 생겼을 것이라 생각하고 더 이상 쓸 수 없는 데이터로 간주하였었다(KAP241424). 21-mer abundance analysis에서도 main peak가 그렇게 잘 드러나지 않았기에 이런 생각을 하는 것도 당연하였다.
그런데 total contig length를 보니 약 1.8 Mb였다. 이는 Methanobrevibacter smithii의 유전체 크기로 알려진 값과 매우 유사하였다. Contig의 수가 지나치게 많이 나올 정도로 sequencing coverage가 충분하지 않다면, total contig length도 매우 낮은 수준이어야 한다. 표준 균주인 ATCC 35061의 유전체를 NCBI에서 다운로드한 뒤 ZGA로 CheckM 점검을 해 보았더니 이것 역시 completeness가 40% 미만으로 나타났다. 이건 말이 되지 않는다. 두 유전체 서열에 대하여 CheckM을 직접 실행해 보았다. ZGA가 설치된 conda 환경을 그대로 이용하였으니 CheckM 버전은 ZGA pipeline의 것과 동일하다.
결과는 놀랍게도 두 균주 모두 completeness가 100%인 것으로 계산되었다. ZGA pipeline 안에서 CheckM을 실행하고 그 메시지를 처리하는 과정에 문제가 있는 것이 틀림없다(나중에 확인한 것이지만 이는 나의 오해였음). 따라서 내가 시퀀싱한 Methanobrevibacter smithii 균주는 아무런 문제가 없으며, 분리 정보만 되찾아 정리한다면 K-BDS에 등록하여도 아무런 품잘 상의 하자가 없으리라는 결론을 내렸다.
수정한 FASTQ 파일을 이용하여 조립을 하면 contig의 수는 겨우 41개에 불과하다(~156x coverage). 이 데이터를 그대로 duplicate하여 분량을 두 배로 만든 상태에서 조립을 하였더니 contig 수가 500개 가까이 된다? 사실 이 현상은 나의 지식이나 경험으로는 잘 이해가 가지 않는다. 일루미나의 경우 대략 150x 정도의 sequencing coverage에서 최적의 결과가 나오게 되고, 투입량을 그 이상으로 올린다고 해도 더 좋아지지는 않는다는 논문을 본 일은 있다. 하지만 투입량을 그 이상 올렸을 때 조립 결과가 심각하게 나빠진다는 것은 이해하기 어렵다.
이번 데이터에 포함된 오류는 단순히 시퀀싱 분량이 2배로 늘어난 것과는 상황이 다르다는 생각이 들었다. 시퀀싱을 한 판 더 실시하여 read를 2배로 얻었다 해도 이전의 것과 모든 염기서열과 quality가 똑같은 read가 발생할 가능성은 극히 적다. 하지만 내가 초래한 실수에서는 그렇지 않다. 완벽하게 똑같은 read가 2개씩 존재하고 있다. 그렇다면 k-mer spectrum에서 가장 왼쪽에 있는 것들(시퀀싱 오류에 의해서 1,2...회 정도의 저빈도로 존재하지만 그 종류는 가장 많은 것)이 둘 씩 만나서 짧은 contig를 다수 만들게 되는 것이 아닐까?
아직 이전 FASTQ 파일로 만든 contig의 길이 분포 및 평균 read depth를 점검해 보지는 않았다. 버릴 데이터라고 생각해서 신경을 쓰지 않았기 때문이다. 만약 이 assembly에서 짧은 low depth contig가 많은 수를 차지한다면, 나의 원인 분석이 그렇게 틀리지는 않았다는 결론을 내릴 수 있을 것이다.
그래서 '망친' 데이터라고 해서 함부로 내다 버려서는 안된다고 생각한다. 조작 실수에 의해 '망친' 것으로 오해를 할 수도 있기 때문이다. 망친 데이터의 모임인 KAP241424는 전반적인 검토를 거친 뒤 데이터와 README 파일을 바꾸어야 한다. 물론 오늘 글의 주제가 된 Methanobrevibacter smithii의 데이터는 한번 더 확인을 거쳐서 '부활'하게 될 것으로 본다.
2025년 3월 5일 업데이트 - CheckM 실행 결과의 고찰
ZGA pipeline 내부에서 CheckM을 실행할 때에는 'checkm taxonomy_wf' subcommand를 쓰되 기본적으로 bacterial marker를 사용한다. 만약 고세균이라면 '--domain archaea'라는 옵션을 주어야 한다. 문제는 바로 여기에 있었다. 여러 genome에 대하여 일괄적으로 ZGA pipeline을 돌렸으니 개별적으로 domain 정보를 주지 않았고, default setting인 bacteria 마커셋을 사용했던 것이다. 그러니 당연히 genome completeness가 40%도 되지 않는 값이 나올 수밖에... Standalone CheckM 실행에서 나는 'checkm lineage_wf' subcommand를 사용했다. 이 경우에는 domain을 소프트웨어가 알아서 판단하여 적용하는 것으로 보인다. 따라서 ZGA pipeline은 아무런 잘못이 없으며, 매뉴얼을 꼼꼼하게 읽지 않은 나의 잘못이었다.
