1. 程式人生 > >生物資訊資料存放型別 {#filetype}

生物資訊資料存放型別 {#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種方案:

  1. Sanger,Phred quality score,值的範圍從0到92,對應的ASCII碼從33到126,但是對於測序資料(raw read data)質量得分通常小於60,序列拼接或者mapping可能用到更大的分數。
  2. Solexa/Illumina 1.0, Solexa/Illumina quality score,值的範圍從-5到63,對應的ASCII碼從59到126,對於測序資料,得分一般在-5到40之間。
  3. Illumina 1.3+,Phred quality score,值的範圍從0到62對應的ASCII碼從64到126,低於測序資料,得分在0到40之間;
  4. Illumina 1.5+,Phred quality score,但是0到2作為另外的標示。
  5. Illumina 1.8+

下面是更為直觀的表示:

wikipedia-fastq-format

來自於 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個欄位包括:

  1. QNAME,比對片段的(template)的編號;

  2. 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
  1. RNAME,參考序列的編號,如果註釋中對SQ-SN進行了定義,這裡必須和其保持一致,另外對於沒有mapping上的序列,這裡是’*‘;

  2. POS,比對上的位置,注意是從1開始計數,沒有比對上,此處為0;

  3. MAPQ,mappint的質量;

  4. 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(錯配,位置是一一對應的);
  1. RNEXT,下一個片段比對上的參考序列的編號,沒有另外的片段,這裡是’*‘,同一個片段,用’=‘;

  2. PNEXT,下一個片段比對上的位置,如果不可用,此處為0;

  3. TLEN,Template的長度,最左邊得為正,最右邊的為負,中間的不用定義正負,不分割槽段(single-segment)的比對上,或者不可用時,此處為0;

  4. SEQ,序列片段的序列資訊,如果不儲存此類資訊,此處為’*‘,注意CIGAR中M/I/S/=/X對應數字的和要等於序列長度;

  5. QUAL,序列的質量資訊,格式同FASTQ一樣。read質量的ASCII編碼。

  6. 可選欄位(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分為兩部分,註釋資訊和變異位點記錄資訊。

註釋資訊通常以#開頭,會描述該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格式。

注意事項:    針對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列是:

  1. chrom, 染色體或scafflold 的名字(eg chr3, chrY, chr2_random, scaffold0671 )

  2. chromStart 染色體或scaffold的起始位置,染色體第一個鹼基的位置是0

  3. chromEnd 染色體或scaffold的結束位置,染色體的末端位置沒有包含到顯示資訊裡面。例如,首先得100個鹼基的染色體定義為chromStart =0 . chromEnd=100, 鹼基的數目是0-99

2、9 個額外的可選列:

  1. name 指定BED行的名字,這個名字標籤會展示在基因組瀏覽器中的bed行的左側。

  2. score 0到1000的分值,如果在註釋資料的設定中將原始基線設定為1,那麼這個分值會決定現示灰度水平(數字越大,灰度越高)

  3. strand 定義鏈的方向,''+” 或者”-”

  4. thickStart 起始位置(The starting position at which the feature is drawn thickly)(例如,基因起始編碼位置)

  5. thickEnd 終止位置(The ending position at which the feature is drawn thickly)(例如:基因終止編碼位置) 

  6. itemRGB 是一個RGB值的形式, R, G, B (eg. 255, 0, 0), 如果itemRgb設定為'On”, 這個RBG值將決定資料的顯示的顏色。

  7. blockCount BED行中的block數目,也就是外顯子數目

  8. blockSize 用逗號分割的外顯子的大小,這個item的數目對應於BlockCount的數目

  9. 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