Introduction to Next Generation Sequencing data analysis

아래 글을 참고하였습니다.

http://genestack-user-tutorials.readthedocs.io/guide/intro-to-ngs.html

NGS 기술은 차세대 시퀀싱 기술이라고 말하며, DNA의 시퀀스 정보를 빠르고 값 싸게 읽어내는 기술이다. NGS 기술이란 광범위한 의미인데, 예를 들어, WGS, RNA-Seq, WES, WGBS, ChIP-Seq 등을 NGS 기술이라고 말한다. WES(Whole Exon Sequencing) 데이터의 예를 들어서 NGS 데이터 분석하는 과정에 대해서 알아보자. WES이란 DNA 상의 Exon 부위의 염기서열을 읽어내는 기술을 말한다.

우선 큰 그림을 살펴보면, 일반적인 NGS 데이터 분석은 raw reads quality control, preprocessing, mapping, post-alignment processing, variant calling, annotation, and prioritization 으로 이루어져 있다.

이제 차례대로 WES 데이터를 처리할 때 어떤 스텝을 거쳐야 하는지 알아보자. 시퀀싱 데이터를 가지고 첫 번째로 해야할일은 바로 시퀀스의 quality를 체크하는 일이다. 이를 quality control(QC)이라고 한다. 일반적으로 read의 길이와 수를 체크하게 되고 contaminating 시퀀스나 낮은 quality의 시퀀스가 있는지 찾아야한다. 이후에는 preprocessing 과정을 거치는데 이는 시퀀스의 퀄리티를 증가시키기 위한 과정이다. QC와 preprocessing 과정은 매우 중요하며, 이 과정이 제대로 되어야지 이후의 분석 결과를 신뢰할 수 있게 된다.

다음 스텝은 mapping이다. 또는 aligning이라고도 불린다. reads를 reference genome이나 reference transcriptome에 정렬하는 것을 뜻한다.  (reference genome은 표준게놈 또는 참조게놈 또는 참조 서열 등으로 불리며, HG19, HG38 등을 예로 들 수 있다.) 이렇게 함으로써 데이터의 뉴클레오타이드 서열을 de novo assembly 없이 연구할 수 있게 된다. 예를 들어 WES 데이터의 경우 WES 데이터의 read를 reference genome에 mapping 한다면, reference와 WES 데이터의 시퀀스 간의 다른 부분(variant)을 알아낼 수 있게 되고, 이 variant의 정확도는 mapping accuracy에 의존적이다.

따라서 그 다음에 해야할 일은 바로 mapping quality를 체크하는 일이다. 데이터의 특정 종류의 bias는 mapping step 이후 나타난다. mapping quality가 만족스럽지 않으면 이를 processing하는 과정을 거쳐야하는데 이를 post-alignment processing 이라고 부른다. 예를 들어, 중복된 mapped read를 제거하는 것을 의미한다. (이는 PCR의 인공물이다.) 이는 매우 중요하며 앞으로의 분석에 큰 영향을 미친다.

시퀀스가 참조 서열에 align 되었다면, 데이터를 experiment-specific한 방법으로 분석을 해야할 필요가 있다. 우리는 여기서 WES 데이터를 참조 서열에 mapping한 후에, variant 분석을 할 것이다. 즉, variant calling과 그 variant가 gene에 미치는 영향 (단백질의 변화, frame shift 등)을 알아볼 것이다. 이 과정에서 시퀀스를 참조 서열과 비교하고, 차이를 확인하여 이 차이가 유전자에 얼마나 큰 영향을 미칠지에 대해 분석하여야한다. 예를 들어, 이것이 snynonymous variant(mRNA가 생성하는 아미노산에 변화가 없는 변이)인 경우에는 이 영향이 미미할 것이다. 하지만 그 variant가 사이즈가 큰 deletion 이라면, 해당 시퀀스를 포함하는 유전자에 큰 영향을 줄 것이라고 예측해볼 수 있다. 분석과는 별개로 데이터를 visualization 해볼 수도 있다. mapped read 데이터를 visualising 함에 있어 가장 표준적인 tool 중 하나는 Genome Browser이다.

ClinVar DB 를 통한 질병 연관 변이 찾기

심근증(cardiomyopathy)과 관련있는 gene은 무엇이고, 그 gene의 어떤 위치의 variant가 pathogenic(질병의 원인이 될 수 있는)한지 알아보고자 한다. 이를 알아보기 위해서는 genotype-phenotype 데이터베이스를 이용하면되는데 이번 포스팅에서는 NCBI의 ClinVar 이라는 데이터베이스를 이용해보도록 하자. 우선 NCBI 홈페이지에 들어가 cardiomyopathy를 검색해본다.

