NGS 분석 기초 3 - NGS데이터 분석. Variants calling(GATK3.8활용)

공부하면서 정리하는 Bioinformatics

 

본 포스팅은 genomics를 공부하며 작성하는 포스팅이기 때문에 잘못된 부분이 있을수 있습니다.

또한, 범문에듀케이션에 출판된 유전체 데이터 분석2 (NGS편, 암과 질병 유전체) 서적을 기반으로 공부하여 작성하였음을 미리 알립니다.

---

 NGS 데이터 처리 단계

 

    1) reference sequence indexing(BWA)

    2) reference sequence에 reads alignment(BWA)

    3) SMA -> BAM (samtools)

    4) sorting BAM (Picard)

    5) 중복 reads 제거 (Picard)

    6) variants calling에 사용될 .dict, .fai 파일 생성 (Picard)

 

지난 포스팅에서 variants calling을 할 준비를 마쳤다.

이어서 변이들을(SNP/InDel) calling하는 과정을 공부한다.

 

7) InDel 분석을 위한 indel realignment

 

-Create intervals

$$JAVA -Xmx1024m -jar $GATK \
> -R $REFERENCE_GENOME \
> -T RealignerTargetCreator \
> -known Mills_and_1000G_gold_standard.indels.hg.chr21.vcf \
> -I dedup_sorted_ex.bam \
> -o realigner.intervals

#-R : 참조유전체서열 파일명
#-T : GATK내의 tool, interval 파일생성을 위해 RealignerTargetCreator사용
#-known : 앞서 알려진 indel data
#-I : input
#-o : output

-Realigninig

 

$$JAVA -Xmx1024m -jar $GATK \
> -T IndelRealigner\
> -R $REFERENCE_GENOME \
> -I dedup_sorted_ex.bam \
> -targetIntervals realigner.intervals \
> -o realign_dedup_sorted_ee.bam

#-R : 참조유전체서열 파일명
#-T : GATK내의 tool, interval 파일생성을 위해 IndelRealigner
#-targetIntervals : RealignerTargetCreator로 생성된 interval파일 이름
#-I : input
#-o : output

 

8) Quality Score Recalibration

시퀀싱 장비에서 얻은 염기서열 퀄리티는 실질적인 sequancing calling error rate를 잘 반영하지 못함.

그래서 리드들의 데이터를 이용하여 시퀀싱 장비가 출력한 phred quality score를 보정해준다.

 

-BaseRecalibrator

$$JAVA -Xmx1024m -jar $GATK \
> -T BaseRecalibrator \
> -R $REFERENCE_GENOME \
> -I realign_dedup_sorted_ex.bam \
> -knownSites dbsnp_138.hg19.chr21.vcf \
> -o recall.grp

#-T: 분석타입
#-knownSites : 알려진 SNP파일 (from dbSNP138)
#-o : recalibration table이름

-recal.grp 테이블을 이용해 스코어보정 후 recalibration table 작성

$$JAVA -Xmx1024m -jar $GATK \
> -T BaseRecalibrator \
> -R $REFERENCE_GENOME \
> -I realign_dedup_sorted_ex.bam \
> -knownSites dbsnp_138.hg19.chr21.vcf \
> -BSQR recal.grp \
> -o post_recal.grp

#-BSQR : recalibration table 이름
#-o : recalibration 후 보정된 recalibration table

-Recalibration 평가 그래프

$$JAVA -Xmx1024m -jar $GATK \
> -T AnalyzeCovariates \
> -R $REFERENCE_GENOME \
> -befor recal.grp \ 
> -after posst_recal.grp\
> -plots recal_plots.pdf

#-befor : 보정 전 계산된 recalibration table
#-after : 보정 후 계산된 recalibration table
#-plots : recalibration 전/후 퀄리티 그래프 생성

-Recalibration 정보를 BAM에 쓰기

보정된 테이블을 이용하여 BAM파일의 퀄리티 스코어를 다시 작성.

$$JAVA -Xmx1024m -jar $GATK \
> -T PrintReads \
> -R $REFERENCE_GENOME \
> -I realign_dedup_sorted_ex.bam \
> -BSQR post_recal.grp \
> --out recal_realign_dedup_sorted_ex.bam

9) SNP calling

snp, indel 추추을 위해서는 두가지 콜러를 사용할수있다.

