[WGS 실습] 변이 호출 과정 검증

검증은 생성한 VCF 파일을 GIAB truth set이랑 비교함.

 

검증에는 rtg-tools 이라는 툴을 활용하는데 설치를 먼저 진행.

conda install -c bioconda rtg-tools

 

잘 설치됐는지 확인.

rtg --version

---

검증

검증 실행 전 truth VCF 인덱스 생성

tabix -p vcf ~/bioinfo/bi_study/wgs_practice/truth/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz

 

위 명령어를 통해 인덱스를 생성하고 밑과 같이 tbi 파일이 생성된 것을 확인.

ls ~/bioinfo/bi_study/wgs_practice/truth/

---

검증 실행 전 에러 발생을 예방하기 위한 reference SDF 형식 변환

rtg format -o ~/bioinfo/bi_study/wgs_practice/ref/chr20.sdf ~/bioinfo/bi_study/wgs_practice/ref/chr20.fa

rtg-tools는 reference를 SDF 형식으로 변환해야 하기 때문에 SDF 형식으로 변환

---

검증 실행

rtg vcfeval \
  -b ~/bioinfo/bi_study/wgs_practice/truth/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz \
  -c ~/bioinfo/bi_study/wgs_practice/variants/SRR062634.vcf.gz \
  -t ~/bioinfo/bi_study/wgs_practice/ref/chr20.sdf \
  -o ~/bioinfo/bi_study/wgs_practice/variants/eval_output

옵션 설명

 

-b : baseline.

truth set VCF. 정답지 역할.

-c : calls.

우리가 만든 VCF. 채점 대상.

-t : template.
reference genome.

-o : output 폴더. 결과 파일들이 여기 저장됨.

---

검증 결과 확인

Selected score threshold using: maximized F-measure
Threshold  True-pos-baseline  True-pos-call  False-pos  False-neg  Precision  Sensitivity  F-measure
----------------------------------------------------------------------------------------------------
    3.000                149            149      28056      79558     0.0053       0.0019     0.0028
     None                149            149      28078      79558     0.0053       0.0019     0.0028

 

결과를 확인했을 때 수치가 낮지만, 결과를 확인해보면

항목	값	의미
True-pos	149	우리가 맞게 찾은 변이
False-pos	28,056	우리가 찾았는데 truth에 없는 것 (과검출)
False-neg	79,558	truth에 있는데 우리가 못 찾은 것
Precision	0.0053	우리가 찾은 것 중 진짜 비율
Sensitivity	0.0019	truth 중 우리가 찾은 비율
F-measure	0.0028	Precision + Sensitivity 종합 점수

해석

1. 실습에서는 chr20 전체가 아니라 100만 read subset만 써서 커버리지가 매우 낮음.

2. truth set은 chr20 전체 변이를 담고 있어서 실습에서 커버 못 한 위치의 변이가 False-neg로 잡히고 있음.

3. 필터링을 안 해서 False-pos가 많음. 저품질 변이들이 다 포함된 원인.

 

실제 WGS 전체 데이터로 제대로 돌리면 F-measure 0.99 이상 나오지만, 파이프라인 흐름을 이어나가는 실습이다보니 결과 자체에 큰 비중을 두지는 않음.

---

이번 실습동안 진행한 파이프라인 흐름

Raw FASTQ (SRR062634_1.fastq.gz, _2.fastq.gz)
         ↓
    [FastQC] - 품질 확인
         ↓
  [Trimmomatic] - 저품질 염기 제거
         ↓
Trimmed FASTQ (_1_paired.fastq.gz, _2_paired.fastq.gz)
         ↓
    [FastQC] - Trimming 후 품질 재확인
         ↓
  [bwa-mem2 mem] - reference(chr20.fa)에 정렬
         ↓
      SAM 파일
         ↓
  [samtools sort] - 좌표 기준 정렬 + BAM 변환
         ↓
      BAM 파일
         ↓
  [samtools index] - BAM 인덱스 생성
         ↓
  [GATK AddOrReplaceReadGroups] - Read Group 추가
         ↓
   BAM (RG 포함)
         ↓
  [GATK HaplotypeCaller] - 변이 호출
         ↓
      VCF 파일
         ↓
  [rtg vcfeval] - GIAB truth set과 비교 검증
         ↓
   TP / FP / FN / F-measure