그러면 아래와 같은 화면이 나오는데 ClinVar 을 클릭하여 들어간다.

이곳에서 왼쪽 메뉴바에서 Pathogenic을 클릭한다. 그러면 cardiomyopathy에 pathogenic한 gene의 variant를 찾을 수 있다.

또 많은 연구자들에 의해 보고된 pathgenic한 variant를 알아보고 싶다면, 아래와 같이 Multiple submitters를 체크한 후에 살펴보면 된다.

또 그 variant의 인해 분자적인 특성이 어떻게 변하는지 궁금할 수 있다. 이 때는 Molecular consequence에서 관심있는 분자적 특성에 체크를 하면된다. 예를 들어, missense(variant로 인해 mRNA의 코돈이 변해서 정상적인 경우와 다른 아미노산을 생성하는 분자적 변이)를 클릭하면 아래와 같은 화면을 볼 수 있다. 아래와 같이 다양한 pathogenic한 variant와 해당 gene을 알 수 있다. 해당 gene의 variant로 인해 missense가 일어나고 이로 인해 cardiomyopathy 질병을 일으키는 것이다.

크게보면 아래와 같이 variant를 나타내고 있다. 이렇게 variant를 표기하는 방법을 HGVS notation이라고 한다.

이를 통해 검색할 수도 있는데 예를 들어,

NM_213653.3(HFE2):c.959G>T (p.Gly320Val)

를 ClinVar 검색창에 검색해보자.

그러면 위와 같은 화면을 볼 수 있다. Links를 보면 해당 Variant가 UniProtKB, OMIM, dbSNP, 100 Genome Browser에 등록되어 있다는 것을 알 수 있다. 이를 통해 이 variant가 Multiple Submitters에 submit 되었다는 것을 알 수 있다. (Number of Submissions에 4가 표시되어 있다.)

ClinVar 데이터 다운로드 받기

https://www.ncbi.nlm.nih.gov/clinvar/에 접속 후 FTP => Go to the FTP site

여기서 vcf_GRCh37, vcf_GRCh38이 variant를 vcf 파일 형태로 보여주는 것이고, xml은 xml 파일 형태로 보여주는 것이다.

들어가면 다양한 버전이 존재하는데 가장 위에 있는 것을 받으면 된다. 

CIGAR Format

출처 - https://genome.sph.umich.edu/wiki/SAM

/* 2017.8.13 by 3months. */

CIAGR Format은 어떠한 Sequence A가 reference genome에 align되어있을 때 이를 나타내 주는 포맷이다.

예를 들어, 아래와 같은 reference sequence와 Read가 있다고 하자. 이 Read를 reference sequence에 align 해보면

RefPos:     1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T  G  A  A  C  T  G  A  C  T  A  A  C
Read: ACTAGAATGGCT

이와 같이 align 된다.

RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T

이 때 이것을 나타내는 CIGAR Format은

POS: 5
CIGAR: 3M1I3M1D5M

이렇게 되는데 reference의 5번 position부터 시작해서 3개의 Matching, 1개의 insertion, 다시 3개 matching, 1개 deletion, 5개 matching 이라는 뜻이다. ref의 14번 position에 read와 시퀀스가 다르지만, align이 이렇게 되기 때문에 matching으로 쳐준다.

Lastz를 통한 CIGAR Format 생성

lastz_test.fa

lastz_test2.fa

의도적으로 lastz_test fasta 파일의 sequence에다가 4개의 sequence를 insert하고 2개의 sequence를 delete해서 lastz_test2.fa 파일을 생성하였다.

그 다음

lastz lastz_test.fa lastz_test2.fa --format=general:cigar > cigar.txt

를 실행 후, cigar.txt 파일을 열어보면

#cigar
208M4I163M2D138M

위와 같은 결과가 나왔음을 볼 수 있다. 결과를 해석하면 208개 matching 4개 insertion, 163개 matching 2개 deletion, 138개 matchin 이다.

2017.10.8 추가

lastz 사용해서 cigar을 얻을 때 다른 옵션을 주면 substitution도 고려가 가능하다.

lastz_test.fa

lastz_test2.fa

Lastz 실행

/* 2017.8.13 by. 3months */

Lastz를 설치하였으면 간단하게 결과를 뽑아보자.

두 개의 fasta 파일을 폴더에 붙여넣고,

해당 폴더에서 아래 명령어를 입력해보자.

