VCF 파일 파싱 및 주석 달기

이 페이지는 다음 방법에 대해 설명합니다.

  • 입력 VCF 파일에서 주석 필드를 파싱합니다.
  • 오픈소스 주석 도구와 Variant Effect Predictor(VEP) 같은 데이터베이스를 사용하여 VCF 파일을 가져오는 동안 주석을 추가합니다.

주석 필드 파싱

Variant Transforms 도구는 VCF 형식의 변이 주석에 지정된 주석 필드 형식을 지원합니다. VEP도 이 형식을 사용합니다. 다음 예시는 VEP에 의해 주석이 추가된 파일의 VCF 헤더를 보여줍니다.

##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|ALLELE_NUM|...>

REF=AACALT=TAC,A인 샘플 변이 레코드의 해당 주석은 다음과 같습니다.

CSQ=-|downstream_gene_variant|MODIFIER|WASH7P|ENSG00000227232|Transcript|ENST00000423562|unprocessed_pseudogene|||2,...,T|upstream_gene_variant|MODIFIER|DDX11L1|ENSG00000223972|Transcript|ENST00000518655|transcribed_unprocessed_pseudogene|||1

각 변이 레코드에 대해 여러 개의 쉼표로 구분된 주석 목록이 제공됩니다. 각 집합의 첫 번째 필드는 대체 대립유전자(ALT)에 해당합니다. 위 예시에는 각 ALT에 대해 하위 집합만 사용됩니다. 실제 파일은 필드가 더 있습니다. 그리고 다른 주석을 사용해서 각 대체 대립유전자를 여러 번 반복할 수 있습니다.

Variant Transforms 도구를 실행할 때 --annotation_fields를 전달하여 주석 목록을 분할하고 해당 alternate_bases 레코드 아래로 이동할 수 있습니다. 예를 들어 아래에서는 각 alternate_bases 아래에 CSQ 필드가 반복되는 것을 확인할 수 있습니다.

이름 유형 모드 설명
alternate_bases.CSQ 레코드 반복 이 대체에 대한 CSQ 주석의 목록입니다.
alternate_bases.CSQ.allele 문자열 null 허용 주석 필드의 ALT 부분입니다.
alternate_bases.CSQ.Consequence 문자열 Null 허용 없음
alternate_bases.CSQ.IMPACT 문자열 Null 허용 없음
alternate_bases.CSQ.SYMBOL 문자열 Null 허용 없음
alternate_bases.CSQ.Gene 문자열 Null 허용 없음
alternate_bases.CSQ.Feature_type 문자열 Null 허용 없음
alternate_bases.CSQ.Feature 문자열 Null 허용 없음
alternate_bases.CSQ.BIOTYPE 문자열 Null 허용 없음
alternate_bases.CSQ.EXON 문자열 Null 허용 없음
alternate_bases.CSQ.INTRON 문자열 Null 허용 없음

대체 대립유전자 일치

위의 예시에서 CSQ 열에 나열된 대체 대립유전자는 변이 레코드의 ALT 문자열과 다릅니다. 예를 들어 AAC->C 변이는 삭제를 의미하는 -로 표시됩니다. 이러한 상황은 VEP가 --minimal 플래그로 데이터세트에서 실행되기 때문에 발생합니다. Variant Transforms는 이 minimal 모드에서 ALT 필드를 일치시키는 것을 지원합니다.

다음과 같은 ALT 일치 모드가 지원됩니다.

  • 기본적으로 ALT 일치 플래그가 설정되지 않은 경우 VCF 형식의 변이 주석 사양이 사용됩니다.

  • 위 예시에서와 같이 --use_allele_num 플래그를 사용할 때 Variant Transforms 도구는 ALLELE_NUM 주석 필드의 값을 조회합니다. 이것은 VEP의 기능입니다.

  • --minimal_vep_alt_matching 플래그를 사용하여 Variant Transforms 도구로 --minimal 플래그를 시뮬레이션할 수 있습니다. 하지만 이로 인해 모호한 일치가 발생할 수 있습니다. 이러한 모호한 일치는 Dataflow Console에서 계산되고 보고됩니다. 생성된 테이블의 새 필드 아래에도 추가됩니다.

    최상의 결과를 얻으려면 --minimal_vep_alt_matching 대신 --minimal을 사용하세요.

VEP를 사용하여 VCF 파일에 자동으로 주석 추가

일반적으로 주석 필드를 파싱하여 개별 필드로 분할하려면 먼저 VCF 파일에 주석을 미리 추가해야 합니다. 하지만 Variant Transforms 도구를 사용할 때는 VCF 파일에 주석을 미리 추가할 필요가 없습니다.