1. HaplotypeCaller - de-novo assembly 알고리즘 기반이라 성능이좋지만 계산속도가 느림.

#haplotypecaller를 이용한 SNP/indel Calling
$$JAVA -Xmx1024m -jar $GATK \
> -R $REFERENCE_GENOME \
> -T HaplotypeCaller \
> -I recal_realign_dedup_sorted_ex.bam \
> -o gatk_raw_hc.vcf\
> --dbsnp dbsnp_138.hg19.chr21.vcf \
> -stand_call_conf 30.0

#--dbsnp : vcf파일의 ID 컬럼에 해당하는 부분에 주석을 달기위한 용도로만 사용됨.
#-stand_call_conf 30.0 : high confidence를 구분하는 최소 phred-scaled Qscore threshold값

2. UnifiedGenotyper - Baysian genotype likelihood model기반. 샘플수가 많은 경우 사용 권장.

#UnifiedGenotyper를 이용한 SNP/indel Calling
$$JAVA -Xmx1024m -jar $GATK \
> -R $REFERENCE_GENOME \
> -T UnifiedGenotyper \
> -I recal_realign_dedup_sorted_ex.bam \
> -o gatk_raw_ug.vcf\
> --dbsnp dbsnp_138.hg19.chr21.vcf \
> -stand_call_conf 30.0 \
> -glm BOTH

#--dbsnp : vcf파일의 ID 컬럼에 해당하는 부분에 주석을 달기위한 용도로만 사용됨.
#-stand_call_conf 30.0 : high confidence를 구분하는 최소 phred-scaled Qscore threshold값
#-glm BOTH : SNP/InDel 모두 추출한다는 뜻. SNP,INDEL,BOTH

 

-SelectedVariants

UnifiedGenotyper로 생성한 gatk_raw_ug.vcf파일중 SNV 또는 InDel로 이뤄진 variant들을 선별하여 추출

#UnifiedGenotyper를 이용한 SNP/indel Calling
$$JAVA -Xmx1024m -jar $GATK \
> -R $REFERENCE_GENOME \
> -T SelectVariants \
> --variant gatk_raw_ug.vcf \
> -o snp_gatk_raw_ug.vcf\ # indel_gatk_raw_ug.vcf
> -selectType SNP #or INDEL 선택

-InDel Filtering

$$JAVA -Xmx1024m -jar $GATK \
> -R $REFERENCE_GENOME \
> -T VariantFiltration \
> --variant indel_gatk_raw_ug.vcf \
> -o filtered_indel_gatk_raw_ug.vcf \ 
> --filterExpression "QD < 2.0" --filterName "QD2" \
> --filterExpression "ReadPosRankSum < 20.0" --filterName "ReadPosRankSum-20" \
> --filterExpression "FS > 200.0" --filterName "FS200" \

-SNP Filtering

$$JAVA -Xmx1024m -jar $GATK \
> -R $REFERENCE_GENOME \
> -T VariantFiltration \
> --variant snp_gatk_raw_ug.vcf \
> -o filtered_snp_gatk_raw_ug.vcf \ 
> --filterExpression "QD < 2.0" --filterName "QD2" \
> --filterExpression "MQ < 40.0" --filterName "MQ40" \
> --filterExpression "FS > 60.0" --filterName "FS60" \

# --filterExpression "QD < 2.0" --filterName "QD2" : Quality by Depth(QUAL/DP)값이 2미만으로
# --filterExpression "ReadPosRankSum < 20.0" : alternative base가 리드의 시작이나 끝에 왜곡되있으면 오류로 간주
# --filterExpression "MQ < 40.0" :stand bias 검출
# --filterExpression "FS > 60.0" : 해당 ALT SNV에서 median quality값이 40미만인경우

-필터가 PASS인 변이 추출

$$JAVA -Xmx1024m -jar $GATK \
> -R $REFERENCE_GENOME \
> -T SelectVariants \
> --variant filtered_SNP_gatk_raw_ug.vcf \
> -select 'vc.isNotFiltered()' \
> -o pass_filtered_SNP_gatk_raw_ug.vcf \ 

#Indel도 마찬가지
# -select 'vc.isNotFiltered()': filter가 pass인 variant추출

 

이로써 GATK3.8을 활용하여 FASTQ파일에서 SNV와 indel을 calling하는 모든 단계를 진행했다.

필터는 입맛대로 옵션값을 바꿀수 있다.