三級分析 / VCF Format === ###### tags: `基因體/三級分析` ###### tags: `生物資訊`, `基因體`, `三級分析`, `VCF` <br> [TOC] <br> ## 基本說明 - **全名**： - Variant Call Format (.vcf 檔) - **用途**： - 用來描述「**基因序列變異**」的文字檔 - 上一步由 HaplotypeCaller 等相關工具找出序列變異 - 接著，將變異結果寫到 .vcf 檔（用來儲存變異的資料） - 下一步供「基因標注/註釋」使用，並確認基因相關功能、是否有相對應的疾病 - **維護者**： - GA4GH - The Global Alliance for Genomics and Health - 全球基因組學與健康聯盟 - **目前版本**： - 4.X - **規格 (Specifications, Spec.)**：(非專業人員勿入) - [SAM/BAM and related specifications](http://samtools.github.io/hts-specs/) ([github](https://github.com/samtools/hts-specs)) - The Variant Call Format (VCF) Version 4.X Specification - [VCF v4.1](http://samtools.github.io/hts-specs/VCFv4.1.pdf) - [VCF v4.2 - 2015/10](http://samtools.github.io/hts-specs/VCFv4.2.pdf) - [VCF v4.3 - 2022/08/22](http://samtools.github.io/hts-specs/VCFv4.3.pdf) - [VCF v4.4 - 2023/01/27](http://samtools.github.io/hts-specs/VCFv4.4.pdf) <br> ## 預備知識 - ### 術語 - **chrom** ← 可以想成用 Chrome 瀏覽器，瀏覽基因 chromosome 染色體 [ˋkromə͵som] - **ref** reference base(s) 「參考基因體」的「鹼基序列」 [ˋrɛfərəns] - **alt** alternate base(s) 「待測樣本」的「鹼基序列」 [ˋɔltɚnɪt] - **contig** - 連續片段 (from ChatGPT) - 片段重疊群 - 將片段序列組裝好的連續體，通常指單一染色體的「鹼基序列」 ![](https://i.imgur.com/uKC2jSA.png) ([圖片來源](![](https://i.imgur.com/ze5iJYd.png) )) - **allele** [əˋlil] 孟德爾遺傳學中的）對偶基因 對偶基因(早期) / 等位基因(晚期) - 何謂「等位基因」？ - [遺傳 02 -- 基因、等位基因與棋盤方格法 (影片)](https://www.youtube.com/watch?v=zRTRg2AkS9g&t=22s) - 在單套的染色體上稱為等位基因 - 在單套染色體上的片段，還能稱為基因嗎？ - 為避免混淆，在單套染色體上的片段，稱為等位基因 - 為什麼要叫它們等位基因？ - 因為他們位置相等 ![](https://i.imgur.com/qTao1k3.png) - 等位基因位於相對位置 - 控制同一性狀的等位基因，會位於同源染色體上相對位置 - [國一下生物2-2(1)基因與遺傳：等位基因 (影片)](https://www.youtube.com/watch?v=-n1nCwf2Mqc) <br> ## VCF 實際案例詳解 - [VCF - D15780_S13_L001.vcf 樣本](https://hackmd.io/JhM5KGz7Q9yIzZ6_eXQxDA) - [VCF - Genomics1152samples_c2_1000.vcf (豬隻小資料)](https://hackmd.io/uYwoH4XfTUuFLdrluGdWkQ) <br> ## VCF 主體結構 ![](https://us.v-cdn.net/5019796/uploads/FileUpload/70/9eec0b6faaa664f7630abddaf15594.png) - ### 圖解： [![](https://i.imgur.com/M2O1oGD.png)](https://i.imgur.com/M2O1oGD.png) (圖片來源：[我們的基因體時代 OUR "GENE"RATION](https://weitinglin.com/2017/05/29/vcfvariant-call-format-基因突變資料儲存格式/)) <br> ## VCF Header ### 快速理解 header - 類似 HTML 的 metadata 概念，用於存放網頁相關的屬性，像是編碼方式、標題、作者、顯示在搜尋引擎上的描述 ```html <meta http-equiv="content-type" content="text/html; charset=utf-8"> ``` ````html <meta charset="utf-8"> <meta name="author" content="tj_tsai"> <meta name="description" content="Introduction of VCF"> ``` - 其中 FORMAT 和 INFO 比較特別 - 意義：用來描述變異位點的標注資訊 - 語法：用來描述一個表格的欄位，很像資料庫的「綱要(schema)」 - [schema 是什麼樣的概念？](https://www.sqlite.org/cli.html) ```sql CREATE TABLE lost_and_found( rootpgno INTEGER, -- root page of tree pgno is a part of pgno INTEGER, -- page number row was found on nfield INTEGER, -- number of fields in row id INTEGER, -- value of rowid field, or NULL c0, c1, c2, c3... -- columns for fields of row ); ``` - 可以把它看成程式語言的「變數宣告」 <br> ### Header Spec - 格式 ```##屬性名稱="自定義值"``` - VCF屬性名稱，常見範例： - `##fileformat=VCFv4.2` (VCF版本) - `##FILTER=<...>` (是否通過所設定的篩選條件) - `##INFO=<...>` (延伸資訊，可註釋更多資訊) - `##FORMAT=<...>` (sample 1~N 的格式資訊) - `##contig=<...>` (定義序列資訊) - 自定義值 - 單一值 ```vcf ##VCF屬性名稱=value ``` - 範例 ``` ##fileformat=VCFv4.2 ##bcftools_normVersion=1.7+htslib-1.9 ##bcftools_normCommand=norm -m- -o aed5ad61ea26e520e12aaf647567d31b06296fdb.norm.vcf /workspace/datasets/germline/HG002/wes-output/output.vcf; Date=Wed Dec 28 11:19:03 2022 ``` - `<...>`: field 清單 (欄位清單) (子屬性表格) > - i.e. key/value(子屬性名稱/子屬性值)清單 > - column vs field > - column 是指一個表格的縱向儲存格之集合; > - field 是指一個結構訊息中的一個屬性 (像是 C 語言的 struct); > - [column 只是顯示一個 field 內容的一種方式](https://learn.microsoft.com/en-us/office/vba/outlook/concepts/forms/difference-between-a-field-and-a-column) > - https://chat.openai.com/chat > ![](https://i.imgur.com/4qfxaFc.png) > > Q: field 和 column 有何區別？ > > A: "Field" 和 "column" 有相似的意思，都指的是資料表裡的一個欄位。通常來說，"field" 和 "column" 可以互換使用。 > > > > 但是，在一些用法上可能有所不同。例如，"field" 可能會指更廣義的欄位概念，不僅僅限於資料表裡的欄位。在某些情況下，"field" 可能指的是表單裡的一個輸入框，或是網頁上的一個輸入欄位。 > > > > 所以，"field" 和 "column" 在某些情況下可能有所不同，但在資料庫領域裡，通常可以把它們當作同義詞使用。 ```vcf ##VCF屬性名稱=<key1=value1,key2=value2,key3=value3,...> ``` - 實際範例 ``` ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> ``` - `ID=子屬性名稱` - `Description="描述子屬性資訊"` - **`Number=[0,1,2,3,...|A|R|G|.]` (子屬性的資料個數)** - `0`: 子屬性值有 0 個 (沒有值) (Type 必然是 Flag) - `1`: 子屬性值有 1 個 (一個, single) - `2`: 子屬性值有 2 個 (一對, pair, 通常用 comma 隔開) - `3`: 子屬性值有 3 個 (多個, tuple, 通常用 comma 隔開) - ... - **`Number=A`: ALT 欄位的額外值** - 對 ALT 欄位進行擴充/註釋 - ALT 有 3 個，對應的子屬性值就有 3 個值 - 例如 - REF:A -> ALT:[T|C|G] - REF:T -> ALT:[A|C|G] - REF:C -> ALT:G - 實際範例 ``` ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed"> ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed"> #CHROM POS REF ALT INFO chr1 14590 G A AC=1;AF=0.5;... chr1 14599 T A AC=1;AF=0.5;... ``` - 資料來源：HG002 + dbsnp - 只有一個sample, 就只有一個ALT <br> - **`Number=R`: REF 欄位的額外值 + ALT 欄位的額外值** - 對 REF 欄位進行擴充/註釋 - ALT 有 2 個，對應的子屬性值就有 3(1+2) 個值 - 實際範例：AD欄位 ``` ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed"> #CHROM POS REF ALT FORMAT sample chr1 14590 G A GT:AD:DP:GQ:PL 0/1:1,1:2:39:39,0,39 chr1 14599 T A GT:AD:DP:GQ:PL 0/1:1,1:2:39:39,0,39 ... chr1 16949 A C GT:AD:DP:GQ:PL 0/1:16,24:40:99:502,0,378 ``` - 資料來源：HG002 + dbsnp - **資料解讀：REF=`G`, ALT=`A`, AD=`1,1`, DP=`2`** - 表示獲得 1 個 G 的讀數(read) 和 1 個 A 的讀數(read) - 定序深度DP=1+1=2 - **資料解讀：REF=`A`, ALT=`C`, AD=`16,24`, DP=`40`** - 表示獲得 16 個 A 的讀數(read) 和 24 個 A 的讀數(read) - 定序深度DP=16+24=40 - 參考資料 - [[Biostars] What Is Ad (Allelic Depth) In 1000Genomes Vcf?](https://www.biostars.org/p/5268/) Allele specific depth. i.e. if ref is 'A' and alt is 'G' and AD is '6,9' you got 6 A reads and 9 G reads. - [[GATK] Allele Depth (AD) is lower than expected](https://gatk.broadinstitute.org/hc/en-us/articles/360035532252-Allele-Depth-AD-is-lower-than-expected) <br> - **`Number=G`: GenoType 欄位的額外值** > see Genotype Ordering - 對 GenoType 欄位進行擴充/註釋 - 實際範例：PL欄位 ``` ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification"> #CHROM POS REF ALT FORMAT sample chr1 14590 G A GT:AD:DP:GQ:PL 0/1:1,1:2:39:39,0,39 chr1 14599 T A GT:AD:DP:GQ:PL 0/1:1,1:2:39:39,0,39 ... chr1 16949 A C GT:AD:DP:GQ:PL 0/1:16,24:40:99:502,0,378 ``` - **資料解讀：GT=`0/1`, PL=`39,0,39`** - P=2 (ploidy, 倍性, 人類為2倍體，來自父母各一條) - 基因型：0/0(同合子), 0/1(雜合子), 1/1(同合子) - 基因似然率 - 0/0: 10^(-39/10) = 0.000126 - 0/1: 10^(-0/10) = 1 => ==GT== - 1/1: 10^(-39/10) = 0.000126 - **資料解讀：GT=`0/1`, PL=`502,0,378`** - 基因似然率 - 0/0: 10^(-502/10) ~ 0 - 0/1: 10^(-0/10) = 1 => ==GT== - 1/1: 10^(-378/10) ~ 0 - 參考資料 - [[biostars] Genotype Likelihoods](https://www.biostars.org/p/76462/#76508) #phred-scaled likelihood scores <br> - **`Number=.`: 子屬性值之個數未知、不定個數** - 實際範例 ``` ##INFO=<ID=dbsnp_CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter (from /workspace/datasets/variants/ncbi/dbsnp/00-All.vcf.gz)"> ##INFO=<ID=dbsnp_TOPMED,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on TOPMed, starting with the reference allele followed by alternate alleles as ordered in the ALT column. The TOPMed minor allele is the second largest value in the list. (from /workspace/datasets/variants/ncbi/dbsnp/00-All.vcf.gz)"> #CHROM POS REF ALT INFO chr1 14590 G A dbsnp_RS=707679;dbsnp_VC=SNV;dbsnp_TOPMED=0.93632931957186544,0.06367068042813455 chr1 14599 T A dbsnp_RS=707680;dbsnp_VC=SNV;dbsnp_CAF=0.8524,0.1476,.;dbsnp_TOPMED=0.90694285168195718,0.09304918450560652,0.00000796381243628 chr1 14604 A G dbsnp_RS=541940975;dbsnp_VC=SNV;dbsnp_CAF=0.8524,.,0.1476;dbsnp_TOPMED=0.88364870030581039,0.00000796381243628,0.11634333588175331 ``` - [rs707679](https://www.ncbi.nlm.nih.gov/snp/rs707679) - dbsnp_TOPMED=0.93632931957186544,0.06367068042813455 - ref: 0.9363 - G>A: 0.0637 - [rs707680](https://www.ncbi.nlm.nih.gov/snp/rs707680) - dbsnp_CAF=0.8524,0.1476,. - ref: 0.8524 - T>A: 0.1476 - T>C: . - [rs541940975](https://www.ncbi.nlm.nih.gov/snp/rs541940975) - dbsnp_CAF=0.8524,.,0.1476 - ref: 0.8524 - A>C: . - A>G: 0.1476 <br> - **Type=`[Integer|Float|Flag|Character|String]`** - Flag 實際範例 ``` ##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP Membership"> #CHROM POS REF ALT INFO chr1 14590 G A AC=1;AF=0.5;... chr1 14599 T A AC=1;AF=0.5;DB;... ``` - 只要填 DB 就表示啟用/參與，沒有填 DB 就表示關閉/無參與 - 就像 html 的 tag 屬性: disabled ```<input type="text" name="lname" disabled>``` - **[##ALT=<ID=type,Description=description>](https://www.internationalgenome.org/wiki/Analysis/Variant%20Call%20Format/VCF%20(Variant%20Call%20Format)%20version%204.0/encoding-structural-variants)** - ID 是用來指出結構變異的類型 - ID 可以允許多個，用冒號隔開 - ID 第一個值(主要類別)，一定是下列其中之一： - DEL - Deletion relative to the reference - 相對於參考基因體的刪除 - INS - Insertion of novel sequence relative to the reference - 相對於參考基因體的新序列插入 - DUP (duplication, 重複) - Region of elevated copy number relative to the reference - 相對於參考基因體的提昇複製數區域 - INV - Inversion of reference sequence - 參考基因體序列的反轉 - CNV - Copy number variable region (may be both deletion and duplication) - 拷貝數變異區（可能同時是插入與重複） - **contig (片段重疊群) 類別之屬性：** - 用來紀錄各個染色體，及其他序列的長度 - 人類 - 第1～22條是體染色體 - 第23條(X, Y)是性染色體 ``` ##contig=<ID=X,length=155270560> ##contig=<ID=Y,length=59373566> ``` - MT(mitochondrion) 粒線體染色體 ``` ##contig=<ID=MT,length=16569> ``` - GLXXXXXX 無法定位的基因體片段序列 [[NCBI] GL000207.1](https://www.ncbi.nlm.nih.gov/nuccore/GL000207.1?report=genbank) ``` ##contig=<ID=GL000207.1,length=4262> ``` [[NCBI] GL000226.1](https://www.ncbi.nlm.nih.gov/nuccore/224183275/) ``` ##contig=<ID=GL000226.1,length=15008> ``` - [NC_007605](https://www.ncbi.nlm.nih.gov/nuccore/NC_007605) (NC:Nucleotide, 核苷酸) 人類疱疹病毒第四型 ``` ##contig=<ID=NC_007605,length=171823> ``` - **unlocalized? unplaced?** - unlocalized 未定位的 - 已知染色體編號，但不知所在位置(locus) - 在 contig 中，用 GLXXXXXX 表示 - unplaced 未放置的 - 不知道要歸類到哪條染色體 - 在 contig 中，用 JHXXXXXX 表示(???) - https://www.biostars.org/p/188628/ > chr_random - are called "unlocalized" sequences. The chromosome is known, but not the location on the chromosome. chrUn_ - are called "unplaced" sequences. > > chr_random - 稱為 "未定位的" 序列。 染色體是已知的，但不是指染色體上的位置（該序列在染色體上的位置是未知的）。chr_Un - 稱為 "未放置的"。他們可能屬於定序基因體，但此時位置尚不清楚。 - [人類染色體版本釋疑 - 序列命名](https://yourgene.pixnet.net/blog/post/117778207) - GRCh37在主要的染色體上的命名方式為 1, 2, 3, 4, ….X, Y, MT。 - hg19的命名則為 chr1, chr2, chr3, chr4, …..chrX, chrY, chrM。 - [[Note] 人類基因體計畫(The Human Genome Project, HGP)](https://hackmd.io/qmzwm4idRyWROeeboqpfSw#%E9%81%8E%E5%8E%BB%E8%A8%88%E7%95%AB%EF%BC%9A) - #### INFO 類別之屬性： - 意義：用來描述變異位點的標注資訊 - 可看成是： - 參與決策的因素 f(x1, x2, x3, ...) = 有無變異？ - 可以想像成： - 投資股票用到很多常見的 - 基本資訊 - 本益比、營業毛利率、稅前淨利率、... - 基本交易指標 - 成交量、委賣張數、委買張數、融資融券比、... - 進階交易指標 - 移動平均線、KD、MACD、RSI、乖離率、... - 甚至很多自創的交易指標 - 夏林指數、布林通道、... - 凡人不懂沒關係，那是專業投資者的心法 - 常見的 INFO 屬性/欄位： ==(不負責任解釋)== - AC (Allele Count) > **Allele count in genotypes, for each ALT allele, in the same order as listed** > 每個 ALT 等位基因的基因型等位基因計數，順序與所列順序相同 > 例如：AC=1, AC=2 > 在人類 (見[染色體倍性](https://zh.m.wikipedia.org/zh-tw/%E6%9F%93%E8%89%B2%E9%AB%94%E5%80%8D%E6%80%A7))， > - 精子和卵子是單倍體 (AC=1) > - 體細胞是雙倍體 (AC=1 表示一條出現變異, AC=2 表示兩條都出現變異) > - XX: > - X: AC=1 or AC=2 > - XY: > - X: AC=1 > - Y: AC=1 - AF (Allele Frequency) (等位基因頻率) (個體變異頻率) (群體變異頻率) > Allele Frequency, for each ALT allele, in the same order as listed > 每個 ALT 等位基因的等位基因頻率，與所列順序相同 > 例如：AC=1;AF=0.500;, AC=2;AF=1 > - 單倍體 > - AC=1;AF=1; > - 雙倍體 > - AC=1;AF=0.500; > - AC=2;AF=1; > - 三倍體 > - AC=1;AF=0.333; > - AC=2;AF=0.666; > - AC=3;AF=1; - AN > Total number of alleles in called genotypes - DB > dbSNP Membership - DP > Total Depth; > Approximate read depth; > some reads may have been filtered; - DS > Were any of the samples downsampled? - MQ > RMS Mapping Quality - QD > Variant Confidence/Quality by Depth - 不常見的 INFO 屬性/欄位： - BaseQRankSum > Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities - Dels > Fraction of Reads Containing Spanning Deletions - ExcessHet > Phred-scaled p-value for exact test of excess heterozygosity - FS > Phred-scaled p-value using Fisher's exact test to detect strand bias - HaplotypeScore > Consistency of the site with two (and only two) segregating haplotypes - HRun > Largest Contiguous Homopolymer Run of Variant Allele In Either Direction - InbreedingCoeff > Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation - MLEAC > Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed - MLEAF > Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT alle in the same order as listed - MQ0 > Total Mapping Quality Zero Reads - MQRankSum > Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities - ReadPosRankSum > Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias - RPA > Number of times tandem repeat unit is repeated, for each allele (including reference) - RU > Tandem repeat unit (bases) - SB > Strand Bias - SOR > Symmetric Odds Ratio of 2x2 contingency table to detect strand bias - STR > Variant is a short tandem repeat - VQSLOD > log10-scaled probability of variant being true under the trained gaussian mixture model - 其他屬性/欄位： - [VCF详解](https://www.cnblogs.com/xiaofeiIDO/p/6610157.html) - [转载-VCF格式详解](http://www.bio-info-trainee.com/863.html) (有中文解釋 INFO 欄位) - [变异信息那些事](https://www.jieandze1314.com/post/cnposts/variant/) (VCF超詳解) - #### Format 類別之屬性： - 意義：用來描述變異位點的標注資訊 <br> - AD (Allelic depths, 定序深度) > Allelic depths for the ref and alt alleles in the order listed > 在所列的順序中，ref 和 alt 等位基因的等位深度 - DP (Read Depth, reads數量) > Read Depth (only filtered reads used for calling) - GQ (Genotype Quality, 基因型品質) > Genotype Quality - GT (Genotype, 基因型) > Genotype - PL (Phred-scaled likelihoods, Phred-scaled 似然值) > Normalized, Phred-scaled likelihoods for AA,AB,BB genotypes where A=ref and B=alt; not applicable if site is not biallelic > AA，AB，BB 基因型的標準化 Phred-scaled 似然值，其中 A = ref 和 B = alt (如果該位置不是雙等位的，則不適用 PL) - #### 延伸屬性 - [VAF:Variant Allel Frequency簡介](https://ppfocus.com/0/sc799a472.html) $$VAF = \frac{allel\_depth}{tatal\_depth} = \frac{AD}{DP}$$ - **VCF**: Variant Allel Frequency 變異等位基因頻率 - **MAF**: Minor Allel Frequency 次等位基因在人羣中的頻率 - **VAF的值的大小有什麼含義呢？** 以二倍體生物爲例，假設所有的細胞中該位點都是雜合的，那麼50%的染色體上包含了ref allel, 另外50%的染色體上包含了alt allel, 則測序結果中該位點的VAF值應該爲0.5。對於germline genotype而言，一個可靠的突變位點其VAF的值應該在0.5附近。 - 對於生殖變異的檢測，認爲其VAF的偏移來源於拷貝數的變化， - 對於體細胞檢測而言，更多的認爲VAF的偏移來源於腫瘤細胞的異質性。 <br> ## VCF Body - ### body 快速理解： - 如果你是一隻病毒，你要怎麼把你身上的病毒，植入人類身體？ - 尋找第2號染色體 -> #CHROM - 尋找能匹配的起始碼 -> POS - 開始植入置換 -> src: REF, dest/target: ALT - ### body 區段講解： - 每一行代表一個變異(variant)資訊 - 每一行是由多個欄位組成，每一個欄位用 tab 區分 <br> - 底下有 5 種常見的變異結構，與對應的表示法： 起始索引為 1 (base = 1) ![](https://weitinglindotcom.files.wordpress.com/2017/05/screenshot39.png) <br> - **[欄位資訊](https://www.itread01.com/content/1548698402.html)**： > 參考資料 > - [dbSNP2.0 VCF File / Fixed Columns](https://www.ncbi.nlm.nih.gov/snp/docs/products/vcf/redesign/#loc_all_info) > - **CHROM** (Chromosome) 表示變異位點是在哪個 contig 裡 call 出來的 - **POS** (Position) 「變異位點」相對於「參考基因體」所在的位置 (起始索引值 = 1) - **ID** (Identifier) 變異體(variant)的ID ``` if (call 出來的 SNP) 位於 (dbSNP 資料庫): 顯示 dbSNP 裡相應的 rs 編號 else: 用 . 表示其為一個新的變異體(novel variant) ``` - [[Spec-v4.2] 1.4.1 Fixed fields](https://samtools.github.io/hts-specs/VCFv4.2.pdf) > 3. ID - identifier: Semicolon-separated list of unique identifiers where available. If this is a dbSNP variant it is encouraged to use the rs number(s). No identifier should be present in more than one data record. If there is no identifier available, then the missing value should be used. (String, no whitespace or semicolons permitted) - **REF** (Reference Base(s)) 「參考基因體」的「鹼基序列」 - **ALT** (Alternate Base(s)) 「待測樣本」的「鹼基序列」 (alternate 代替者，就是指實驗樣本、待測樣本) - [What is REF and ALT ? Does all variations are mutations ?](https://www.biostars.org/p/372865/) - Ref = the allele in the reference genome. - Alt = any other allele found at that locus. > People use **mutation** in different ways and there seems to be no real consensus. Some people call any variant with a frequency lower than a certain value a **mutation**, and those with a higher frequency **polymorphisms**. Personally, I stick to the word **variant** for everything, and only use **mutation** when referring to the act of mutation, such as when talking about a somatic mutation or a de novo mutation (although they are both still variants). > > 人們以不同的方式使用**突變**，似乎沒有真正的共識。有些人將頻率低於某個值的任何變異稱為**突變**，而將頻率較高的變異稱為**多態性**。就我個人而言，我對所有事物都堅持使用**變體**一詞，並且僅在提及突變行為時才使用**突變**，例如在談論體細胞突變或從頭突變時（儘管它們仍然是變體）。 - **QUAL** (Quality) > Quality, Phred格式(Phred_scaled)的品質分數，可以理解為所 call 出來的變異位點的品質分數。表 示在該位點存在variant的可能性；該值越高，則variant的可能性越大； > 計算方法： > > ① Q = -10 log10(P\) > Q 表示品質分數 > P 表示這個位點發生錯誤的概率 > > ② Phred 值 Q = -10 log10(1-p) > p 為 variant 存在的概率; > > 通過計算公式可以看出值為10的表示錯誤概率為0.1，該位點為variant的概率為90%。同理，當Q=20時，錯誤率就控制在了0.01。 - [[wiki] Phred quality score](https://en.wikipedia.org/wiki/Phred_quality_score) - **FILTER** (Filter Status) > 使用上一個 QUAL 值來進行過濾的話，是不夠的。 > > 理想情況下，QUAL這個值應該是用所有的錯誤模型算出來的，這個值就可以代表正確的變異位點了，但是事實是做不到的。因此，還需要對原始變異位點做進一步的過濾。無論你用什麼方法對變異位點進行過濾，過濾完了之後，在FILTER一欄都會留下過濾記錄， > > 如果是通過了過濾標準，那麼這些通過標準的好的變異位點的 FILTER 一欄就會註釋一個 PASS， > > 如果沒有通過過濾，就會在 FILTER 這一欄提示除了 PASS 的其他資訊。 > > 如果這一欄是一個 “.” 的話，就說明沒有進行過任何過濾。 - **INFO** (Additional Information) - [詳見 header 說明](#INFO-類別之屬性：) - **FORMAT** - [詳見 header 說明](#Format-類別之屬性：) - ***your_sample_name*** - 由原始的 bam 檔的 @RG 標籤提供，如： ``` @RG ID:D15780_S13_L001 SM:D15780_S13_L001 PL:Illumina ``` - @RG: [Read Groups](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) - `ID` = Read group identifier - `PU` = Platform Unit - `SM` = Sample - `PL` = Platform/technology used to produce the read - `LB` = DNA preparation library identifier - @RG / SM = Sample > The name of the sample sequenced in this read group. > 在該片段序列群組中定序的樣本名稱 > > GATK tools treat all read groups with the same SM value as containing sequencing data for the same sample, > GATK 工具將具有相同 SM 值的所有片段序列群組，視為包含相同樣本的定序資料 > (講單講，相同的 SM 值都是同一個樣本) > > and this is also the name that will be used for the sample column in the VCF file. > 這也是將用於 VCF 檔案中的樣本欄位的名稱 > > Therefore it is critical that the SM field be specified correctly. When sequencing pools of samples, use a pool name instead of an individual sample name. > 因此，正確指定 SM 欄位至關重要。對樣本池進行定序時，使用池的名稱而不是個別的樣本名稱。 - @PG [program group](https://gatk.broadinstitute.org/hc/en-us/articles/360036713171-AddOrReplaceReadGroups-Picard-) <br> ## VCF 相關工具 - ### VCF 可視化工具 - [Integrative Genomics Viewer (IGV)](https://software.broadinstitute.org/software/igv/viewing_vcf_files) > The Integrative Genomics Viewer (IGV) is a high-performance visualization tool for interactive exploration of large, integrated genomic datasets. It supports a wide variety of data types, including array-based and next-generation sequence data, and genomic annotations. > 整合基因體檢視器 (IGV) 是一種高效能視覺化工具，用於互動探索大型整合基因體資料集。它支援各種資料類型，包含基於陣列和次世代的序列資料，以及基因體標注。<br> ![](https://software.broadinstitute.org/software/igv/sites/cancerinformatics.org.igv/files/images/vcfwgenotypes.jpg) - 下載 - [Dowanlod / IGV for Linux](https://data.broadinstitute.org/igv/projects/downloads/2.8/IGV_Linux_2.8.0.zip) - 啟動 ```bash IGV_Linux_2.8.0$ bash igv_hidpi.sh ``` - 預設 port 為 60151 <br> - 其他工具: - on-going - ### VCF 比對工具 - [How to compare 2 VCF files](https://www.biostars.org/p/100901/) - http://snpeff.sourceforge.net/SnpSift.html#concordance ```bash $ java -Xmx1g -jar SnpSift.jar concordance -v genotype.vcf sequencing.vcf > concordance.txt ``` - http://vcftools.sourceforge.net/ (old) - https://vcftools.github.io/index.html (new) - on-going <br> ## [FORMAT 欄位](#Format-類別之屬性：) ### FORMAT 範例 1 ``` 1 899282 rs28548431 C T [CLIPPED] GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26 ``` - Chrom = 1, POS = 899282 - 變異點位於 1 號染色體的第 899282 個鹼基 - (起始索引值為1) - ID 欄位為 rs28548431, 非空資料 - 表示該變異的資訊是「 已知的 (已經確認的SNP) 」 - 存在於 dbSNP 資料庫中，rs 編號為 rs28548431 - REF: 參考基因體的鹼基為 C (當成群體) - ALT: 待測樣本的鹼基為 T (當成個體) <br> ``` C T ... GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26 ``` - **GT (Genotype, 基因型)** - 如果等位基因可以區分來自父源或母源 該基因型可以表示為 phased genetype (已定相基因型) <br> - [人類基因組的Phasing原理是什麼？](https://kknews.cc/zh-tw/science/9lzm4pl.html) ![](https://i.imgur.com/pYH5TOA.jpg) ![](https://i.imgur.com/njuIPI9.jpg) <br> - [目前高通量定序結果，大都是 unphased genotype](https://www.jianshu.com/p/5cf18078a5a7) <br> - **unphased genetype (混淆基因型, 未定相基因型)** - 對於二倍體來說：有兩個等位基因，表示為 -/- - [1/0 和 0/1 寫法是等價的](https://www.jianshu.com/p/5cf18078a5a7) - 對於三倍體來說：有三個等位基因，表示為 -/-/- - 對於四倍體來說：有四個等位基因，表示為 -/-/-/- - **phased genetype (已定相基因型)** - 對於二倍體來說：有兩個等位基因，表示為 -|- - [1|0 和 0|0 表示不同基因](https://www.jianshu.com/p/5cf18078a5a7) - 對於三倍體來說：有三個等位基因，表示為 -|-|- - 對於四倍體來說：有四個等位基因，表示為 -|-|-|- <br> ``` C T ... GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26 ``` - 已知群體是 C，個體是 T，根據定序資料，推測 GT 為 0/1 - phased genetype 有 4 種：0|0, 0|1, 1|0, 1|1 - 但次世代定序結果為 unphased genetype，結果只有 3 種： - 0/0 (C/C) - 0/1 (C/T) 和 1/0 (T/C) 算一種 ... (無法區分來源) - 1/1 (T/T) - 其中 0/0 (C/C), 1/1 (T/T) - 等位基因都一樣 - 同型合子(homozygous) / 純合子 - 其中 0/1 (C/T) 和 1/0 (T/C) - 等位基因不一樣 - 異型合子(heterozygous) / 雜合子 - 0/1 (C/T) 和 1/0 (T/C) 信心度比較高，所以 GT 為 0/1 - 兩個 ALT 的例子：https://www.biostars.org/p/96580/ - 0/2的例子：https://www.internationalgenome.org/wiki/Analysis/vcf4.0/ <br> ``` C T ... GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26 ``` - AD (Allelic depths, 定序深度) - Allelic depths for the ref and alt alleles in the order listed 在所列的順序中，REF 和 ALT 等位基因的等位深度 - AD = 1,3 - 1 表示涵蓋到 REF 鹼基(=C)的片段序列(read)有 1 個 - 3 表示涵蓋到 ALT 鹼基(=T)的片段序列(read)有 3 個 - DP (Read Depth, reads 數量) - DP = 4 - sum(AD) = sum([1, 3]) = 4 <br> ``` C T ... GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26 ``` - PL (Phred-scaled likelihoods, Phred-scaled 似然值) > Normalized, Phred-scaled likelihoods for AA,AB,BB genotypes where A=ref and B=alt; not applicable if site is not biallelic > AA，AB，BB 基因型的標準化 Phred-scaled 似然值，其中 A = ref 和 B = alt (如果該位置不是雙等位的，則不適用 PL) - PL=103,0,26，分別對應於 - 0/0 (C/C) 的 PL 值 = 103 - 0/1 (C/T) 和 1/0 (T/C) = 0 - 1/1 (T/T) = 26 - [Math notes: How PL is calculated in HaplotypeCaller](https://gatkforums.broadinstitute.org/gatk/discussion/5913/math-notes-how-pl-is-calculated-in-haplotypecaller) ``` # Alleles Reference: A Read: T # Conditional probabilities calculated by HC P(AA | Data) = 0.000001 P(AT | Data) = 0.000100 P(TT | Data) = 0.010000 ``` - PL 為延伸指標（非原本的定序資料之統計值） - 由 GATK 的變異偵測工具所計算，如 HaplotypeCaller - 其中 - AA: A/A=0/0 - AT: A/T=0/1, T/A=1/0 ... (無法區分來源) - TT: T/T=1/1 - 在實際的資料中，這三種基因型的機率總和應為 1 <br> - 想把小數點表示，轉成 **「比較可讀」** 的形式 - 把它映射到 0 ~ 100 - 所以要取 log - AA: log(0.000001) = log(10e-6) = -6 - AT: log(0.000100) = log(10e-4) = -4 - TT: log(0.010000) = log(10e-2) = -2 - 再乘上 -10 - AA: -6 x -10 = 60 - AT: -4 x -10 = 40 - TT: -2 x -10 = 20 - 原始 PL 公式：PL = -10 log(P\) - 減去最小的 PL 值 = min([60, 40, 20]) - AA: 60 - 20 = 40 - AT: 40 - 20 = 20 - TT: 20 - 20 = 0 - PL 值越低，表示「是該基因型」的機率越高 - 標準化 PL 公式： - PL = -10 log(P\) - min(PL list) - [Phred-scaled Quality Scores](https://software.broadinstitute.org/gatk/documentation/article.php?id=4260) <br> ``` C T ... GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26 ``` - PL=103,0,26，分別對應於 - 0/0 (C/C) 的 PL 值 = 103 - 0/1 (C/T) 和 1/0 (T/C) = 0 - 1/1 (T/T) = 26 - PL 對應的可能性機率： ```math PL = -10 log(P) -PL/10 = log(P) 10^(-PL/10) = 10^log(P) 10^(-PL/10) = P P = 10^(-PL/10) ``` | 基因型分配 | PL | P | | -------- | -------- | -------- | | 0/0 (C/C) | 103 | ```10^(-103/10) = 5.01e-11 ≒ 0%``` | | 0/1 (C/T)<br>1/0 (T/C) | 0 | ```10^(-0/10) = 1 = 100%``` | | 1/1 (T/T) | 26 | ```10^(-26/10) = 0.002512 = 0.25%``` | - 可以很肯定，該位點有變異，因為 0/0 趨近於 0 - 但基因型分配，可能不是正確的 - 可能是雜合 het (heterozygous, 雜合子, 異型合子) - 也可能是 hom-var (homozygous with the variant allele, 帶有變異等位基因的純合子) <br> ``` C T ... GT:AD:DP:GQ:PL 0/1:1,3:4:26:103,0,26 ``` - GQ (Genotype Quality, 基因型品質) - PL=103,0,26，分別對應於 - 0/0 (C/C) 的 PL 值 = 103 - 0/1 (C/T) 和 1/0 (T/C) = 0 - 1/1 (T/T) = 26 - 基因型的最可能分配(assignment)形式為 - PL 值最低的那一型（應該就是多數決） - 0/1 (C/T) 和 1/0 (T/C) - 所以，GT = 0/1 - GQ 定義為 min(第二低的 PL 值, 上限值 99) - GQ = 26 - 為何要取第二低的 PL 值？ https://software.broadinstitute.org/gatk/documentation/article.php?id=1268 > The values of the PLs are normalized so that the most likely PL is always 0, so the GQ ends up being equal to the second smallest PL, unless that PL is greater than 99. > > PL 的值是被正規化的，使得最可能的 PL 總是為 0。因此 GQ 最後定為第二個最小的 PL，除非 PL 大於 99。(亦即，若第二個最小的 PL 若大於 99，則 PL 為 0) <br> ### [[VCFv4.3][Spec] 1.6.2 Genotype fields](https://samtools.github.io/hts-specs/VCFv4.3.pdf) | Field | Number | Type | Description | |-------|--------|------|-------------| | AD | R | Integer | Read depth for each allele | | ADF | R | Integer | Read depth for each allele on the forward strand | | ADR | R | Integer | Read depth for each allele on the reverse strand | | DP | 1 | Integer | Read depth | | EC | A | Integer | Expected alternate allele counts | | FT | 1 | String | Filter indicating if this genotype was “called” | | GL | G | Float | Genotype likelihoods | | GP | G | Float | Genotype posterior probabilities | | GQ | 1 | Integer | Conditional genotype quality | | GT | 1 | String | Genotype | | HQ | 2 | Integer | Haplotype quality | | MQ | 1 | Integer | RMS mapping quality | | PL | G | Integer | Phred-scaled genotype likelihoods rounded to the closest integer | | PP | G | Integer | Phred-scaled genotype posterior probabilities rounded to the closest integer | | PQ | 1 | Integer | Phasing quality | | PS | 1 | Integer | Phase set | - GT - `/` : genotype unphased 基因型未定相 - `|` : genotype phased 基因型已定相 <br> <hr> <br> ## VCF v.s. gVCF - [VCF和GVCF格式说明](https://www.cnblogs.com/timeisbiggestboss/p/9134733.html) ![](https://i.imgur.com/BLETrz8.png) - GVCF和VCF的最大区别是在于GVCF文件会记录所有的点，包括哪些没有突变的点。 - 在GVCF模式下，那些没有变异的点会形成一个未变异块，non-var block record。 - [GVCF - Genomic Variant Call Format](https://gatk.broadinstitute.org/hc/en-us/articles/360035531812-GVCF-Genomic-Variant-Call-Format) - VCF: Variant Call Format - GVCF: Genomic Variant Call Format <br> <hr> <br> ## 基因註釋類型 - [GWASLab – GWAS實驗室](https://gwaslab.com/author/cloufield/page/6/) - Variant Normalization 變異的標準化 - 使用ANNOVAR 對Variants進行功能註釋 Annotation POST-GWAS analysis - **基於基因的註釋 Gene-based annotation** 主要針對SNP或CNV是否引起蛋白編碼改變進行註釋，可以靈活選用 RefSeq genes, UCSC genes, ENSEMBL genes, GENCODE genes, AceView genes等多種不同來源的基因定義系統。 - **基於區域的註釋 Region-based annotation** 針對基因組某一特定區域的變異進行註釋，例如44個物種的保守區域，預測的轉錄因子結合位點，GWAS hit， ENCODE H3K4Me1/H3K4Me3/H3K27Ac/CTCF sites，ChIP-Seq peaks, RNA-Seq peaks等 - **基於篩選的註釋 Filter-based annotation** 使用某一特定的數據庫進行篩選註釋，例如註釋變異的rs id，1000基因組項目中的MAF，或是ExAC、gnomAD等，再例如SIFT/ PolyPhen/ LRT/ MutationTaster/ MutationAssessor/ FATHMM/ MetaSVM/ MetaLR 分數等。 <br> <hr> <br> ## 參考資料 - [VCF格式](https://www.jianshu.com/p/957efb50108f) - [变异信息那些事（上）/ VCF（Variant Call Format）是什么？](https://www.jianshu.com/p/3045aadbd723) - [SAM/BAM and related specifications](http://samtools.github.io/hts-specs/) - [[我們的基因體時代 OUR "GENE"RATION]:<br>VCF(Variant Call Format) 基因突變資料儲存格式](https://weitinglin.com/2017/05/29/vcfvariant-call-format-%E5%9F%BA%E5%9B%A0%E7%AA%81%E8%AE%8A%E8%B3%87%E6%96%99%E5%84%B2%E5%AD%98%E6%A0%BC%E5%BC%8F/) - [VCF详解](https://www.cnblogs.com/xiaofeiIDO/p/6610157.html) - [[ITRead01] 生物基因資料檔案——vcf格式詳解](https://www.itread01.com/content/1548698402.html) - [[GATK] VCF - Variant Call Format](https://software.broadinstitute.org/gatk/documentation/article?id=11005) - GATK: - [VCF - Variant Call Format](https://software.broadinstitute.org/gatk/documentation/article?id=11005) - [What is a VCF and how should I interpret it?](https://software.broadinstitute.org/gatk/documentation/article.php?id=1268) - [Phred-scaled Quality Scores](https://software.broadinstitute.org/gatk/documentation/article.php?id=4260) - [Math notes: How PL is calculated in HaplotypeCaller](https://gatkforums.broadinstitute.org/gatk/discussion/5913/math-notes-how-pl-is-calculated-in-haplotypecaller) - [VCF (Variant Call Format) version 4.0 | 1000 Genomes](https://www.internationalgenome.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-40) - [我們的基因體時代-AI, Data和生物資訊 Day18-基因變異的檔案格式VCF](https://ithelp.ithome.com.tw/articles/10267022) ![](https://i.imgur.com/L6JGvpN.png) #AA, AC, AF, AN, BQ, CIGAR, DB, DP, END, H2, H3, MQ, MQ0, NS, SB, SOMATIC, VALIDATED, 1000G