Variant Transforms는 VCF 파일이 BigQuery에 로드되는 동안 자동으로 주석을 추가할 수 있습니다. 이렇게 하려면 Variant Transforms 도구로 VEP를 실행해야 하며, 그러기 위해서는 미리 구성된 Docker 이미지가 포함된 Compute Engine 가상 머신 인스턴스를 만들고 실행해야 합니다.

Variant Transforms 도구로 VEP 실행

Variant Transforms 도구를 VEP와 함께 실행하고 VCF 파일이 BigQuery에 로드되는 동안 자동으로 주석을 추가하려면 적어도 다음 플래그를 사용하여 Variant Transforms 도구를 실행해야 합니다.

  • --run_annotation_pipeline
  • --max_num_workers MAX_NUM_WORKERS 또는 --num_workers NUM_WORKERS
  • --annotation_output_dir gs://CLOUD_STORAGE_BUCKET

    CLOUD_STORAGE_BUCKET 변수에 대해 선택하는 버킷은 프로젝트에 속한 것이어야 합니다.

VEP를 실행할 때 Variant Transforms 도구에 전달할 수 있는 플래그의 전체 목록이 아래에 나와 있습니다.

플래그 기본값 설명
--run_annotation_pipeline False 자동 주석을 사용합니다.
--max_num_workers 또는 --num_workers 자동으로 결정됨 전자는 Cloud Dataflow가 자동 확장할 작업의 최대 작업자 수, 후자는 작업에 처음 할당된 작업자 수입니다.
--annotation_output_dir 없음 VEP 출력 파일을 위한 Cloud Storage 경로입니다. 입력 파일의 계층 구조가 이 위치에 복제되지만 새 파일이 _vep_output.vcf와 함께 추가됩니다. 디렉터리가 이미 존재하는 경우 Variant Transforms 도구가 실패합니다.
--vep_image_uri gcr.io/gcp-variant-annotation/vep_91 Cloud Genomics VEP Docker 이미지입니다. Google은 기본 이미지를 유지합니다(VEP 버전 91).
--vep_cache_path gs://gcp-variant-annotation-vep-cache/vep_cache_homo_sapiens_GRCh38_91.tar.gz(GRCh38 참조 염기서열에 맞춘 인간 게놈에 권장됨) 압축된 버전의 VEP 캐시를 위한 Cloud Storage 경로입니다. 이 파일은 --build_vep_cache.sh 스크립트를 사용하여 만들 수 있습니다.
--vep_info_field CSQ_VT 새 주석을 포함하는 새 INFO 필드의 이름입니다.
--vep_num_fork 2 단일 파일에서 VEP를 실행할 때 사용되는 로컬 프로세스 수입니다. 자세한 내용은 VEP 문서의 --fork [num_forks] 플래그를 참조하세요.
--shard_variants VEP를 사용하여 VCF 파일에 주석을 추가하기 전에 입력 파일을 더 작은 임시 VCF 파일로 분할할지 여부를 결정합니다. 입력 파일이 작은 경우(예: 각 VCF 파일에 포함된 변이가 50,000개 미만인 경우) 이 플래그를 True로 설정하면 컴퓨팅 리소스가 낭비될 수 있습니다.
--number_of_variants_per_shard 2,000 shard_variantsTrue이면 각 샤드에 기록된 최대 변이 수입니다. 샘플이 많은 데이터세트가 있는 경우 더 작은 값을 사용하는 것이 좋습니다. 더 작은 값을 제공하면 파이프라인이 완료되는 데 시간이 더 오래 걸릴 수 있습니다. 대부분의 경우 기본값이 권장됩니다.

문제해결

Variant Transforms 도구로 VEP를 실행하는데 도구가 실패하면 프로세스가 실행된 가상 머신의 로그를 제공하는 logs 디렉터리가 생성됩니다. logs 디렉터리는 --annotation_output_dir 플래그에 사용되는 경로에 지정됩니다. 표준 Variant Transforms 도구 로그 외에 이 로그 파일을 검사하여 문제의 원인을 파악할 수 있습니다.

VEP 플래그 --check_ref는 Variant Transforms 도구로 VEP를 실행할 때 자동으로 사용 설정됩니다. 따라서 logs 디렉터리의 파일에는 참조와 일치하지 않는 변이 목록이 포함됩니다. 올바른 VEP 캐시를 사용하는 경우에는 불일치가 자동으로 삭제되고 BigQuery 출력 표에 나타나지 않기 때문에 불일치가 나타나서는 안 됩니다.

샘플 쿼리

중요한 변이 추출

다음 샘플은 헤모글로빈의 베타 체인을 인코딩하는 HBB 유전자에서 중요한 변이를 추출하는 방법을 보여줍니다.

#standardSQL
SELECT reference_name, start_position, reference_bases, ALT.alt, CSQ.*
FROM vcf_imports_external.gnomad_genomes_chr_hg19 AS T, T.alternate_bases AS ALT, ALT.CSQ AS CSQ
WHERE CSQ.SYMBOL = "HBB" AND CSQ.IMPACT = "HIGH"
ORDER BY start_position