lastz lastz_test.fa lastz_test2.fa

결과

그러면 콘솔창에 이런 결과가 나온다.

또는 아래와 같이 옵션을 주어서 실행하면 maf 파일로 내보낼 수 있다.

lastz lastz_test.fa lastz_test2.fa --format=maf > lastz_result.maf

lastz_result.maf 파일을 열어보면 아래와 같이 되어있다.  (이것이 무엇을 뜻하는지는 다음에 알아보자..)

##maf version=1 scoring=lastz.v1.04.00
# lastz.v1.04.00 --format=maf
#
# hsp_threshold      = 3000
# gapped_threshold   = 3000
# x_drop             = 910
# y_drop             = 9400
# gap_open_penalty   = 400
# gap_extend_penalty = 30
#        A    C    G    T
#   A   91 -114  -31 -123
#   C -114  100 -125  -31
#   G  -31 -125  100 -114
#   T -123  -31 -114   91

a score=47093
s ref 0 511 + 511 GAAGAGCAGCACACTCAGATCTTCTTACACCATCTCTGGGGGCCGACTGGCAGAGAAGGCTAGGGACTTGGCCAATATCATGCAGCCTTCATGCCTTATTCGGGAAAGGTGCTCCAAATAGGAGGAATATTTGGAAGAAATCCCAATAGGGATTTGGAGACATGCCATTTGGGAATGCAAACCCAGCCCCTCCT--------AGTGCTCCACCCTG----CCCATGTCCCCATTGAACCGACAACCTGCAGAGAGAGCAGGATGCCCACAGCAAGCCCCAAGGACAAAGGCCACAGAGTCAGGGCACCCAAATCCTAACCTAACCTTGCCTGTGCTTCCTGTGTACTGTGGCCATTTTACTTTCCCCCTGTTCCCTCCTCTGTAAATGAGGCCCTTGGACTGCATCAGTGGTTTTTACACTGAGCTCCCTGGAGCTCTGGGGGCCTCAGGGACTATTGTGGAGGCCTTGCTGCATTTGAGAGTGGTGTCTTGTACCCTACTCAACCCATTTTTTCTGTTTTAC

s ref 0 521 + 521 GAAGAGCAGCACACTCAGATCTTCTTACACCATCTCTGGGGGCCGACTGGCAGAGAAGGCTAGGGACTTGGCCAATATCATGCAGCCTTCATGCCTTATTCGGGAAAGGTGCTCCAAATAGGAGGAATATTTGGAAGAAATCCCAATAGGGATTTGGAGACATGCCATTTGGGAATGCAAACCCAGCCCCTCCTAGAGAGAGAGTGCTCCACCCTGCCCCCCCATGTCCCCATTGAACCGACAACCTGCAGAGAGAGCAGGATGCCCACAGCAAGCCCCAAGGACAAAGGCCACAGAGTCAGGGCACCCAAATCCTAACCTAACCTTGCCTGTGCTTCCTGTGTACTGTGGCCATTTTACTTTCCCCCTGTTCCCTCCTCTGT--ATGAGGCCCTTGGACTGCATCAGTGGTTTTTACACTGAGCTCCCTGGAGCTCTGGGGGCCTCAGGGACTATTGTGGAGGCCTTGCTGCATTTGAGAGTGGTGTCTTGTACCCTACTCAACCCATTTTTTCTGTTTTAC

아래 링크에서 lastz 실행과 관련된 다양한 예제를 볼 수 있다.

http://www.bx.psu.edu/miller_lab/dist/README.lastz-1.02.00/README.lastz-1.02.00a.html#examples

Lastz 설치

Lastz는 Large-Scale Genome Alignment Tool의 약자로 말그대로 큰 genome 시퀀스를 alignment하는 툴이다.

이 링크에 들어간 후 https://github.com/lastz/lastz

프로젝트 파일을 git을 통해 clone하거나 zip파일로 다운받은 후, 압축을 해제하여 README.lastz.html 이 파일을 웹 브라우저에서 연다.

이런 느낌의 html 문서가 나오는데 Installation 파트에 가면 설치 과정을 자세히 볼 수 있다.

다른 생명정보학 툴과 마찬가지로 리눅스 계열 운영체제에서 동작한다.

lastz\src 폴더에 들어가서

make

make install

이 명령어를 통해 두 개의 실행가능한 프로그램이 운영체제에 설치가된다. lastz, lastz_D인데 차이는 lastz는 integer score, lastz_D는 floating point score를 쓴다고 한다.

