상용 및 오픈소스를 이용한 NGS 분석 서비스의 Mapping Rate 비교

아래 글에서 언급했던 논문의 서플을 보면 다음과 같이 NGS 분석 서비스를 비교한 테이블이 나온다. 테이블의 첫번째 나오는 것이 논문 저자들이 만든 HugeSeq이라는 파이프라인에 대한 것으로 다음과 같은 항목에 대해서 각각의 서비스들을 비교했다. 크게 항목을 나누자면 웹기반의 서비스와 직접 다운로드해서 설치하는 것으로 나눌 수 있다. 우선 웹 기반의 서비스중 대표격인 DNAnexus는 Align, SNP Call, Indel Call을 지원하지만 아직 SV (Structural Variant) Call은 지원하지 않고 있으며 commercial한 서비스이다. 즉 커머셜하면서 웹 기반으로 제공되는 서비스 중 최고를 달리고 있다.

테이블 1. NGS 데이터 분석 서비스 비교

	Alignment	SNP Calling	Indel Calling	SV Calling	Availablity	Data Size Limit	License	Functional Annotation
HugeSeq	yes	yes	yes	yes	downloadable	no	public open-source	yes
SOAP	yes	yes	yes	yes	downloadbale	no	public open-source	no
GATK	partial	yes	yes	no	downloadbale	no	public open-source	yes
Galaxy	yes	yes	yes	no	web-based downloadable	yes in online version	public open-source	yes
DNAnexus	yes	yes	yes	no	web-based	no	commercial	yes
Spiral Genetics	yes	yes	no	no	web-based	no	commerical	yes

자 그렇다면 HugeSeq과 DNAnexus, Spiral Genetics의 서비스를 각각 이용해서 동일 샘플을 분석한 결과를 한번 살펴보도록 하자. 샘플 데이터는 Illumina HiSeq 2000으로 시퀀싱된 데이터를 사용했으며, 총 read의 갯수는 1,389,159,180 (약 13억개)로 90 bp의 read length / paired-end 로 된 데이터이다.

BWA(HugeSeq)의 경우 Reference에 매핑된  read는 1,297,227,533개로 매핑률은 93.4% 로 이중 paired-end의 쌍을 만족하는 confidently하게 매핑된 것은 92%를 보이고 있다. 이처럼 일반적인 파이프라인(BWA를 사용하는 HugeSeq과 같은 공개툴을 이용한 경우)의 별다른 옵션 없이 수행했을 경우, 매핑률이 가장 높았으며, DNAnexus, Spiral Genetics의 순이다. 매핑률과 발굴된 SNP의 갯수는 별 상관 관계가 없이 DNAnexus가 가장 높고 다음으로 BWA, Spiral Genetics가 그 뒤를 잇고 있다.

테이블 2. 한국인 샘플의 reference 매핑률 (hg19 기준) - Total reads : 1,389,159,180 (694,579,590 x 2)

	Mapping Rate	Properly Mapping	Mapped Reads	SNPs
BWA	93.40%	1,277,997,460	1,297,227,533	3,552,087
DNAnexus	89.77%	1,167,367,656	1,247,083,485	4,073,143
Spiral Genetics	80.02%	1,031,585,051	1,111,657,780	2,923,278

그건 그렇다고 치고... 그렇다면 평균 10% 이상의 read들은 quality도 좋은데 왜 reference genome (hg19)에 매핑되지 않았을까? 돈 많이 들여서 시퀀싱 했는데, 10%는 한번 써보지도 못하고 버리다니,  요즘 다들 한국인 표준레퍼런스/참조레퍼런스를 만든다고들 떠들고 있는데 이걸 이용해서 align을 수행한다면 어떤 결과가 나올까?  1) 한국인 레퍼런스를 사용한다고 매핑률이 올라갈까? 2) 찾은 SNP에는 어떤 차이가 있을까?

현재 한국인에 대해서는 두편의 논문이 발표된 상황으로 각각 한국생명공학연구원의 KOBIC과 서울대학교 GMI에서 나온 결과로, consensus sequence를 fasta 포맷의 데이터로 제공하는 것은 현재 SJK(KOBIC)데이터뿐이므로 이 데이터를 이용하여 alignment를 수행한 후 한국인 SNV 데이터를 몽땅 모아서 발굴한 SNV를 filter out 시켜본다. 물론 이때 dbSNP와 1000 genomes에서 찾은  SNV도 모두 filter out 한다. 왜? 그냥 한번 해보는거다 -.-;;

SJK
ftp://bioftp.org/BiO/Store/Genome/KOREF_KoreanReferenceGenome/fasta/

다운로드된 파일은 각 염색체별로 나누어 압축되어 있기 때문에 하나의 FASTA 포맷의 파일로 만든다. 이때 간단하게 gunzip의 -c 옵션을 통해서 STDOUT으로 압축푼 결과를 출력하면 간단하게 명령어 하나로 끝

gunzip -c chromosome_* > sjk_ref.fa

Reference가 준비되었으면, BWA를 이용할 수 있도록 다음과 같이 reference genome에 대한 index를 생성한다.

bwa index -a bwtsw sjk_ref.fa (생성되는 파일 sjk_ref.fa.bwa, sjk_ref.fa.ann)
samtools faidx sjk.ref.fa (생성되는 파일 sjk_ref.fa.fai)

samtools의 faidx를 레퍼런스에 대한 index가 생성되며 각 컬럼은 contig이름(chr1,chr2...), contig의 길이(즉 염색체의 size), 위치(ref파일에서의 위치, 각라인당 길이, 각라인당 bytes 사이즈를 나타낸다. 이 인덱스 정보를 통해 해당 reference를 빠르게 찾아갈 수 있도록 해준다. 

자, 이제 align 수행에는 시간이 좀 걸리니까... 어떠한 결과가 나올지 궁금하겠지만... 다음 블로깅에서 만나자.