쿼리를 실행하면 다음이 반환됩니다.

reference_name start_position reference_bases alt allele 결과 IMPACT SYMBOL Gene Feature_type 기능 BIOTYPE EXON INTRON
1 11 5246947 G GC C frameshift_variant 높음 HBB ENSG00000244734 텍스트 변환 ENST00000335295 protein_coding 3/3
2 11 5246957 T C C splice_acceptor_variant 높음 HBB ENSG00000244734 텍스트 변환 ENST00000335295 proten_coding 2/2
3 11 5246957 T C C splice_acceptor_variant&non_coding_transcript_variant 높음 HBB ENSG00000244734 텍스트 변환 ENST00000475226 retained_intron 1/1
4 11 5247805 C T T splice_donor_variant 높음 HBB ENSG00000244734 텍스트 변환 ENST00000335295 protein_coding 2/2
5 11 5247805 C T T splice_donor_variant&non_coding_transcript_variant 높음 HBB ENSG00000244734 텍스트 변환 ENST00000475226 retained_intron 1/1
6 11 5247991 CAAAG C - frameshift_variant 높음 HBB ENSG00000244734 텍스트 변환 ENST00000335295 protein_coding 2/3
7 11 5247991 CAAAG C - frameshift_variant 높음 HBB ENSG00000244734 텍스트 변환 ENST00000380315 protein_coding 4/4
8 11 5248003 G A A stop_gained 높음 HBB ENSG00000244734 텍스트 변환 ENST00000335295 protein_coding 2/3
9 11 5248003 G A A stop_gained 높음 HBB ENSG00000244734 텍스트 변환 ENST00000380315 protein_coding 4/4
10 11 5248199 T A A stop_gained 높음 HBB ENSG00000244734 텍스트 변환 ENST00000335295 protein_coding 1/3
11 11 5248199 T A A stop_gained 높음 HBB ENSG00000244734 텍스트 변환 ENST00000380315 protein_coding 3/4

gnomAD 데이터세트에는 호출이나 샘플이 없지만 이 쿼리를 확장하여 중요한 변이가 있는 샘플을 선택할 수 있습니다.

데이터를 공개 BigQuery 데이터베이스와 결합

다음 샘플은 경로 데이터베이스 등 BigQuery에서 사용할 수 있는 다른 데이터베이스와 데이터를 결합하는 방법을 보여줍니다. ISB-CGC를 통해 사용할 수 있는 데이터베이스는 다음과 같습니다.

다음 샘플은 GO(Gene Ontology) 생물학적 프로세스(GO : 0006302)에 정의된 대로 이중 가닥 파열 복원 경로에서 중요한 변이를 추출하는 방법을 보여줍니다.

#standardSQL
SELECT
  reference_name, start_position, reference_bases, ALT.alt,
  CSQ.Consequence, CSQ.Impact, CSQ.SYMBOL
FROM
  `vcf_imports_external.gnomad_genomes_chr_hg19` AS T,
  T.alternate_bases AS ALT, ALT.CSQ AS CSQ
WHERE
  # Note: Matching based on symbol is "best effort" as the names may not be
  # standardized across sources.
  CSQ.SYMBOL IN (SELECT DB_Object_Symbol
                 FROM `isb-cgc.genome_reference.GO_Annotations`
                 WHERE GO_ID = 'GO:0006302')
  AND CSQ.IMPACT = "HIGH"
ORDER BY
  start_position

쿼리를 실행하면 다음이 반환됩니다.

reference_name start_position reference_bases alt 결과 IMPACT SYMBOL
1 5 917449 G T splice_acceptor_variant&non_coding_transcript_variant 높음 TRIP13
2 12 1022568 A C stop_gained 높음 RAD52
3 12 1022568 A C stop_gained 높음 RAD52
4 12 1022568 A C stop_gained 높음 RAD52
5 12 1023125 G A stop_gained 높음 RAD52
6 12 1023125 G A stop_gained 높음 RAD52
7 12 1023125 G A stop_gained 높음 RAD52
8 12 1023167 G A stop_gained 높음 RAD52
9 12 1023167 G A stop_gained 높음 RAD52
10 12 1023167 G A stop_gained 높음 RAD52
11 12 1023217 G T stop_gained 높음 RAD52
12 12 1023217 G T stop_gained 높음 RAD52
13 12 1023217 G T stop_gained 높음 RAD52
14 12 1025654 TGA T frameshift_variant 높음 RAD52
15 12 1025699 T G splice_acceptor_variant&non_coding_transcript_variant 높음 RAD52
16 12 1036004 GC G frameshift_variant 높음 RAD52
17 12 1038977 C CT frameshift_variant 높음 RAD52