GATK4 (Genome Analysis Tool Kit)로 넘어오면서, NGS 시퀀싱 분석을 위한 파이프라인이 많이 개선 및 간소화된 것 같습니다. 덕분에 저도 최근 연구실에 구축되어 있던 파이프라인도 새롭게 뜯어고쳤는데, 이 참에 전반적인 분석을 위한 코드를 정리해볼까 합니다. 추후에 비슷한 파이프 라인을 구축하고자 하는 분들께 도움이 되었으면 하고 피드백도 환영합니다.
관련 포스팅 보기>
I. 들어가기에 앞서 (준비물)
우선 기본적으로 모든 작업은 리눅스 (linux) 환경에서 이루어집니다. 리눅스 환경에서 기본 명령어를 숙지하고 필요한 프로그램들을 설치합니다. 기본적으로 아래와 같은 프로그램들을 다운로드 및 설치합니다. 자세한 설치 방법의 경우, 해당 웹사이트에서 자세히 제공하고 있습니다. 저의 경우, 현재 최신버젼인 BWA 0.7.5, PICARD 2.20, GATK 4.1.3을 기준으로 설치 및 작업을 진행하였습니다.
BWA (Burrows-Wheeler Aligner): http://bio-bwa.sourceforge.net/
Samtools: http://www.htslib.org
위의 프로그램이 설치가 끝나면, path 등록을 해서 리눅스 상에서 명령어를 입력하면 잘 실행되는지 확인합니다. 마지막으로 Reference Genome 파일을 다운로드 합니다. hg19 또는 hg38 버젼에 맞게 파일을 모두 다운로드하여 위치시키면 기본적인 준비들이 끝납니다. (파일은 UCSC genome browser에서 다운 가능)
II. Data pre-processing: Fastq to Bam
아래 코드는 GATK best practice의 Data pre-processing step에 해당하는 내용을 정리한 것입니다. 전체 Work-flow는 아래와 같습니다.
GATK Best practice: Data pre-processing for variant discovery
1. Map to Reference: Fastq to Bam
bwa index [hg19_reference.fa] bwa mem -t 8 -M -R '@RG\tID:foo\tSM:bar\tLB:library1' [hg19_reference.fa] [sample01_R1.fq.gz] [sample01_R2.fq.gz] > [sample01.sam] samtools view -bt [hg19_reference.fa.fai] [sample01.sam] -o [sample01.bam] samtools sort [sample01.bam] -o [sample01.sorted.bam] samtools index [sample01.sorted.sam]
Reference sequence를 다운로드 받았으면, 기본적으로 indexing을 한번 진행해 준다. 이후, 위의 명령어는 Paired-end Read의 2가지 fastq 파일을 1개의 Sam 파일로 mapping하는데 필요한 명령어이다. []안의 부분을 적당히 변경하여 이용할 수 있다. @RG에 해당하는 부분은 기본적인 read에 붙은 바코드 정보를 말한다. 실제 명령어에서는 대괄호는 [] 입력하지 않는다. Fastq 파일은 압축 형태를 바로 이용할 수 있으므로, 압축을 모두 풀 필요는 없다. -t 는 thread option으로 컴퓨터 성능에 따라 적절히 지정해준다. samtools를 이용하여, sam > bam 변환 및 동시에 sorting, indexing을 진행한다.
2. MarkDuplicate + SortSam
java -jar picard.jar MarkDuplicates I=[sample01.sorted.bam] O=[sample01.dedup.bam] M=[sample01.markdup.metrics.txt] java -jar picard.jar SortSam I=[sample01.dedup.bam] O=[sample01.sortsam.bam]
MarkDuplicate와 SortSam을 이용하여 중복되는 Read를 걸러주고, Read를 정렬시켜서 이후에 Variant call을 위한 기본 작업을 진행한다. 이 부분은 과거 Samtools를 이용하여, 주로 이용하였는데, SortSam으로 대체해서 진행이 가능하다.
3. Base Quality Score Recalibration (BQSR)
gatk --java-options 'Xmx16g' 'Xms8g' BaseRecalibrator -I [sample01.sortsam.bam] -R [hg19_reference.fa] --known-sites [dbsnp_151.GRCh37.vcf] -O [sample01.recal_data.table] gatk --java-options 'Xmx16g' 'Xms8g' ApplyBQSR -I [sample01.sortsam.bam] -R [hg19_reference.fa] --bqsr-recal-file [sample01.recal_data.table] -O [sample01.final.bam]
BQSR는 시퀀싱을 통해 생성된 염기들의 Quality Score에 발생하는 bias를 보정해서, 잘못된 mapping을 보정해주는 QC 과정의 일종이다. 위의 명령어를 성공적으로 실행시키면, 최종적으로 분석 가능한 Bam 파일이 생성된다. 생성된 Bam 파일은 이후에 변이를 Call하여 VCF 파일을 생성하는데 이용된다.
[References]
GATK Best practice: Data pre-processing for variant discovery
NCI GDC DNA-Seq Analysis pipeline: https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/DNA_Seq_Variant_Calling_Pipeline/
Samtools (Mapping to Reference): http://www.htslib.org/workflow/#mapping_to_variant
두마디 정밀의료 
대단하십니다
좋아요좋아요
직접 해보시면 별로 어렵지 않습니다 ^^;
좋아요좋아요
안녕하세요. 스텝별로 따라하고 있었습니다. 좋은 자료 올려주셔서 감사합니다. 다름이 아니라 BQSR 명령어를 치면
-bash: gatk: command not found
라고 나오는데 에러를 검색해봐도 답을 찾지 못해 여쭤봅니다.
그리고
–known-sites [dbsnp_151.GRCh37.vcf]
에서 dbsnp_151.GRCh37.vcf 파일을 다운로드받아야 하는 것인지요? 이 부분 없이 하면 의미가 없는지도 여쭙습니다. 감사합니다.
좋아요좋아요
1) bash: gatk: command not found > gatk 가 제대로 설치가 안된 것 같습니다. gatk라고 쳤을 때, gatk가 돌아갈 수 있도록 설치부터 하셔야 합니다. 자세한 설치 과정은 gatk 홈페이지를 참고 하세요. 2) dbsnp database를 따로 다운 받아야 합니다. 구글링 하시면 dbSNP database에서 다운 받으실 수 있습니다. 말그대로 알려진 site에서 training 시키는 DB이기 때문에, 다른 DB를 이용해도 되나, 필터링 모델의 정확도가 떨어지게 됩니다. 얼마나 차이가 나냐고 물으시면 저도 안해봐서 잘 모르겠습니다.
좋아요좋아요