Long read sequencing 결과물 다루기

외부 시퀀싱 업체에서 long read 기반의 미생물 유전체 해독을 의뢰하는 경우 대부분 de novo assembly를 하여 완성 수준의 contig 서열을 만들어 제공하게 된다. 그러나 raw data를 직접 다루거나 자체적으로 조립을 실시하여 여러 가지 측면에서 평가를 하는 것이 중요할 때가 있다. Canu 또는 SPAdes와 같이 널리 쓰이는 assembler는 short/long read data의 단독 조립 및 hybrid assembly 기능도 제공한다. 호주 멜버른 대학에서 제공하는 Long read assembly workshop의 자료를 검토해 보는 것도 본 과정을 이해하는데 도움이 될 것이다.

시퀀싱 업체에서 제공하는 zip 파일의 압축을 풀면 각 SMRT cell마다 생성되는 디렉토리 하위에 Analysis_Results라는 서브디렉토리가 보일 것이다. 그 내부로 들어가면 각각 세 개씩의 subreads.fasta/subreads.fastq/bax.h5와 하나의 bas.h5 파일이 존재한다. Bas.h5/bax.h5 파일은 PacBio RSII 장비의 primary analysis pipeline을 거쳐서 만들어지는 주된 출력물이다(bas.h5 Reference Guide). 이들 HDF5 파일은 SMRT Analysis 또는 SMRT Link 프로그램을 설치하여 PacBio data의 de novo assembly나 resequencing analysis를 실시할 목적이 아니라면 직접 이용할 일은 많지 않다.

외부에서 입수한 PacBio raw data에는 일반 연구자에게 친숙한 fasta/fastq 파일이 없이 HDF5 파일만 들어있는 경우가 있다. 이러한 상황에서는 pbh5tools에서 제공하는 bash5tools.py(사용법)을 사용하여 필터링을 거친 read를 추출할 수 있다. bash5tools.py는 시스템에 설치된 python 환경에 존재한다. --readType으로 선택할 수 있는 것은 ccs, subreads 및 unrolled이다.

$ cd Analysis_Results
$ bash5tools.py *bas.h5 --outFilePrefix myreads --outType fasta --readType subreads --minReadScore 0.75

진핵생물의 유전체의 경우 20x 정도의 시퀀싱 데이터만 있어도 canu를 사용하면 현존하는 hybrid method를 능가하는 결과를 얻을 수 있다. 그러나 최소 30~60x 정도의 데이터로 조립을 시작하는 것이 바람직하다. 다음은 25x의 대장균 PacBio 시퀀싱 데이터를 다운로드하여 조립하는 명령어이다. Nanopore data는 -nanopore 옵션을 사용한다(샘플 데이터 다운로드 링크). Input read는 FASTA/Q 무엇이든 관계가 없으며, gz/bz2/xz로 압축이 되어 있어도 좋다. 모든 결과 파일은 -d 옵션으로 지정된 디렉토리 아래에 -p 옵션으로 지정한 문자열을 접두사로 택하여 만들어진다. 입력 파라미터로서 target genome의 추정되는 크기를 제공해야 한다.

$ curl -L -o pacbio.fastq http://gembox.cbcb.umd.edu/mhap/raw/ecoli_p6_25x.filtered.fastq
$ canu -p ecoli -d ecoli-pacbio  genomeSize=4.8m -pacbio pacbio.fastq

조립 결과물 중 <prefix>.contigs.fasta는 unique 혹은 repetitive element를 모두 망라한 서열을 수록한다. 반면 <prefix>.unitigs.fasta는 alternate path에서 분리된 서열을 저장한다. 서열의 description 항목에는 다음과 같은 metadata가 기록되어서 결과를 해석하는데 도움을 준다.

>tig00000001 len=2449075 reads=6437 class=contig suggestRepeat=no suggestBubble=no suggestCircular=yes
>tig00000002 len=27634 reads=884 class=contig suggestRepeat=no suggestBubble=no suggestCircular=no

Contig와 unitig에 대한 보다 정확한 의미는 Celera Assembler Technology 위키 문서를 참조하라. ReasearchGate에 올라온 질문과 답(Contig vs. Unitig), 그리고 여기에 링크된 PacBio hybrid assembly 매뉴얼도 도움이 된다.

위에서 소개한 명령어는 raw long read의 correction, trimming, 그리고 조립을 한 번의 명령어로 진행하는 사례를 보인 것이지만, -correct, -trim 및 -assemble 옵션과 함께 세부적인 파라미터를 조정하여 단계별로 진행하는 것도 가능하다.

UniCycler는 세균의 유전체 조립을 위하여 만들어진 프로그램이다. 일루미나 데이터만 있다면 SPAdes-optimiser로 작동하는 반면, long-read data만 있는 상황에서는 minisam+Racon 파이프라인으로 작동한다. 그러나 일루미나 데이터와와 long read를 함께 사용하는 hybrid assembler로서 작동할 때 가장 정확한 결과를 생성한다.

