生物資訊資料存放型別 {#filetype}
##前言
各行各業都有在自己的標準體系,生物資訊學資料分析也不例外,各個廠商出品的晶片系列,還有各種NGS組學分析,都會涉及到不同的分析步驟,有著豐富多樣的中間檔案。其中一些常用的檔案就被規定成檔案格式。 檔案格式那麼多,都可以瞭解一二,當然,不需要背誦它們所有的細節,不過對下面我們單獨拿出來詳細介紹的還是儘量要耳熟能詳。
簡單來說,測序得到的是帶有質量值的鹼基序列(fastq格式);參考基因組是(fasta格式);用比對工具把fastq格式的序列回帖到對應的fasta格式的參考基因組序列,就可以產生sam格式的比對檔案;把sam格式的文字檔案壓縮成二進位制的bam檔案可以節省空間;如果對參考基因組上面的各個區段標記它們的性質,比如哪些區域是外顯子,內含子,UTR等等,這就是gtf/gff格式; 如果只是為了單純描述某個基因組區域,就是bed格式檔案,記錄染色體號以及起始終止座標,正負鏈即可;如果是記錄某些位點或者區域鹼基的變化,就是VCF檔案格式;如果僅僅是為了追蹤參考基因組的各個區域的覆蓋度,測序深度就可以用bigwig/wig格式。上面是對這些生信中的資料格式進行簡介,下面將對這些資料格式進行詳述。
FASTQ
簡介
FASTQ用於儲存生物序列(通常是核酸序列)和其測序質量資訊的標準格式。 其序列以及質量資訊都是使用一個ASCII字元標示,最初由Sanger開發。 目的是將FASTA序列與質量資料放到一起,目前已經成為高通量測序結果的實施標準。
一、定義和示例
FASTQ檔案中每個序列通常有四行:
第一行是序列標識以及相關的描述資訊,以‘@’開頭 第二行是序列 第三行以‘+’開頭,後面是序列標示符、描述資訊,或者什麼也不加,但是“+”不能少。 第四行,是質量資訊,和第二行的序列相對應,每一個序列都有一個質量評分,根據評分體系的不同,每個字元的含義表示的數字也不相同。
一個簡單的示例如下:
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65
格式口訣:
*‘@’開頭引標識*
*序列老二對老四*
*老三沒“+”不好使*
*老四質量分數最多事*
二、序列標識
上面說到第一行是序列標識以及相關的描述資訊,以‘@’開頭。可以像上面的示例那麼簡單,但如果是正規測序儀下機的真實資料,通常會很複雜。比如:
@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
這個序列標識以及相關描述資訊以冒號分割,每一個欄位資訊如下:
欄位 | 解釋 |
---|---|
EAS139 | the unique instrument name |
136 | the run id |
FC706VJ | the flowcell id |
2 | flowcell lane |
2104 | tile number within the flowcell lane |
15343 | ‘x’-coordinate of the cluster within the tile |
197393 | ‘y’-coordinate of the cluster within the tile |
1 | the member of a pair, 1 or 2 (paired-end or mate-pair reads only) |
Y | Y if the read fails filter (read is bad), N otherwise |
18 | 0 when none of the control bits are on, otherwise it is an even number |
ATCACG | index sequence |
當然,上面的表格介紹的只是其中一個測序儀下機資料,如果是其它機器,產商可以自由定義識別符號格式,因為fastq格式的第一行只需要以@符號開頭即可。
不過,也有一些時候fastq資料並不是測序儀直接下機的,而且他人上傳到了NCBI的SRA中心,我們下載下來解壓後一般就沒有了測序儀相關的標識,例子如下:
@SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACC
+SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IG9IC
三、質量編碼格式
質量評分指的是一個鹼基的錯誤概率的對數值。 其最初在Phred拼接軟體中定義與使用,其後在許多軟體中得到使用。 其質量得分與錯誤概率的對應關係見下表:
PHRED QUALITY SCORE | PROBABILITY OF INCORRECT BASE CALL | BASE CALL ACCURACY |
---|---|---|
10 | 1 in 10 | 90 % |
20 | 1 in 100 | 99% |
30 | 1 in 1000 | 99.9% |
40 | 1 in 10000 | 99.99% |
50 | 1 in 100000 | 99.999% |
Phred quality scores Q are defined as a property which is logarithmically related to the base-calling error probabilities P.
Q=-10lgP
除了Phred質量得分換算標準,還有就是Solexa標準:是把P換成p/(1-p)
對於每個鹼基的質量編碼標示,不同的軟體採用不同的方案,目前有5種方案:
- Sanger,Phred quality score,值的範圍從0到92,對應的ASCII碼從33到126,但是對於測序資料(raw read data)質量得分通常小於60,序列拼接或者mapping可能用到更大的分數。
- Solexa/Illumina 1.0, Solexa/Illumina quality score,值的範圍從-5到63,對應的ASCII碼從59到126,對於測序資料,得分一般在-5到40之間。
- Illumina 1.3+,Phred quality score,值的範圍從0到62對應的ASCII碼從64到126,低於測序資料,得分在0到40之間;
- Illumina 1.5+,Phred quality score,但是0到2作為另外的標示。
- Illumina 1.8+
下面是更為直觀的表示:
來自於 wikipedia:
四、檔案字尾
沒有特別的規定,通常使用.fq, .fastq, .txt等。 但是要注意,這個檔案格式主要指的是文字檔案裡面的每行每列的內容規則,並不是我們常見的計算機領域的mp3,mp4,avi,xls,doc等等。
其它注意事項:
- 雙端測序一般有兩個檔案(也可通過某種規則把兩個檔案合併成一個)。
- 第一個檔案與第二個檔案的行數完全一樣,且測序序列的排列順序完全一致。
- 在第一個檔案中,描述資訊的結尾是“/1”,表示是雙端測序的一端;第二個檔案中同樣位置/行數的相對應的測序序列的描述資訊則以“/2”結尾,表示是雙端測序的另一端。(2.2.2的表2-5中有敘述)
參考連結:
https://en.wikipedia.org/wiki/FASTQ_format
FASTA
簡介
FASTA格式用於表示核苷酸序列或氨基酸序列的格式。 在這種格式中鹼基對或氨基酸用單個字母來編碼,且允許在序列前新增序列名及註釋。 fasta序列格式是blast組織資料的基本格式,無論是資料庫還是查詢序列,大多數情況都使用fasta序列格式。 它要比上一小節介紹的FASTQ格式簡明很多。
一、定義和示例
總的來說,Fasta格式開始於一個識別符號:">",然後是一行描述,下面是的序列,直到下一個">",表示下一條。
下面是一個來源於NCBI的fasta格式序列:
>gi|187608668|ref|NM_001043364.2| Bombyx mori moricin (Mor), mRNA
AAACCGCGCAGTTATTTAAAATATGAATATTTTAAAACTTTTCTTTGTTTTTA
TTGTGGCAATGTCTCTGGTGTCATGTAGTACAGCCGCTCCAGCAAAAATACCT
ATCAAGGCCATTAAGACTGTAGGAAAGGCAGTCGGTAAAGGTCTAAGAGCCAT
ATCAAGGCCATTAA
Fasta格式首先以大於號“>”開頭,接著是序列的識別符號“gi|187608668|ref|NM_001043364.2|”,然後是序列的描述資訊。 換行後是序列資訊,序列中允許空格,換行,空行,直到下一個大於號,表示該序列的結束。 下面簡單給一個表格說明序列來源的資料庫與對應的識別符號
Database Name資料庫名稱 | Identifier Syntax 識別符號 |
---|---|
GenBank | ```gb |
EMBL Data Library | ```emb |
DDBJ, DNA Database of Japan | ``` dbj |
NBRF PIR | ```pir |
Protein Research Foundation | ``` prf |
SWISS-PROT | ```sp |
Brookhaven Protein Data Bank | ```pdb |
Patents | ```pat |
GenInfo Backbone Id | ``` bbs |
General database identifier | ``` gnl |
NCBI Reference Sequence | ``` ref |
Local Sequence identifier | ```lcl |
通常情況下序列的識別符號不會像上面的例子那樣複雜,再複雜的識別符號也是有規則的,上面的識別符號是NCBI定義的,可以去其官網瞭解詳情。
自己總結了一條速記:
*大於號來表開頭,*
*描述緊跟在後頭,*
*ABCDEFG換行組成序列呦~*
二、序列中字母代表的含義
FASTA格式支援的核苷酸程式碼如下:
核苷酸程式碼 | 意義 |
---|---|
A | Adenosine |
C | Cytosine |
G | Guanine |
T | Thymidine |
U | Uracil |
R | G A (puRine) |
Y | T C (pYrimidine) |
K | G T (Ketone) |
M | A C (aMino group) |
S | G C (Strong interaction) |
W | A T (Weak interaction) |
B | G T C (not A) (B comes after A) |
D | G A T (not C) (D comes after C) |
H | A C T (not G) (H comes after G) |
V | G C A (not T, not U) (V comes after U) |
N | A G C T (aNy) |
X | masked |
- | gap of indeterminate length |
FASTA格式支援的氨基酸程式碼如下:
氨基酸程式碼 | 意義 |
---|---|
A | Alanine |
B | Aspartic acid or Asparagine |
C | Cysteine |
D | Aspartic acid |
E | Glutamic acid |
F | Phenylalanine |
G | Glycine |
H | Histidine |
I | Isoleucine |
K | Lysine |
L | Leucine |
M | Methionine |
N | Asparagine |
O | Pyrrolysine |
P | Proline |
Q | Glutamine |
R | Arginine |
S | Serine |
T | Threonine |
U | Selenocysteine |
V | Valine |
W | Tryptophan |
Y | Tyrosine |
Z | Glutamic acid or Glutamine |
X | any |
* | translation stop |
- | gap of indeterminate length |
參考連結 https://en.wikipedia.org/wiki/FASTA_format
注意事項:
對於自己構建的序列資料庫(序列不是來源與NCBI或其他資料),可以採用“gnl|database|identifier”或者“lcl|identifier”格式,以保證可以使用blast的所有功能。 database或者identifier是需要指定的資料庫的標識和序列標識,指定的名稱可以用大小寫字母、數字、下劃線“_”、破折號“-”或者點號“.”。 注意名稱是區分大小寫的,同時不能出現空格,空格表示序列識別符號結束。
資料庫中的序列識別符號必須保證唯一,許多時候格式資料庫是formatdb報告錯誤,就是因為標示符重複,還有一點需要強調的是序列不能為空,否則也會報錯。
SAM格式
簡介
SAM格式主要應用於測序序列mapping到基因組上的結果表示,當然也可以表示任意的多重比對結果。 SAM是一種序列比對格式標準,由sanger制定,是以TAB為分割符的文字格式。 SAM的全稱是sequence alignment/map format。
一、定義和示例
SAM分為兩部分,註釋資訊(header section )和比對結果部分 (alignment section)。 通常是把FASTQ檔案格式的測序資料比對到對應的參考基因組版本得到的。 註釋資訊並不是SAM檔案的重點,是該SAM檔案產生以及被處理過程的一個記錄,規定以@開頭,用不同的tag表示不同的資訊,主要有:
- @HD,說明符合標準的版本、對比序列的排列順序;
- @SQ,參考序列說明;
- @RG,比對上的序列(read)說明;
- @PG,使用的程式說明;
- @CO,任意的說明資訊。
一個簡單的SAM檔案例子如下:
@HD VN:1.0 SO:unsorted
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@PG ID:Bowtie VN:1.0.0 CL:"bowtie genome/hg19 -q reads/SRR3101251.fastq -m 1 -p 4 -S"
SRR3101251.1 0 chr19 9486878 255 49M * 0 0 NTACTCCCACTACTCTCAGATTCAAGCAATCCTCCCACCCTAGCCCACC #1=DDDFFHHHHHIHHIJJJHIJIIJIHIFHJIIJJJJJJJIIJJJJJJ XA:i:1 MD:Z:0A48 NM:i:1
SRR3101251.5 16 chr2 240279787 255 49M * 0 0 CCTGAATCCATCAGAGCAGCCGGGCTGTGACACTCACTGTCATGATGTT JIJJIHIIIIJJJJJJJJJGHJJJJIIHJHICJIGCHHHHHFFFFFCCC XA:i:0 MD:Z:49 NM:i:0
SRR3101251.6 4 * 0 0 * * 0 0 NATTCCCACCTATGAGTGAGAATATGCGGTGTTTGGTTTTTTGTTCTTG #1=DDDFFHHHHHJJJGHIJJJJJJJJJJCGGIIJJIIJJJIJHJIIJJ XM:i:1
前四行是註釋資訊,其後是比對結果,下面對比對結果進行解釋,它是SAM格式檔案的精華部分。
二、比對結果詳解
SRR035022.2621862 163 16 59999 37 22S54M = 60102 179 CCAACCCAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCGACCCTCACCCTCACCC >AAA=>?AA>@@[email protected][email protected][email protected]@A?A>[email protected]?AAAAB??ABAB?79A?AAB;[email protected][email protected]<=8:8 XT:A:M XN:i:2 SM:i:37 AM:i:37 XM:i:0 XO:i:0 XG:i:0 RG:Z:SRR035022 NM:i:2 MD:Z:0N0N52 OQ:Z:[email protected];[email protected]
所以在我們的例子中,每一個欄位的說明如下:
QNAME SRR035022.2621862
FLAG 163
RNAME 16
POS 59999
MAQ 37
CIGAR 22S54M
MRNM =
MPOS 60102
ISIZE 179
SEQ CCAACCCAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCGACCCTCACCCTCACCC
QUAL >AAA=>?AA>@@[email protected][email protected][email protected]@A?A>[email protected]?AAAAB??ABAB?79A?AAB;[email protected][email protected]<=8:8
TAG XT:A:M
TAG XN:i:2
TAG SM:i:37
TAG AM:i:37
TAG XM:i:0
TAG XO:i:0
TAG XG:i:0
TAG RG:Z:SRR035022
TAG NM:i:2
TAG MD:Z:0N0N52
TAG OQ:Z:[email protected];CCCBBC?CCCACCACA
比對結果部分(alignment section),每一行表示一個片段(segment)的比對資訊,包括11個必須的欄位(mandatory fields)和一個可選的欄位,欄位之間用tag分割。
必須的欄位有11個,順序固定,不可自行改動,根據欄位定義,可以為’0‘或者’*‘,這是11個欄位包括:
-
QNAME,比對片段的(template)的編號;
-
FLAG,位標識,template mapping情況的數字表示,每一個數字代表一種比對情況,這裡的值是符合情況的數字相加總和; (picard專門有一個工具解讀sam的flag:http://broadinstitute.github.io/picard/explain-flags.html)
1 The read is one of a pair read是pair中的一條(read表示本條read,mate表示pair中的另一條read)
2 The alignment is one end of a proper paired-end alignment pair一正一負完美的比對上
4 The read has no reported alignments 這條read沒有比對上
8 The read is one of a pair and has no reported alignments mate沒有比對上
16 The alignment is to the reverse reference strand 這條read反向比對
32 The other mate in the paired-end alignment is aligned to the reverse reference strand mate反向比對
64 The read is the first (#1) mate in a pair 這條read是read1
128 The read is the second (#2) mate in a pair 這條read是read2
-
RNAME,參考序列的編號,如果註釋中對SQ-SN進行了定義,這裡必須和其保持一致,另外對於沒有mapping上的序列,這裡是’*‘;
-
POS,比對上的位置,注意是從1開始計數,沒有比對上,此處為0;
-
MAPQ,mappint的質量;
-
CIGAR,簡要比對資訊表示式(Compact Idiosyncratic Gapped Alignment Report),其以參考序列為基礎,使用數字加字母表示比對結果,比如3S6M1P1I4M,前三個鹼基被剪下去除了,然後6個比對上了,然後打開了一個缺口,有一個鹼基插入,最後是4個比對上了,是按照順序的;
“M”表示 match或 mismatch;
“I”表示 insert;
“D”表示 deletion;
“N”表示 skipped(跳過這段區域);
“S”表示 soft clipping(被剪下的序列存在於序列中);
“H”表示 hard clipping(被剪下的序列不存在於序列中);
“P”表示 padding;
“=”表示 match;
“X”表示 mismatch(錯配,位置是一一對應的);
-
RNEXT,下一個片段比對上的參考序列的編號,沒有另外的片段,這裡是’*‘,同一個片段,用’=‘;
-
PNEXT,下一個片段比對上的位置,如果不可用,此處為0;
-
TLEN,Template的長度,最左邊得為正,最右邊的為負,中間的不用定義正負,不分割槽段(single-segment)的比對上,或者不可用時,此處為0;
-
SEQ,序列片段的序列資訊,如果不儲存此類資訊,此處為’*‘,注意CIGAR中M/I/S/=/X對應數字的和要等於序列長度;
-
QUAL,序列的質量資訊,格式同FASTQ一樣。read質量的ASCII編碼。
-
可選欄位(optional fields),格式如:TAG:TYPE:VALUE,其中TAG有兩個大寫字母組成,每個TAG代表一類資訊,每一行一個TAG只能出現一次,TYPE表示TAG對應值的型別,可以是字串、整數、位元組、陣列等。
三、SAM要處理好的問題
- 非常多序列(read),mapping到多個參考基因組(reference)上
- 同一條序列,分多段(segment)比對到參考基因組上
- 無限量的,結構化資訊表示,包括錯配、刪除、插入等比對資訊
要注意的幾個概念,以及與之對應的模型:
- reference
- read
- segment
- template(參考序列和比對上的序列共同組成的序列為template)
- alignment
- seq
參考連結:
https://en.wikipedia.org/wiki/SAM_(file_format)
BAM 格式
簡介
本質上就是二進位制壓縮的SAM檔案,大部分生物資訊學流程都需要這個格式,為了節省儲存空間以及方便索引。
# BiocInstaller::biocLite('Rsamtools')
library(Rsamtools)
test_bam_file <- 'data/CHIP-seq.bam'
#fileter bam
filter <- FilterRules(list(MinWidth = function(x) width(x$seq) > 35))
res <- scanBam(test_bam_file, filter=filter)[[1]]
sapply(res, head)
從上面的例子(需要在R裡面執行)可以看到BAM檔案需要用特殊的方法來讀取,可以是R裡面的Rsamtools包,也可以是linux環境下安裝好的samtools軟體,因為它是二進位制檔案,不能像普通的文字檔案那樣來開啟。
我們用R裡面的head函式查看了該BAM檔案的前6行,比對的flag分別是16 0 16 16 0 0
,說明有3條序列沒有成功比對到基因組。width資訊說明該序列長度都是36bp。序列的鹼基以及對應的鹼基質量也如上所述。
VCF
簡介
Variant Call Format(VCF)是一個用於儲存基因序列突變資訊的文字格式。 可以表示單鹼基突變, 插入/缺失, 拷貝數變異和結構變異等。 通常是對BAM檔案格式的比對結果進行處理得到的。 BCF格式檔案是VCF格式的二進位制檔案。我們就不再介紹BCF格式啦。
提到vcf就必須提到千人基因組計劃,因為千人計劃組才產生的vcf。生信菜鳥團有一篇部落格《居然可以下載千人基因組計劃的所有資料bam,vcf資料》 http://www.bio-info-trainee.com/1339.html 專門講了千人計劃的資料下載。
當然,所有的資料格式定義,都推薦大家看原汁原味的英文介紹,那個是金標準,我們的翻譯只是為了促進大家的理解,如果有模稜兩可的地方,以英文原文為準:https://samtools.github.io/hts-specs/VCFv4.2.pdf
一、定義和示例
如上圖所示,vcf記錄的即為各型別的變異。例如:點突變,拷貝數變異,插入,缺失等結構變異。
VCF分為兩部分,註釋資訊和變異位點記錄資訊。
註釋資訊通常以#開頭,會描述該VCF版本,目前以4.2居多,然後會一行行記錄變異位點資訊裡面會出現的所有TAG。
下面這個是NCBI的dbSNP資料庫裡面的人類的vcf檔案的部分擷取:
##fileformat=VCFv4.0
##fileDate=20160601
##source=dbSNP
##dbSNP_BUILD_ID=147
##reference=GRCh37.p13
##phasing=partial
##variationPropertyDocumentationUrl=ftp://ftp.ncbi.nlm.nih.gov/snp/specs/dbSNP_BitField_latest.pdf
##INFO=<ID=RS,Number=1,Type=Integer,Description="dbSNP ID (i.e. rs number)">
##INFO=<ID=RSPOS,Number=1,Type=Integer,Description="Chr position reported in dbSNP">
##INFO=<ID=RV,Number=0,Type=Flag,Description="RS orientation is reversed">
#CHROM POS ID REF ALT QUAL FILTER INFO
1 10177 rs201752861 A C . . RS=201752861;RSPOS=10177;dbSNPBuildID=137;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP
1 10177 rs367896724 A AC . . RS=367896724;RSPOS=10177;dbSNPBuildID=138;SSR=0;SAO=0;VP=0x050000020005170026000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP;VLD;G5A;G5;KGPhase3;CAF=0.5747,0.4253;COMMON=1
1 10228 rs143255646 TA T . . RS=143255646;RSPOS=10229;dbSNPBuildID=134;SSR=0;SAO=0;VP=0x050000020005000002000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP
1 10228 rs200462216 TAACCCCTAACCCTAACCCTAAACCCTA T . . RS=200462216;RSPOS=10229;dbSNPBuildID=137;SSR=0;SAO=0;VP=0x050000020005000002000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP
1 10230 rs775928745 AC A . . RS=775928745;RSPOS=10231;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000020005000002000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP
1 10231 rs200279319 C A . . RS=200279319;RSPOS=10231;dbSNPBuildID=137;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP
1 10234 rs145599635 C T . . RS=145599635;RSPOS=10234;dbSNPBuildID=134;SSR=0;SAO=0;VP=0x050100020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;SLO;R5;ASP
1 10235 rs540431307 T TA . . RS=540431307;RSPOS=10235;dbSNPBuildID=142;SSR=0;SAO=0;VP=0x050000020005040024000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP;VLD;KGPhase3;CAF=0.9988,0.001198;COMMON=0
1 10247 rs796996180 T C . . RS=796996180;RSPOS=10247;dbSNPBuildID=146;SSR=0;SAO=0;VP=0x050100020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;SLO;R5;ASP
限於文章篇幅限制,我只是截取了該VCF檔案的部分註釋資訊,很明顯可以看到註釋資訊剛剛開始的幾行其實是沒有規則的,只需要以##開頭即可,描述一些必備資訊,包括參考基因組版本,得到該VCF檔案的命令是什麼等等。
後面的都是以INFO=<ID=······>的形式來介紹一個個TAG,這些TAG都是會在VCF的正文,變異位點記錄裡面用到的。而且每個tag都很容易理解,就是對應的英文描述而已。
接下來我們看看比較複雜的正文部分,就是變異位點記錄資訊。
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample1
1 858691 . TG T 222 . INDEL;IDV=37;IMF=0.486842;DP=76;VDB=0.110516;SGB=-0.693139;MQSB=1;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=12,24,14,22;MQ=60 GT:PL 0/1:255,0,255
1 858801 . A G 222 . DP=59;VDB=0.728126;SGB=-0.692717;RPB=0.748623;MQB=1;MQSB=1;BQB=0.963908;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=14,17,11,12;MQ=60 GT:PL 0/1:255,0,255
1 859404 . C G 222 . DP=81;VDB=0.0896228;SGB=-0.693132;RPB=0.849598;MQB=1;MQSB=1;BQB=0.486963;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=17,15,18,16;MQ=60 GT:PL 0/1:255,0,255
1 859690 . C G 222 . DP=75;VDB=0.0662538;SGB=-0.69312;RPB=0.959181;MQB=1;MQSB=1;BQB=0.962588;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=20,15,18,14;MQ=60 GT:PL 0/1:255,0,255
1 859701 . C G 222 . DP=74;VDB=0.274853;SGB=-0.693127;RPB=0.97201;MQB=1;MQSB=1;BQB=0.717302;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=19,15,19,14;MQ=60 GT:PL 0/1:255,0,255
1 859913 . A G 222 . DP=67;VDB=0.756546;SGB=-0.693139;RPB=0.950685;MQB=1;MQSB=1;BQB=0.662934;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=18,10,19,17;MQ=60 GT:PL 0/1:255,0,255
1 860416 . G A 222 . DP=79;VDB=0.673886;SGB=-0.693144;RPB=0.11919;MQB=1;MQSB=1;BQB=0.992984;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=18,15,24,15;MQ=60 GT:PL 0/1:255,0,255
每一行代表一個Variant的資訊。
CHROM 和 POS:代表參考序列名和variant的位置;如果是INDEL的話,位置是INDEL的第一個鹼基位置。
ID:variant的ID。比如在dbSNP中有該SNP的id,則會在此行給出(這個需要自己下載dbSNP資料庫檔案進行註釋才有的)。 若沒有或者註釋不上,則用’.'表示其為一個novel variant。
REF 和 ALT:參考序列的鹼基 和 Variant的鹼基。
QUAL:Phred格式(Phred_scaled)的質量值,表示在該位點存在variant的可能性;該值越高,則variant的可能性越大; 計算方法:Phred值 = -10 * log (1-p) p為variant存在的概率; 通過計算公式可以看出值為10的表示錯誤概率為0.1,該位點為variant的概率為90%。
FILTER:使用上一個QUAL值來進行過濾的話,是不夠的。GATK能使用其它的方法來進行過濾,過濾結果中通過則該值為”PASS”;若variant不可靠,則該項不為”PASS”或”.”。
INFO: 這一行是variant的詳細資訊,內容很多,以下再具體詳述。
FORMAT 和 sample1 :這兩行合起來提供了 sample1 這個sample的基因型的資訊。’sample1′代表這該名稱的樣品,是由SAM/BAM檔案中的@RG下的 SM 標籤決定的。 (當然並不是所有的VCF都是由一個BAM檔案產生,比如資料庫dbSNP提供的vcf檔案,就沒有樣本資訊啦)
1、第8列的INFO
該列資訊最多了,都是以 “TAG=Value”, 並使用”;”分隔的形式 。 其中 很多的註釋資訊在VCF檔案的頭部註釋中給出。以下是這些TAG的解釋
AC,AF 和 AN:AC(Allele Count) 表示該Allele的數目;AF(Allele Frequency) 表示Allele的頻率; AN(Allele Number) 表示Allele的總數目。對於1個diploid sample而言:則基因型 0/1 表示sample為雜合子,Allele數為1(雙倍體的sample在該位點只有1個等位基因發生了突變),Allele的頻率為0.5(雙倍體的 sample在該位點只有50%的等位基因發生了突變),總的Allele為2; 基因型 1/1 則表示sample為純合的,Allele數為2,Allele的頻率為1,總的Allele為2。
DP:reads覆蓋度。是一些reads被過濾掉後的覆蓋度。
Dels:Fraction of Reads Containing Spanning Deletions。進行SNP和INDEL calling的結果中,有該TAG並且值為0表示該位點為SNP,沒有則為INDEL。
FS:使用Fisher’s精確檢驗來檢測strand bias而得到的Fhred格式的p值。該值越小越好。一般進行filter的時候,可以設定 FS < 10~20。
HaplotypeScore:Consistency of the site with at most two segregating haplotypes
InbreedingCoeff:Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hard-Weinberg expectation
MLEAC:Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), 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
MQ:RMS Mapping Quality
MQ0:Total Mapping Quality Zero Reads
MQRankSum:Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities
QD:Variant Confidence/Quality by Depth
RPA:Number of times tandem repeat unit is repeated, for each allele (including reference)
RU:Tandem repeat unit (bases)
ReadPosRankSum:Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias
STR:Variant is a short tandem repeat
2、9和10列代表基因型
GT:AD:DP:GQ:PL 0/1:173,141:282:99:255,0,255
GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26
看上面最後兩列資料,這兩列資料是對應的,前者為格式,後者為格式對應的資料。這些TAG也是可以在VCF的標頭檔案找到的
GT:樣品的基因型(genotype)。兩個數字中間用’/'分開,這兩個數字表示雙倍體的sample的基因型。0 表示樣品中有ref的allele; 1 表示樣品中variant的allele; 2表示有第二個variant的allele。因此: 0/0 表示sample中該位點為純合的,和ref一致; 0/1 表示sample中該位點為雜合的,有ref和variant兩個基因型; 1/1 表示sample中該位點為純合的,和variant一致。
AD 和 DP:AD(Allele Depth)為sample中每一種allele的reads覆蓋度,在diploid中則是用逗號分割的兩個值,前者對應ref基因型,後者對應variant基因型。
DP(Depth)為sample中該位點的測序深度。
GQ:基因型的質量值(Genotype Quality)。Phred格式(Phred_scaled)的質量值,表示在該位點該基因型存在的可能性;該值越高,則Genotype的可能性越大;計算方法:Phred值 = -10 * log (1-p) p為基因型存在的概率。
PL:指定的三種基因型的質量值(provieds the likelihoods of the given genotypes)。這三種指定的基因型為(0/0,0/1,1/1),這三種基因型的概率總和為1。和之前不一致,該值越大,表明為該種基因型的可能性越小。 Phred值 = -10 * log (p) p為基因型存在的概率。
上圖可以幫助很好的理解vcf格式。
注意事項: 針對vcf格式有如bcftools等軟體進行處理。
生信菜鳥團的部落格和生信技能樹裡面都介紹了很多處理vcf的軟體。
強烈推薦去看《生信技能樹》中的:我的基因組28 https://en.wikipedia.org/wiki/Variant_Call_Format
GTF和GFF
簡介
GFF全稱為general feature format,這種格式主要是用來註釋基因組。 GTF全稱為gene transfer format,主要是用來對基因進行註釋。 GTF和GFF格式是Sanger研究所定義,是一種簡單的、方便的對於DNA、RNA以及蛋白質序列的特徵進行描述的一種資料格式. 比如序列的哪裡到哪裡是基因,是轉錄本,是外顯子,內含子或者CDS等等,已經成為序列註釋的通用格式,許多軟體都支援輸入或者輸出gff格式。
一、定義和示例
gff由tab鍵隔開的9列組成,以下是各列的說明:
Column 1: “seqid”
序列的編號,編號的有效字元[a-zA-Z0-9.:^*[email protected]!+_?-|]
Column 2: “source”
註釋資訊的來源,比如”Genescan”、”Genbank” 等,可以為空,為空用”.”點號代替
Column 3: “type”
註釋資訊的型別,比如Gene、cDNA、mRNA等,或者是SO對應的編號
Columns 4 & 5: “start” and “end”
開始與結束的位置,注意計數是從1開始的。結束位置不能大於序列的長度
Column 6: “score”
得分,數字,是註釋資訊可能性的說明,可以是序列相似性比對時的E-values值或者基因預測是的P-values值。”.”表示為空。
Column 7: “strand”
序列的方向, +表示正義鏈, -反義鏈 , ? 表示未知.
Column 8: “phase”
僅對註釋型別為 “CDS”有效,表示起始編碼的位置,有效值為0、1、2。
Column 9: “attributes”
以多個鍵值對組成的註釋資訊描述,鍵與值之間用”=“,不同的鍵值用”;“隔開,一個鍵可以有多個值,不同值用”,“分割。注意如果描述中包括tab鍵以及”,=;”,要用URL轉義規則進行轉義,如tab鍵用 %09代替。鍵是區分大小寫的,以大寫字母開頭的鍵是預先定義好的,在後面可能被其他註釋資訊所呼叫。
預先定義的鍵包括:
ID 註釋資訊的編號,在一個GFF檔案中必須唯一;
Name 註釋資訊的名稱,可以重複;
Alias 別名
Parent Indicates 該註釋所屬的註釋,值為註釋資訊的編號,比如外顯子所屬的轉錄組編號,轉錄組所屬的基因的編號。值可以為多個。
Target Indicates: the target of a nucleotide-to-nucleotide or protein-to-nucleotide alignment.(核苷酸對核苷酸或蛋白質至核苷酸比對的靶點。)
Gap:The alignment of the feature to the target if the two are not collinear (e.g. contain gaps).(如果兩者不共線(例如包含間隙),則該特徵與目標的對準)
Derives_from:Used to disambiguate the relationship between one feature and another when the relationship is a temporal one rather than a purely structural “part of” one. This is needed for polycistronic genes.(用於消除一個特徵與另一個特徵之間的關係,當關系是一個時間的關係,而不是純粹的結構“一部分”時。 這是多順反子基因所必需的。)
Note: 備註
Dbxref :資料庫索引
Ontology_term: A cross reference to an ontology term.
下面是一個簡單的例項:
##gff-version 3
##sequence-region ctg123 1 1497228
ctg123 . gene 1000 9000 . + . ID=gene00001;Name=EDEN
ctg123 . TF_binding_site 1000 1012 . + . Parent=gene00001
ctg123 . mRNA 1050 9000 . + . ID=mRNA00001;Parent=gene00001
ctg123 . mRNA 1050 9000 . + . ID=mRNA00002;Parent=gene00001 ctg123 . mRNA 1300 9000 . + . ID=mRNA00003;Parent=gene00001 ctg123 . exon 1300 1500 . + . Parent=mRNA00003
ctg123 . exon 1050 1500 . + . Parent=mRNA00001,mRNA00002
ctg123 . exon 3000 3902 . + . Parent=mRNA00001,mRNA00003
ctg123 . exon 5000 5500 . + . Parent=mRNA00001,mRNA00002,mRNA00003
ctg123 . exon 7000 9000 . + . Parent=mRNA00001,mRNA00002,mRNA00003
ctg123 . CDS 1201 1500 . + 0 ID=cds00001;Parent=mRNA00001
ctg123 . CDS 3000 3902 . + 0 ID=cds00001;Parent=mRNA00001
ctg123 . CDS 5000 5500 . + 0 ID=cds00001;Parent=mRNA00001
ctg123 . CDS 7000 7600 . + 0 ID=cds00001;Parent=mRNA00001
ctg123 . CDS 1201 1500 . + 0 ID=cds00002;Parent=mRNA00002
ctg123 . CDS 5000 5500 . + 0 ID=cds00002;Parent=mRNA00002
ctg123 . CDS 7000 7600 . + 0 ID=cds00002;Parent=mRNA00002
ctg123 . CDS 3301 3902 . + 0 ID=cds00003;Parent=mRNA00003
ctg123 . CDS 5000 5500 . + 2 ID=cds00003;Parent=mRNA00003
ctg123 . CDS 7000 7600 . + 2 ID=cds00003;Parent=mRNA00003
ctg123 . CDS 3391 3902 . + 0 ID=cds00004;Parent=mRNA00003
ctg123 . CDS 5000 5500 . + 2 ID=cds00004;Parent=mRNA00003
ctg123 . CDS 7000 7600 . + 2 ID=cds00004;Parent=mRNA00003
可以看到第9列其實是可以無限擴充套件的,只需要以封號進行分割即可。
二、GTF 與GFF的差異
GTF檔案以及GFF檔案都由9列資料組成,這兩種檔案的前8列都是相同的(一些小的差別),它們的差別重點在第9列。
GTF檔案的第9列同GFF檔案不同,雖然同樣是標籤與值配對的情況,但標籤與值之間以空格分開,而不是GFF裡面的=號 且每個特徵之後都要有分號,(包括最後一個特徵). 下面看一個GTF的例項:
17 havana five_prime_utr 7687377 7687427 . - . gene_id "ENSG00000141510"; gene_version "16"; transcript_id "ENST00000503591"; transcript_version "1"; gene_name "TP53"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; havana_gene "OTTHUMG00000162125"; havana_gene_version "10"; transcript_name "TP53-003"; transcript_source "havana"; transcript_biotype "protein_coding"; havana_transcript "OTTHUMT00000367399"; havana_transcript_version "2"; tag "cds_end_NF"; tag "mRNA_end_NF"; transcript_support_level "5";
17 havana five_prime_utr 7677325 7677427 . - . gene_id "ENSG00000141510"; gene_version "16"; transcript_id "ENST00000503591"; transcript_version "1"; gene_name "TP53"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; havana_gene "OTTHUMG00000162125"; havana_gene_version "10"; transcript_name "TP53-003"; transcript_source "havana"; transcript_biotype "protein_coding"; havana_transcript "OTTHUMT00000367399"; havana_transcript_version "2"; tag "cds_end_NF"; tag "mRNA_end_NF"; transcript_support_level "5";
17 havana five_prime_utr 7676595 7676622 . - . gene_id "ENSG00000141510"; gene_version "16"; transcript_id "ENST00000503591"; transcript_version "1"; gene_name "TP53"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; havana_gene "OTTHUMG00000162125"; havana_gene_version "10"; transcript_name "TP53-003"; transcript_source "havana"; transcript_biotype "protein_coding"; havana_transcript "OTTHUMT00000367399"; havana_transcript_version "2"; tag "cds_end_NF"; tag "mRNA_end_NF"; transcript_support_level "5";
速記口訣
* GTF,GFF 9列資料來組成*
*結果前8列都相同*
*GFF 9列間以tab 來分割,第9列用等號*
*GTF不服氣,標籤與值用空格分離,特徵之間也用分號分開;要問特立獨行誰最牛,GTF排第一!*
有趣小故事:
從前雙胞胎兄弟,GTF,GFF。GFF是哥哥,GTF是弟弟,性別相同,都是由用來註釋基因的,都是9個器官~~哥哥呢,比較隨和,哥哥器官就用tab分開,然後第九個一對就劃等號!
GTF是弟弟,天天心裡想”WTF“! 不服氣,就要特立獨行!然後不服氣,9個器官有八個都與哥哥相同,第9個子集偏不,就用了空格把它分成兩塊~~各個器官之間也是用分號分開!!預示子集跟各個不一樣,自己的獨立性!
目前兩種檔案可以方便的相互轉化,比如:使用Cufflinks軟體的 的gffread。
參考連結 http://genome.ucsc.edu/goldenPath/help/customTrack.html#GTFhttp://blog.sina.com.cn/s/blog_8a4f556e0102yd3l.html
GenePred
簡介
GenePred格式主要用於基因瀏覽器中基因預測的track。 如果有可變剪下的情況,那表格的每一行就是一個 transcript 的全部資訊。
一、定義和示例
每一行的具體解釋如下
table genePred
"A gene prediction."
(
string name; "Name of gene"
string chrom; "Chromosome name"
char[1] strand; "+ or - for strand"
uint txStart; "Transcription start position"
uint txEnd; "Transcription end position"
uint cdsStart; "Coding region start"
uint cdsEnd; "Coding region end"
uint exonCount; "Number of exons"
uint[exonCount] exonStarts; "Exon start positions"
uint[exonCount] exonEnds; "Exon end positions"
)
如果覺得抽象,我們可以用示例來進行一下對比。小編在這裡首先將模式植物擬南芥的gtf檔案轉化為gpd格式。head -n 1
看一下gpd檔案第一行的樣式。
#gpd
AT1G01010.1 Chr1 + 3630 5899 3759 5630 6 3630,3995,4485,4705,5173,5438, 3913,4276,4605,5095,5326,5899
再來grep
一下gtf檔案中有AT1G01010.1資訊的那些行是什麼樣
#gtf
Chr1 Araport11 5UTR 3631 3759 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 exon 3631 3913 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 start_codon 3760 3762 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 CDS 3760 3913 . + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 exon 3996 4276 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 CDS 3996 4276 . + 2 transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 exon 4486 4605 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 CDS 4486 4605 . + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 exon 4706 5095 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 CDS 4706 5095 . + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 exon 5174 5326 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 CDS 5174 5326 . + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 CDS 5439 5630 . + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 exon 5439 5899 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 stop_codon 5628 5630 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
Chr1 Araport11 3UTR 5631 5899 . + . transcript_id "AT1G01010.1"; gene_id "AT1G01010"
我們可以看出,在GTF中,AT1G01010.1這個transcript共有6個CDS,那麼對應到相應gpd檔案AT1G01010.1這一行的第8列就是6,而第9列和第10列則是這6個CDS對應的起始和終止位置。
細心的朋友可能會發現,GTF檔案中CDS起始位置在GenePred table中統統少了1,這其實就是兩種檔案的起始位置從1開始還是從0開始計數的區別。
GTF的產生和流行有其歷史的原因。但是從技術角度來講,這個檔案格式是個非常差勁的格式。
GTF格式非常冗餘。以人類轉錄組為例,Gencode V22的GTF檔案為1.2G,壓縮之後只有40M。大家知道壓縮軟體的壓縮比是和軟體的冗餘程度。很少有檔案能夠壓縮到1/30的大小。可見GTF格式多麼冗餘。GTF格式(及其早期版本GFF等)有很好的替代格式。從資訊量上來講:GTF 等價於 GenePred (Bed12) + Gene_Anno_table。
GenePred是Jimmy Kent建立UCSC genome browser的時候建立的檔案格式。UCSC的檔案格式定義是非常smart的,包括之後我可能會講到的2bit,bigwig格式。
二、GTF vs GenePred
- 從檔案大小上來看,壓縮前:GTF(1.2G) >> Genepred(23M) + Gene_Anno_table (2.8M)。壓縮後:GTF(40M) >> GenePred(7.8M) +Gene_Anno_table (580K)
- 從可讀性來講,GTF是以gene interval 為單位(行),每行可以是gene,transcript,exon,intron,utr等各種資訊,看起來什麼都在裡面,很全面,其實可讀性非常差,而且容易產生各種錯誤。GenePred格式是以transcript為單位,每行就是一個transcript,簡潔直觀。
- 從程式處理的角度來講:以GTF檔案作為輸入的程式,如果換成以GenePred格式為輸入,程式設計的難度會降低一個數量級,運算時間會快很多,程式碼的可讀性強很多。
GTF 轉換成GenePred:
gtfToGenePred -genePredExt -ignoreGroupsWithoutExons -geneNameAsName2 test.gtf test.gpd
Gene_Anno_table: 其實就是把GTF的所有transcript行的第9列轉換變成一個表格。
BED
簡介
BED 檔案格式提供了一種靈活的方式來定義的資料行,用於描述註釋的資訊。 跟GTF/GFF格式一樣,也可以用來描述基因組特徵。但沒有GTF/GFF格式那麼正規,通常用來描述 任何人為定義的區間。 但沒有GTF/GFF格式那麼正規,通常用來描述任何人為定義的區間。 所以BED格式最重要的就是染色體加上起始終止座標這3列。
一、定義和示例
BED行有3個必須的列和9個額外可選的列。 每行的資料格式要求一致。
1、必須包含的3列是:
-
chrom, 染色體或scafflold 的名字(eg chr3, chrY, chr2_random, scaffold0671 )
-
chromStart 染色體或scaffold的起始位置,染色體第一個鹼基的位置是0
-
chromEnd 染色體或scaffold的結束位置,染色體的末端位置沒有包含到顯示資訊裡面。例如,首先得100個鹼基的染色體定義為chromStart =0 . chromEnd=100, 鹼基的數目是0-99
2、9 個額外的可選列:
-
name 指定BED行的名字,這個名字標籤會展示在基因組瀏覽器中的bed行的左側。
-
score 0到1000的分值,如果在註釋資料的設定中將原始基線設定為1,那麼這個分值會決定現示灰度水平(數字越大,灰度越高)
-
strand 定義鏈的方向,''+” 或者”-”
-
thickStart 起始位置(The starting position at which the feature is drawn thickly)(例如,基因起始編碼位置)
-
thickEnd 終止位置(The ending position at which the feature is drawn thickly)(例如:基因終止編碼位置)
-
itemRGB 是一個RGB值的形式, R, G, B (eg. 255, 0, 0), 如果itemRgb設定為'On”, 這個RBG值將決定資料的顯示的顏色。
-
blockCount BED行中的block數目,也就是外顯子數目
-
blockSize 用逗號分割的外顯子的大小,這個item的數目對應於BlockCount的數目
-
blockStarts- 用逗號分割的列表, 所有外顯子的起始位置,數目也與blockCount數目對應.
3、一個簡單的示例如下:
track name=pairedReads description="Clone Paired Reads" useScore=1
chr22 1000 5000 cloneA 960 + 1000 5000 0 2 567,488, 0,3512
chr22 2000 6000 cloneB 900 - 2000 6000 0 2 433,399, 0,3601
bed格式有相應的軟體來處理這類格式的檔案,如bedtools。
-
注意:用於在GBrowse上展示相關注釋的bed格式通常第一行有一個關於track的描述資訊。
速記:
-
bed不是床,缺了主要3列就得黃~
-
9列可選列,看了不會胡略略~
二、BED 與GFF的差異
BED檔案中起始座標為0,結束座標至少是1,;
GFF中起始座標是1而結束座標至少是1。
參考連結:
http://5527lok.blog.163.com/blog/static/64751582011530113134590/ http://blog.sciencenet.cn/home.php?mod=space&uid=442719&do=blog&id=721452 https://www.plob.org/article/3748.htmlhttp://ju.outofmemory.cn/entry/193943 https://en.wikipedia.org/wiki/General_feature_format
MAF
簡介
MAF格式通常用於記錄體細胞突變。 MAF本來並不是一個常見的文字檔案格式,只是因為癌症研究實在是太熱門了,對它的理解也變得需求旺盛起來了。
一、MAF的說明
這些檔案應該使用下面描述的突變註釋格式(MAF)進行格式化。另外下文中有檔案命名規範。
以下幾種型別的體細胞突變會在MAF檔案中出現:
- 錯義突變及無義突變 *剪接位點,其定義為剪接位點2 bp以內的SNP *沉默突變 *與基因的編碼區、剪接位點或遺傳元件目標區域重疊的引物。 *移碼突變 *調控區突變
大部分MAF提交提交的是原始資料。這些原始資料中在體細胞中標記的位點與已知的變異型別相重合的。為避免有可能出現的細胞系汙染,MAF規定了一定的下細胞過濾標準。根據現行政策,可開放獲取MAF資料應滿足:
- 包括所有已驗證的體細胞突變名稱 *包括與編碼區域或剪接位點重疊的所有未驗證的體細胞突變名稱 *排除所有其他型別的突變(即非體細胞突變、不在編碼區域或剪接位點的未驗證體的細胞突變以及未在dbSNP、COSMIC或OMIM中註釋為體細胞的dbSNP位點)
我們提交給DCC MAF存檔的資料包括兩種:Somatic MAF(named .somatic.maf)的開放訪問資料以及不經過篩選的包含原始資料的Protected MAF(named.protected.maf)。所有資料將使用MAF標準進行格式化。
二、MAF檔案欄位
MAF檔案的格式是以製表符分隔的列。這些列都在表1中進行了詳細註釋,每個MAF檔案中都必須按照相應格式進行處理,DCC將驗證每列的順序是否符合標註。每列的標題和值有時需要區分大小寫。列中允許出現空值(即空白單元格)或列舉值。驗證器將查詢是否存在一個符合規範(例如#version 2.4)的標題,如果沒有,驗證將會失敗。除了出現在表1中的列外,也可以選擇附加其他列。可選列不經過DCC驗證,可以按任何順序進行。
MAF檔案可能有兩種格式 ,可能是47列,或者120列,第一行一般都是 標頭檔案,註釋著每一列的資訊,的確,資訊量有點略大。如下:
1 Hugo_Symbol
2 Entrez_Gene_Id
3 Center
4 NCBI_Build
5 Chromosome
6 Start_Position
7 End_Position
8 Strand
9 Consequence
10 Variant_Classification
11 Variant_Type
12 Reference_Allele
13 Tumor_Seq_Allele1
14 Tumor_Seq_Allele2
15 dbSNP_RS
16 dbSNP_Val_Status
17 Tumor_Sample_Barcode
18 Matched_Norm_Sample_Barcode
19 Match_Norm_Seq_Allele1
20 Match_Norm_Seq_Allele2
21 Tumor_Validation_Allele1
22 Tumor_Validation_Allele2
23 Match_Norm_Validation_Allele1
24 Match_Norm_Validation_Allele2
25 Verification_Status
26 Validation_Status
27 Mutation_Status
28 Sequencing_Phase
29 Sequence_Source
30 Validation_Method
31 Score
32 BAM_File
33 Sequencer
34 t_ref_count
35 t_alt_count
36 n_ref_count
37 n_alt_count
38 HGVSc
39 HGVSp
40 HGVSp_Short
41 Transcript_ID
42 RefSeq
43 Protein_position
44 Codons
45 Hotspot
46 cDNA_change
47 Amino_Acid_Change
1 Hugo_Symbol
2 Entrez_Gene_Id
3 Center
4 NCBI_Build
5 Chromosome
6 Start_Position
7 End_Position
8 Strand
9 Variant_Classification
10 Variant_Type
11 Reference_Allele
12 Tumor_Seq_Allele1
13 Tumor_Seq_Allele2
14 dbSNP_RS
15 dbSNP_Val_Status
16 Tumor_Sample_Barcode
17 Matched_Norm_Sample_Barcode
18 Match_Norm_Seq_Allele1
19 Match_Norm_Seq_Allele2
20 Tumor_Validation_Allele1
21 Tumor_Validation_Allele2
22 Match_Norm_Validation_Allele1
23 Match_Norm_Validation_Allele2
24 Verification_Status
25 Validation_Status
26 Mutation_Status
27 Sequencing_Phase
28 Sequence_Source
29 Validation_Method
30 Score
31 BAM_File
32 Sequencer
33 Tumor_Sample_UUID
34 Matched_Norm_Sample_UUID
35 HGVSc
36 HGVSp
37 HGVSp_Short
38 Transcript_ID
39 Exon_Number
40 t_depth
41 t_ref_count
42 t_alt_count
43 n_depth
44 n_ref_count
45 n_alt_count
46 all_effects
47 Allele
48 Gene
49 Feature
50 Feature_type
51 One_Consequence
52 Consequence
53 cDNA_position
54 CDS_position
55 Protein_position
56 Amino_acids
57 Codons
58 Existing_variation
59 ALLELE_NUM
60 DISTANCE
61 TRANSCRIPT_STRAND
62 SYMBOL
63 SYMBOL_SOURCE
64 HGNC_ID
65 BIOTYPE
66 CANONICAL
67 CCDS
68 ENSP
69 SWISSPROT
70 TREMBL
71 UNIPARC
72 RefSeq
73 SIFT
74 PolyPhen
75 EXON
76 INTRON
77 DOMAINS
78 GMAF
79 AFR_MAF
80 AMR_MAF
81 ASN_MAF
82 EAS_MAF
83 EUR_MAF
84 SAS_MAF
85 AA_MAF
86 EA_MAF
87 CLIN_SIG
88 SOMATIC
89 PUBMED
90 MOTIF_NAME
91 MOTIF_POS
92 HIGH_INF_POS
93 MOTIF_SCORE_CHANGE
94 IMPACT
95 PICK
96 VARIANT_CLASS
97 TSL
98 HGVS_OFFSET
99 PHENO
100 MINIMISED
101 ExAC_AF
102 ExAC_AF_Adj
103 ExAC_AF_AFR
104 ExAC_AF_AMR
105 ExAC_AF_EAS
106 ExAC_AF_FIN
107 ExAC_AF_NFE
108 ExAC_AF_OTH
109 ExAC_AF_SAS
110 GENE_PHENO
111 FILTER
112 CONTEXT
113 src_vcf_id
114 tumor_bam_uuid
115 normal_bam_uuid
116 case_id
117 GDC_FILTER
118 COSMIC
119 MC3_Overlap
120 GDC_Validation_Status
重要標準
列表中每列的順序最好與索引列相同。
標有“Case Sensitive“的列所有標題都需要區分大小寫。
標有”Null“的列表示允許具有空值。
標有“Enumerated column”的列表示具有指定的值,比如“Enumerated value” 是"No"表示該列沒有指定的值;其他值表示允許列出的具體值; “Set”表示該列的值來自指定的一組已知值(例如HUGO基因符號)
三、MAF檔案檢查
DCC 檔案驗證器將檢查MAF檔案的完整性。如果MAF檔案中的任何一項出現錯誤,驗證將失敗:
1. 列標題文字(包括大小寫)和順序必須與表1完全一致
2. 表1中列出的列標題下的值不是空值時必須具有相應的值
3. 表1中指定為“Case Sensitive”的值必須區分大小寫。
4. 如果列標題在規範中列出為具有列舉值(即“Enumerated”列的“Yes”),則這些列中的值必須來自“Enumerated”下列出的列舉值。
5. 如果規範中列標題具有設定值(即“Enumerated”列的“Set”),則那些列下的值必須來自該域的列舉值(例如,HUGO基因符號)。
6. 所有Allele-based列必須包含 –(deletion)或由以下大寫字母組成的字串:A,T,G,C。
7. 如果Validation_Status == "Untested" 那麼Tumor_Validation_Allele1, Tumor_Validation_Allele2, Match_Norm_Validation_Allele1, Match_Norm_Validation_Allele2 可以是空值(由 Validation_Status決定).
a) 如果Validation_Status == "Inconclusive" 那麼Tumor_Validation_Allele1, Tumor_Validation_Allele2, Match_Norm_Validation_Allele1, Match_Norm_Validation_Allele2 可以是空值(由 Validation_Status決定)
8. 如果Validation_Status == Valid, 那麼Validated_Tumor_Allele1 and Validated_Tumor_Allele2一定要是 A, C, G, T, 或“-“中的一種
a) 如果Validation_Status == "Valid" 那麼Tumor_Validation_Allele1, Tumor_Validation_Allele2, Match_Norm_Validation_Allele1, Match_Norm_Validation_Allele2 不可以是空值
b) 如果Validation_Status == "Invalid" 那麼Tumor_Validation_Allele1, Tumor_Validation_Allele2, Match_Norm_Validation_Allele1, Match_Norm_Validation_Allele2 不可以是空值Tumor_Validation_Allelle1 == Match_Norm_Validation_Allele1
Tumor_Validation_Allelle2 == Match_Norm_Validation_Allele2 (出現錯誤時,增加以替代8a)
9. 檢查Mutation_Status的等位基因值:
檢查Validation_status的等位基因值:
a) 如果Mutation_Status == "Germline" and Validation_Status == "Valid", 那麼Tumor_Validation_Allele1 == Match_Norm_Validation_Allele1 Tumor_Validation_Allele2 == Match_Norm_Validation_Allele2.
b) 如果Mutation_Status == "Somatic" ,Validation_Status == "Valid", 那麼Match_Norm_Validation_Allele1 == Match_Norm_Validation_Allele2 == Reference_Allele 且(Tumor_Validation_Allele1 or Tumor_Validation_Allele2) != Reference_Allele
c) 如果Mutation_Status == "LOH" and Validation_Status=="Valid", 那麼Tumor_Validation_Allele1 == Tumor_Validation_Allele2 Match_Norm_Validation_Allele1 != Match_Norm_Validation_Allele2 and Tumor_Validation_Allele1 == (Match_Norm_Validation_Allele1 or Match_Norm_Validation_Allele2).
10. 檢查 Start_position <= End_position
11. 根據Variant_Type檢查Start_position和End_position:
a) 如果Variant_Type == "INS", 那麼(End_position - Start_position + 1 == length (Reference_Allele) 或End_position - Start_position == 1) 且length(Reference_Allele) <= length(Tumor_Seq_Allele1 and Tumor_Seq_Allele2)
b) 如果Variant_Type == "DEL", 那麼End_position - Start_position + 1 == length (Reference_Allele),
且length(Reference_Allele) >= length(Tumor_Seq_Allele1 and Tumor_Seq_Allele2)
c) 如果Variant_Type == "SNP", 那麼length(Reference_Allele and Tumor_Seq_Allele1 and Tumor_Seq_Allele2) == 1 且(Reference_Allele and Tumor_Seq_Allele1 and Tumor_Seq_Allele2) != "-"
d) 如果Variant_Type == "DNP", 那麼length(Reference_Allele 和Tumor_Seq_Allele1 and Tumor_Seq_Allele2) == 2 且(Reference_Allele and Tumor_Seq_Allele1 andTumor_Seq_Allele2) !contain "-"
e) 如果Variant_Type == "TNP", 那麼length(Reference_Allele and T