아래 커맨드를 입력했을 때 아웃풋이 나오지 않으면 설치가 제대로 된 것이다.

make test

그리고

./lastz

를 통해 lastz 프로그램이 동작하는지 확인해보자.

편하게 lastz 프로그램을 사용하기 위해 환경변수를 등록한다.

export PATH=[설치된경로]/lastz/src:$PATH

이제 어느곳에서나 lastz 명령어를 통해 프로그램을 실행할 수 있다.

아래는 lastz 명령어 실행시 나타나는 메시지. 타겟 파일 이름을 지정하라고 나온다.

You must specify a target file
lastz-- Local Alignment Search Tool, blastZ-like
  (version 1.04.00 released 20170312)
usage: lastz target [query] [options]
  (common options;  use --help for a more extensive list)
  target, query          specifiers or files, containing sequences to align
                         (use --help=files for more details)
......

Platypus 설치 확인

platypus 설치하기 링크 : 3months.tistory.com/155

platypus를 설치하였으면 아래 명령어를 통해 variant calling을 하여 플라티푸스가 잘 설치되었는지 확인해보자.

python Platypus.py callVariants --bamFiles=input.bam --refFile=ref.fa --output=VariantCalls.vcf

(이 명령어를 실행하기 전에 전에 samtools를 활용해 bam 파일과 fasta 파일의 인덱스 파일을 만들어 놓아야한다.)

실행 에러

간단한 sample 데이터로 위 명령어를 실행하니

ImportError: libhts.so.2: cannot open shared object file: No such file or directory

라는 에러메시지가 출력되면서 실행되지 않았다. 찾아보니 htslib의 라이브러리 경로가 제대로 지정되지 않았기 때문이라고 한다.

https://github.com/andyrimmer/Platypus

위 링크를 보고 해결하였는데 htslib를 설치하면 /usr/local 폴더에 htslib 가 설치되는데 여기서 /usr/local/lib와 /usr/local/include 폴더의 라이브러리 경로를 잡아주어야 한다.

export C_INCLUDE_PATH=/usr/local/include
export LIBRARY_PATH= /usr/local/lib
export LD_LIBRARY_PATH=/usr/local/lib

bash 에서 위 명령어를 실행하여 라이브러리 경로를 잡아주고 다시 실행하니 잘 실행이 되는 것을 볼 수 있다.

그러면 실제 데이터로 variant calling을 해보자.

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq/

위 링크에서 아무거나 골라서 bam 파일과 인덱스 파일을 다운받는다.

그리고 hg19의 chromosome1 fasta file을 다운받는다. http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/

그리고 chromosome1의 variant 만 확인하기 위해 bam 파일에서 chromosome1의 시퀀스만 뽑아낸다.

samtools view -b sample.bam chr1 > sample_ch1.bam

위 명령어를 통해 bam 파일에서 ch1 부분만 뽑아낼 수 있다.  그리고 이 bam 파일에서 인덱스 파일을 생성한다.

samtools index sample_ch1.bam

그러면 sample_ch1.bam.bai 파일이 생성되며 준비가 끝났다.

또한 hg19의 chr1 fasta 파일도 인덱스 파일을 생성해준다.

samtools faidx chr1.fasta

이제 플라티푸스 명령어를 돌리기 위한 준비가 끝났다.

python Platypus.py callVariants --bamFiles=sample_ch1.bam --refFile=chr1.fa --output=VariantCalls.vcf

실행되는 모습

생성된 VariantCalls.vcf 파일

bam 파일 인덱싱

BAM 파일은 FASTA 파일처럼 sequence 를 저장하는 파일인데 binary로 저장이된다. (저장공간의 효율을 위해 그런것같다.) 시퀀스를 저장하고 있는 모양은 SAM 파일과 같은데, 이를 단지 binary 형태로 저장한 것이다.

bam 파일 인덱싱 관련 정보는 https://www.biostars.org/p/15847/ 이 곳에서 얻었는데, bam 파일은 companion file 인 index 파일을 갖을 수 있다. 이것은 외부에 저장되는 bam 파일의 인덱스 테이블이다. 이 테이블을 통해 프로그램이 bam 파일의 원하는 위치로 점프할 수 있다.

