User Tools

Site Tools


Sidebar

This is the sidebar. Without it, the main text is too wide!

ont_sequencing_data_analysis

This is an old revision of the document!


Data analysis for Oxford Nanopore sequencing

Oxford Nanopore Technologies(ONT)에서 제공하는 모든 공식 문서는 http://community.nanoporetech.com/을 참조하라(로그인 필요). PoreCamp에서도 매우 유용한 자료를 제공한다. 2016년에 영국에서 열렸던 PoreCamp의 교육자료에서는 데이터 분석 방법과 필수 소프트웨어에 대해서 많은 영감을 얻을 수 있다.

  • MinION Mk1B (product code MIN-101B)
  • Flow cell: Spot-On Flow cell MkI (product code FLO-MIN106 R9)
  • 구동용 컴퓨터 I: Xeon E5520 @2.27, 16 GB memory, Ubuntu 14.04.5 LTS. (uname -a 명령 시 4.4.0-134-generic #160~14.04.1-Ubuntu). USB 3.0 인터페이스 카드는 Via 805를 사용한 ipTIME PX300 Plus) MinKNOW v2.0.0으로 잘 동작.
  • 구동용 컴퓨터 II: microbe server(Xeon E5-2640 x 2, SuperMicro X9DRI-F)에 같은 OS와 카드를 설치하였으나 인식이 되지 않음. /etc/default/grub 파일에서 GRUB_CMDLINE_LINUX를 다음과 같이 변경하고 sudo update-grub을 한 뒤 재부팅을 하면 된다고 하였으나(참고 1 참고 2) 여전히 여기에 꽂힌 USB 장비를 읽지 못함. 그러나 USB 3.0 카드를 다른 위치에 꽂으니 정상 작동함. 세상에 이럴 수가…
GRUB_CMDLINE_LINUX="iommu=soft"
GRUB_CMDLINE_LINUX_DEFAULT="quiet splash pciehp.pciehp_force=1"

EPI2ME에서는 더 이상 basecall을 진행하지 않으므로 Albacore를 설치하여 local basecalling을 하라는 고객 센터의 알림이 있었다. 단, MinKNOW 설정 창의 Basecalling(Live or None 중에 선택)에서 나타나는 Live basecalling은 Albacore를 뜻하는 것이 아니라 1D sequencing protocol을 위해 내장된 것이다. 파이썬 2.7과 3.5(Albacore)를 번갈아 이용해야 하므로 pyenv를 활용하는 것을 권장한다. 혹은 anaconda도 좋다.

새로 시작하기

이 글은 2019년 5월 2일부터 작성하기 시작한다. 다음 섹션(처음에 해야 할 것) 이후의 글은 더 이상 유효하지 않다. MinKNOW를 이용하면 albacore를 별도로 설치하지 않아도 basecalling이 진행되고, 새로운 pomoxis에서는 read 간의 overlap을 de novo assembly 전에 산출할 필요가 없다.

Best long read mapper?

한동안 BWA-mem이 널리 쓰여 왔으나, 이제는 pomoxis에 포함되어 있는 minimap2가 최선인 것으로 보인다(참고 글

처음에 해야 할 것

프로그램 설치

https://community.nanoporetech.com/downloads (Nanoporetech login required)

  • MinKNOW
  • GridION software
  • EPI2ME Desktop Agent
  • EPI2ME Command Line Agent
  • Albacore (local basecalling for MinKNOW)

Configuration test

Platform QC

Flow cell을 수령하면 5일 이내에 platform QC를 실시하여 포어가 어떤 상태인지를 확인해야 한다. 보관된 flow cell에 라이브러리 로딩 직전에도 하는 것이 좋다.

MUX란 무엇인가?

ONT 사이트의 Interpreting the sequencing run 동영상 강좌에서 platform QC 및 mux scan의 설명, 그리고 sequencing run의 해석 방법에 대한 설명을 들어보라.

Official genome assembly, consensus, and polishing tools

2017년 7월 30일 Nanopore Community를 통한 공지(링크)

fast5 data file의 이해

하나의 포어에 DNA 분자가 통과하면서 발생하는 정보는 하나의 fast5 파일에 수록된다. 이는 HDF5라고도 불리는 hierarchical data format의 한 변형으로 이해하면 된다. 가장 낮은 레벨의 정보는 이벤트이다. 이벤트는 DNA 분자에 대한 aggregated signal measurement로서, MinKNOW가 기록하는 것은 이러한 pre-basecalled fast5 file이다.

만약 ONT cloud에서 basecall을 했다면(즉 Metrichor를 이용했다면) 업로드된 fast5 파일은 처리를 거쳐서 pass와 fail 디렉토리로 나뉘어 다운로드된다. pass에는 바코드에 의해 성공적으로 분리되고 사전에 정의된 threshold를 넘는 mean base quality score를 넘는 read들이 쌓인다.

QC using poretools

portools를 사용하려면 basecall 정보가 담긴 fast5이 필요하므로, MinKNOW가 생성한 fast5 파일을 그대로 사용하면 안된다. 아래에서 설명한 바와 같이 Albacore를 사용하여 basecall을 한 fast5 파일을 먼저 만들도록 한다. 처음에는 이와 같은 사실을 몰라서 왜 첫 러닝한 fast5에 아무런 정보가 없는 것일까하고 무척 고민을 많이 하였다.

poretools에서는 유용한 diagnostic plot을 만들어내는 기능이 있다. 그러나 십중팔구 다음과 같은 에러 메시지가 나올 가능성이 크다. 참고로 Tkinter는 파이썬 모듈, _tkinter는 C 모듈이다.

    import _tkinter # If this fails your Python may not be configured for Tk
ImportError: No module named _tkinter

이 문제를 해결하는 것은 상당히 까다롭다. 우분투는 python-tk와 tk-dev, CentOS는 tkinter와 tk-devel을 설치하면 된다고는 하지만 시스템의 기본 파이썬 버전이 실제 사용할 것과 다른 경우(pyenv 등을 사용하는 환경) 이를 매끄럽게 설치하기가 무척 어렵다. 며칠동안 고생하여 얻은 결론은, yum으로 설치한 tkinter는 시스템 고유의 python과 같은 버전일 수밖에 없다는 것이다. pyenv로 제아무리 파이썬 2.7.x를 설치하였다고 해도, tkinter는 여전히 2.6.x이다. 그러면 tkinter-2.7.5-48.el7.x86_64.rpm를 설치하면 안되나? 이것은 CentOS 7용이고, python 2.7.5-48이 필요하다. CentOS 6에서 파이썬을 2.7 버전으로 올리면 yum이 문제가 된다고 하였다.

Successful poretools installation under pristine-installed CentOS 7 (1611)

System python (version 2.7.5) was used.

# yum install tkinter tcl-devel tk-devel hdf5-devel
# git clone https://github.com/arq5x/poretools
# cd poretools
# python setup.py install

Basecalling (local - using Albacore)

MinKNOW가 생성한 fast5 파일은 fast5/0/에 존재하는 상황이라서 -r(–recursive) 인수가 필요하다. basecall 결과는 fast5와 fastq 전부 씌여지게 하였다. fastq 파일에는 여러 read가 한꺼번에 수록된다(-q or –reads_per_fastq_batch로 설정 가능; 기본값은 4000).

$ read_fast5_basecaller.py -l # or --list_workflows
$ read_fast5_basecaller.py -i fast5 -r -t 2 -s fast5_albacore/ -f FLO-MIN106 -k SQK-LSK208 -o fast5,fastq

인수 설명

long form을 보면 무엇을 의미하는지 명확히 알 수 있다. 필수 옵션(값을 동반해야 함)은 굵은 글씨로 표현하였다.

  • -i or –input: input path. fast5 파일이 있는 곳.
  • -t –worker_threads
  • -s or –save_path: basecalled .fast5 파일과 다른 정보가 저장되는 곳. 미리 만들 필요는 없다.
  • -f –flowcell
  • -k or –kit
  • -c or –config
  • -o or –output_format: fast5와 fastq를 같이 만들려면 -o fast5,fastq처럼 콤마만으로 연결하면 된다. 이 인수가 없으면 fastq 파일만 만들어진다.
  • -r or –recursive: –input_path 하위의 서브디렉토리에 fast5 파일이 존재할 경우 사용.
  • –barcoding
  • -n or –files_per_batch_folder
  • -q or –reads_per_fastq_batch

임시 디렉토리의 파일 확장자 일괄적으로 바꾸기

tmp 디렉토리 아래에 저장되는 raw fast5 파일의 이름은 .fast5.tmp로 끝나기 때문에 read_fast5_basecaller.py 스크립트에 그대로 공급할 수 없다. 다음 방법을 이용하여 recursive하게 일괄 변경하면 된다. fast5 파일이 백만 개를 넘으면 이것도 시간이 꽤 많이 걸린다.

$ find tmp -name "*.fast5.tmp" -type f | while read file
do
mv ${file} ${file%.tmp}
done

Mapping

bwa 0.7.12-r1039, samtools 1.4.1 조건에서 다음을 실시하였다. BWA-MEM 'ont2d' 옵션은 버전 0.7.11부터 포함되었다. samtools view의 옵션에서 -S는 입력이 SAM임을 명시하기 위함이다. 예전 버전의 samtools에서는 이를 꼭 지정해 주어야 하지만 요즘 버전은 자동으로 인식하므로 필요가 없다. 그러나 하위 호환성을 위해서 그대로 둔 것이다.

$ bwa index reference.fasta
$ bwa mem -x ont2d reference.fasta sample2D.fasta | samtools view -bS - | samtools sort -o sample2D.sorted.bam
$ samtools index sample2D.sorted.bam
$ samtools stats sample2D.sorted.bam > sample2D.stats.txt
$ head -n 40 samples2D.stats.txt
$ grep '^COV' samples2D.stats.txt > samples2D.coverage.txt
R> library(ggplot2)  
R> cov=read.table("/path/to/your/sample2D.coverage.txt", sep="\t")  
R> cov[1,]  
R> ggplot(cov, aes(x=V3, y=V4)) + geom_bar(stat='identity') + xlab('coverage') + ylab('count')

BAM file로부터 다음의 질문에 대답을 할 수 있어야 한다.

  • How many reads were mapped?
  • What was the average length of the reads?

Assembly의 교정

Racon

Racon is intended as a standalone consensus module to correct raw contigs generated by rapid assembly methods which do not include a consensus step. Racon의 실행에는 MHAP/PAF/SAM format의 overlap 정보가 필요하다. 이것은 반드시 염기서열 수준의 alignment가 필요함을 의미하는 것은 아니다.

$ racon [options ...] <sequences> <overlaps> <target sequences>

Mapping file이 필요하다. Canu로 만든 assembly를 교정하는 사례를 생각해 보자. overlap 정보를 어떻게 만들 것인가?

  1. bwa를 사용? 그렇다면 raw read 혹은 corrected (trimmed) reads?
  2. minimap2를 사용? 그렇다면 raw read 혹은 corrected (trimmed) reads?
$ bwa index CANU.unitigs.fasta
$ bwa mem -x ont2 -t 16 CANU.unitigs.fasta ... (작성 예정)

Medaka

https://github.com/nanoporetech/medaka https://nanoporetech.github.io/medaka/

python 3.5 환경이 필요하므로, conda py35 environment로 전환하여 설치하였다.

$ cd /data/apps/medaka/ # 설치한 뒤
$ source venv/bin/activate

NanoOK

NanoOK(나눜이라고 읽는다)이란 ONT sequencing read의 전처리, 매핑, post-analysis 및 종합적인 QC를 실시하는 도구이다. 설치에 의존성이 필요하므로 Docker 이미지를 쓰는 것이 편리하다.

데이터 파일의 준비(Albacore 기준)

workspace 디렉토리는 basecall 디렉토리를 경유하지 않고 MinION_run_directory_1/ 바로 아래에 위치시켜도 무방하다.

Main_directory/
|
+- MinION_run_directory_1/ (SampleDir)
|  |
|  +- basecall/ (Albacore 결과)
|     |
|     +- workspace/
|         |
|         +- 0/ (basecalled fast5 파일은 여기에)
|
+- references/ (reference fasta 파일은 이 하위에)
|
+- MinION_run_directory_2/

Docker 이미지 파일을 이용한 실행 방법을 기준으로 설명한다.

# docker run -i -t -v /path/to/Main_directory:/usr/nanopore richardmleggett/nanook bash
# cd /usr/nanopore/MinION_run_directory_1/references
# lastdb -Q 0 ref ref.fasta
# cd ..
# nanook extract -s MinION_run_directory_1 –f basecall/workspace
  <= fast5 read의 위치를 -f 인수로 직접 지정하려면 SampleDir(MinION_run_directory_1)을 기준으로 한 상대 경로를 주어야 함.
# nanook align -s MinION_run_directory_1 –r references/ref.fasta
# nanook analyse -s MinION_run_directory_1 –r references/ref.fasta -passonly
  <= MinION_run_directory_1/latex_last_passonly/ 아래의 pdf file을 확인하라.
# exit

Troubleshooting

  • MinION이 너무 뜨겁거나 차가울 때 해결 방법(링크)
ont_sequencing_data_analysis.1556771187.txt.gz · Last modified: 2019/05/02 13:26 by hyjeong