일반적인 hybrid assembly에서는 short read(일루미나)를 long read에 정렬하여 오류를 정정한 뒤 이를 이용하여 overlap-layout-consesus 기법의 조립을 수행하게 된다. 그러나 UniCycler의 hybrid assembly에서는 일루미나 데이터를 사용하여 먼저 SPAdes로 graph 형태의 assembly를 만든 뒤 여기에 long read를 더하여 repeat을 해소하고 complete genome sequence를 얻어내게 된다. 특히 circlator와 pilon을 마지막 단계에 실행하므로 별도의 post-process를 거치지 않아도 된다.

사용하는 thread의 수는 -t 또는 --threads 옵션으로 지정하지 않으면 8을 기본으로 택한다. 사용하는 thread의 수는 -t 또는 –threads 옵션으로 지정하지 않으면 8을 기본으로 택한다. 실행 모드는 conservative, normal 및 bold의 세 가지가 있는데, –mode 옵션을 통하여 지정한다. 특별히 설정하지 않으면 normal mode로 작동한다. |실행 모드는 conservative, normal 및 bold의 세 가지가 있는데, --mode 옵션을 통하여 지정한다. 특별히 설정하지 않으면 normal mode로 작동한다.

$ unicycler -1 short_1.fastq -2 short_2.fastq -l long.fasta -o OUT_DIR -t 16

신뢰할 수 있는 GFA 포맷의 long read assembly가 있다면 --existing_long_read_assembly 옵션으로 지정하면 된다.

2019년 Nature Biotechnology를 통하여 공개된 Flye는 원핵생물 유전체 조립 성능 비교를 위한 벤치마킹에서 Miniasm/Minipolish와 더불어 가장 우수한 성능을 나타낸 것으로 알려져 있다(참고 문헌).

-t 또는 --threads 파라미터는 지정하지 않을 경우 1이 주어지므로 컴퓨터의 상황에 맞게 적절히 설정한다. --genome size는 유전체 크기의 0.5에서 2배까지 여유롭게 설정하여도 결과에는 큰 영향을 미치지 않으니 정확하게 값을 주지 않아도 된다. 거대 유전체를 조립하는 경우 메모리 사용량을 줄이기 위해 긴 read 위주로 초기에 사용할 데이터를 제한할 수 있다(--asm-coverage). 보통 40x coverage면 무난하며, 조립 마지막 단계에서는 모든 데이터를 다 사용하게 되므로 걱정할 필요는 없다. --plasmids 옵션은 short unassembled read를 extra contig로 처리하는 별도의 단계를 거치게 한다.

$ flye --threads 16 --nano-raw barcode01.fastq --out-dir results --asm-coverage 40 --genome-size 4m --plasmids

Flye는 assembly_info.txt라는 tab-delimited text 파일에 contig의 길이와 coverage 및 원형 구조 등에 대한 정보를 기록한다. Circ. = Y로 표시된 contig는 이미 circularization이 되었기 때문에 양 말단에 중복되는 염기서열이 존재하지 않는다.

$ cat assembly_info.txt 
#seq_name	length	cov.	circ.	repeat	mult.	alt_group	graph_path
contig_2	3646351	115	Y	N	1	*	2
contig_1	5973	486	Y	Y	4	*	1

LastGraph(Velvet), FASTA(SPAdes), 또는 GFA 등 de novo assembler가 만들어 내는 graph의 시각화에는 Bandage를 활용한다. 이 프로그램은 Mac/Linux/Windows용을 같이 제공한다.

$ Bandage load assembly_graph.gfa

인수 없이 Bandage라고만 입력하면 Bandage GUI가 작동하게 되며, File → Load graph에서 그래프를 로드한 뒤 ‘Drawu graph’ 버튼을 클릭하면 화면에 시각화된 그래프가 표현된다.

일루미나 시퀀싱 데이터에 비하여 indel 형태의 오류를 상당한 수준으로 내포하고 있는 long read는 일반적인 read mapping 소프트웨어를 이용하여 reference 서열에 매핑하기가 어렵다. PacBio에서는 자사의 long read에 특화된 long read aligner인 BLASR을 제공하고 있다. BWA가 제공하는 세 가지 알고리즘인 BWA-backtrack, BWA-SW 및 BAW-MEM 중에서 BWA-MEM은 PacBio와 Nanopore의 sequencing read의 mapping 용으로 쓰여 왔는데, 최근에는 성능이 더 좋은 minimap2로 대체되었다.

# Map short accurate genomic reads
$ minimap2 -ax sr ref.fa read1.fq read2.fq > aln.sam
# Map long and noisy genomic reads
$ minimap2 -ax map-pb ref.fa pacbio-reads.fq > aln.sam # for Pacbio subreads
$ minimap2 -ax map-ont ref.fa ont-reads.fq > aln.sam # for Oxford Nanopore reads
# Convert SAM to BAM
$ samtools view –b –S –o aln.bam aln.sam # samtools version # 1.9 (for '-o' option)