왜 platypus 파일에서 variant calling 할 때 인덱싱된 bam 파일을 제공하라고 하는지 알 것 같다. 인덱스 테이블을 통해 원하는 시퀀스로 점프할 수 있게 하는것 이다. 이 index 파일은 bam 파일과 똑같은 이름에 끝에만 .bai 를 붙인다. (bam 파일 이름이 ch1.bam 이면 인덱스 파일은 ch1.bam.bai 이다. 이러한 naming을 통해 해당 인덱스 파일이 어느 bam 파일로부터 나왔는지를 알 수 있다.

bam 파일 인덱싱하는 법은 쉬운데 samtools 설치 후

samtools index bamfile.bam

를 실행하면 된다.

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq/

이곳에서 bam 파일을 아무거나 다운 받은 후 테스트 해보자. 제일 위에 있는 것을 다운 받은 후 파일명을 bamfile.bam으로 바꾸었다. 근데 위 command를 실행하니 오류가 떴다. 검색해 보니 먼저 bam파일을 sorting을 한 후 indexing 을 해야하는 것 같다.

samtools sort bamfile.bam

samtools index bamfile.bam

를 실행한다. 참고로 samtools 명령어 모음집은 아래 링크에서 찾을 수 있다.

http://www.htslib.org/doc/samtools.html

/* 2017.8.7 */

공부와 메모를 위한 포스팅입니다.

platypus 설치 링크

http://www.well.ox.ac.uk/platypus

https://github.com/andyrimmer/Platypus

위의 매뉴얼대로하면 된다.

압축파일을 다운받고 빌드하면 설치 완료.

또한 두 번째 github 링크 주소에서 나오듯이 htslib를 설치하여야 한다. (dependancy가 있다.)

platypus는 haplotype variant calling 하는 프로그램.

command line에서 실행할 때는 파이썬 프로그램을 command line에서 실행하는 것처럼 하면 된다.

platypus 실행

가장 간단하게 platypus를 실행하는 방법은

python Platypus.py callVariants --bamFiles=input.bam --refFile=ref.fa --output=VariantCalls.vcf

플라티푸스 폴더에 들어가서 저 명령어를 실행하는 것인데,

이 때, 1개 이상의 BAM 파일을 인풋으로 제공할 수 있으며, reference file은 fast파일로 제공하면된다.

또한 BAM 파일은 samtools나 이와 비슷한 프로그램으로 indexing이 되어있어야하며,

FASTA 파일도 마찬가지로 "samtools faidx" 나 이와 비슷한 프로그램으로 indexing 되어있어야한다.

아웃풋은 1개의 .vcf 파일과 log 파일이다.

log 파일의 최하단에 'Finished variant calling' 메시지를 보고 제대로 작동하였음을 확인할 수 있다.

그러면 이제 samtools로 bam 파일과 fasta 파일을 indexing 하는 것을 해보자.

http://www.htslib.org/download/

위 웹사이트에서 samtools를 다운받고 설치하는 도중에 에러가 출력되면서 설치가 제대로 되지 않는 문제가 있었다.

문제

cd htslib-1.5 && make lib-static
make[1]: Entering directory `/mnt/c/samtools-1.5/htslib-1.5'
gcc -g -Wall -O2 -I.  -c -o cram/cram_io.o cram/cram_io.c
cram/cram_io.c:60:18: fatal error: lzma.h: 그런 파일이나 디렉터리가 없습니다
 #include <lzma.h>
                  ^
compilation terminated.
make[1]: *** [cram/cram_io.o] 오류 1
make[1]: Leaving directory `/mnt/c/samtools-1.5/htslib-1.5'
make: *** [htslib-1.5/libhts.a] 오류 2

해결

apt-get install -y liblzma-dev

윈도우 1주년 업데이트(레드스톤 버전)부터 지원되는 윈도우 bash shell 기능을 통해 윈도우 우분투에 gatk를 설치하여보겠습니다.

1. 윈도우 커맨드를 열고 "bash" 명령어 입력을 통해 윈도우 bash shell 실행

윈도우 우분투를 이용하는 방법은 https://blogs.msdn.microsoft.com/eva/?p=7633 참고해주세요.

2. oracle java 1.8 버전 설치

sudo add-apt-repository ppa:webupd8team/java
sudo apt-get update
sudo apt-get install oracle-java8-installer

설치 후 java -version을 입력하고 아래와 같은 문구가 나오면 제대로 설치된 것입니다.

java version "1.8.0_131"
Java(TM) SE Runtime Environment (build 1.8.0_131-b11)
Java HotSpot(TM) 64-Bit Server VM (build 25.131-b11, mixed mode)

3. GATK jar 파일 다운로드

https://software.broadinstitute.org/gatk/download/

4. 위의 jav 파일이 있는 곳으로 cd를 통해 이동 후, 아래 명령어 입력

java -jar GenomeAnalysisTK.jar -h