사실 컴퓨터의 코어 수가 충분하다면 basecalling에 걸리는 시간은 충분히 줄이는 것이 가능하다. MinION을 구동하는데 사용한 낡은 Xeon E5520(2.27 GHz)에서 basecalling을 그대로 진행하느라 이렇게 된 것이다. 이 CPU는 코어 4 개를 갖고 있어서 동시 실행 스레드는 8 개에 불과하다. 따라서 local basecalling을 해서는 소용이 없고 /var/lib/MinKNOW/data/reads/tmp 아래에 저장된 백만 개 이상의 raw fast5 파일을 코어가 많은 다른 서버로 복사한 뒤 albacore(실제로는 read_fast5_basecaller.py 스크립트)를 수행해야 한다. 복잡한 디렉토리 구조 하에 파일이 존재하므로 rsync를 사용하는 것이 유리하다. 다행히 read_fast5_basecaller.py는 recursive하게 실행되므로 fast5 파일을 한 곳에 옮기지 않아도 된다.
그러나 파일의 확장자(.fast5.tmp)는 .fast5로 바꾸어야 한다. 파일 이름의 맨 끝에 있는 .tmp를 어떻게 recursive하게 한꺼번에 제거할 것인가?
$ find tmp -name "*.fast5.tmp" -type f | while read filedomv ${file} ${file%.tmp}done세번째 줄은 ${file%%.tmp}라 해도 이번 경우에서는 똑같은 결과가 된다. Bash string manipulation examples - length, substring, find and replace(링크)에서 풍부한 사례를 맛볼 수 있다.
24 시간의 러닝 타임 동안 3만 개가 조금 못되는 fast5 파일(이 경우에는 basecall이 완료된)과 8개의 fastq 파일이 생겼다. 이것만으로도 충분한 분량이므로, PacBio 데이터를 다루면서 가장 익숙하게 사용해온 de novo assembler인 canu를 사용하여 조립을 해 보았다. Contig는 총 2 개가 생성되었으며 각각의 크기는 2G, G + r이다. G는 예상되는 genome size이고 r은 redundant한 서열이다. 이는 circle 형태의 DNA를 대상으로 시퀀싱 라이브러리를 만들어서 염기서열을 만든 뒤 조립하면 흔히 발생하는 현상이다. Circlator를 실행하니 작은 contig는 큰 것에 병합이 되었다. 하지만 redundant한 서열은 없어지지 않았다. 남은 contig(size = 2G)를 가지고 gepard로 자체 dot plot을 그려 보았다.
이것은 무슨 의미인가? 서열 단위가 tandem하게 반복됨을 뜻한다(=====>======>). 원래 circlator는 이를 알아서 처리하여 하나의 반복 단위(즉 1G = genome size에 해당하도록)로 만들어 주어야 한다. 하지만 기대와 달리 작은 contig를 큰 contig에 병합하는 것 말고는 원형화('circularize')를 시키지 못했음을 의미한다. cross_match 결과는 다음과 같다.
약 49 kb의 단위가 두 차례 반복되지만 서열이 100% 동일하지는 않다. Phage genome이 원래 이렇게 생겼을 가능성은 극히 적고, 아마도 long read sequencing 고유의 error 문제일 것으로 생각이 든다.
그렇다면 다른 assembler를 사용하면 한번에 깔끔한 결과가 얻어질까? ONT의 공식 소프트웨어라고 할 수 있는 Pomoxis를 써 보기로 하였다. 나는 GitHub의 클론을 만들어 설치하지 않고 각 구성 프로그램(minimap2, miniasm 및 racon)을 별도로 설치하였다. 실제 실행 방법은 University of Washington의 Genome Assembly - minimap/miniasm/racon overview 페이지(링크)를 참조하였는데, 여기에 약간의 오류가 있으니 주의하기 바란다.
그러면 각 단계를 살펴보도록 하자.
1. minimap2로 long read를 자체 비교하기
minimap2 -t 16 -x ava-ont ont-reads.fastq ont-reads.fastq > ovlp.paf
Self-overlap을 검출하기 위해서 하나의 read file이 연속해서 인수로 주어졌다. PAF 파일은 base-level alignment가 없는 approximate mapping 결과물이다. Reference sequence에 read를 매핑하거나 spliced alignment를 할 경우, 그리고 결과물을 SAM으로 얻으려면 각각 이에 맞는 옵션을 주어야 한다. Read file이 여러 개라면 하나로 합쳐야 한다. PacBio read를 사용하는 경우 -x ava-pb 옵션을 적용한다.
2. miniasm으로 조립하기
miniasm -f ont-reads.fastq ovlp.paf > miniasm_reads.gfa
3. FASTA file 추출하기
awk '$1 ~/S/ {print ">"$2"\n"$3}' miniasm_reads.gfa > miniasm_reads.fasta
4. Contig 서열에 다시 read를 매핑하기(minimap2)
이는 racon 단계에서 필요한 overlap 정보 파일(HMAP/PAF/SAM)을 얻기 위함이다.
minimap2 -t 16 miniasm_reads.fasta ont-reads.fastq > ovlp2.paf
5. racon으로 consensus 추출
racon -t 16 ont-reads.fastq ovlp2.paf miniasm_reads.fasta > consensus.fasta
miniasm 결과는 길이 48328 bp의 contig 하나였고, racon 교정 후에는 49546 bp가 되었다. 혹시 이 contig 서열의 양 끝에서는 더 이상 정돈할 것이 없을까? 이것이 남은 숙제이다.
댓글 없음:
댓글 쓰기