关于高通量数据格式

• 原始测序下机数据格式：fasta / fastq
• 比对(mapping)到参考基因组后数据格式：bam / sam
• 注释数据格式：gff / gtf / bed
• 变异数据格式：vcf

关于Fasta

　　Fasta格式也称为Pearson格式，是一种基于文本用于表示核苷酸序列或氨基酸序列的格式。在这种格式中碱基对或氨基酸用单个字母来编码，且允许在序列前添加序列名及注释。

Fasta格式

　　序列文件的第一行是由大于号">"或分号";"打头的任意文字说明（习惯常用">"作为起始），用于序列标记。从第二行开始为序列本身，只允许使用既定的核苷酸或氨基酸编码符号（参见下表）。通常核苷酸符号大小写均可，而氨基酸常用大写字母。使用时应注意有些程序对大小写有明确要求。文件每行的字母一般不应超过80个字符。

FASTA格式的氨基酸序列实例：

>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken

ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA DIDGDGQVNYEEFVQMMTAK*

关于Fastq

　　FASTQ是基于文本的，保存生物序列（通常是核酸序列）和其测序质量信息的标准格式。其序列以及质量信息都是使用一个ASCII字符标示，最初由Sanger开发，目的是将FASTA序列与质量数据放到一起，目前已经成为高通量测序结果的事实标准。

　　Fastq的格式，，FASTQ文件中每个序列通常有四行：

第一行，序列标识以及相关的描述信息，以‘@’开头；

第二行是序列；

第三行以‘+’开头，后面是序列标示符、描述信息，或者什么也不加；

第四行，是质量信息，和第二行的序列相对应，每一个碱基或者氨基酸都有一个质量评分，根据评分体系的不同，每个字符的含义表示的数字也不相同。

 

　　原始数据简单处理的工具：seqtk工具，得到一些基本的统计数据，如：

1. 将fastq 文件转换成fasta 文件：seqtk seq ‐A input.fastq > output.fasta
2. 转换为fasta文件并进行质量值 或者 N碱基过滤标记：seqtk seq ‐aQ64 ‐q20 in.fq > out.fa  或  seqtk seq ‐aQ64 ‐q20 ‐n N in.fq > out.fa
3. 从文件中提取部分序列（name.lst列出）：seqtk subseq in.fq name.lst > out.fq
4. 截取两端（测序质量一般较差）的序列：seqtk trimfq ‐b 5 ‐e 10 in.fa > out.fa
5. 提取一定区域内的序列（reg.bed给出区域）：seqtk subseq in.fa reg.bed > out.fa

Sam&bam文件

　　SAM的全称是sequence alignment/map format。SAM是一种序列比对格式标准， 由sanger制定，是以TAB为分割符的文本格式。主要应用于测序序列mapping到基因组上的结果表示，当然也可以表示任意的多重比对结果。文件后缀名.sam

　　当测序得到的fastq文件map到基因组之后，我们通常会得到一个sam或者bam为扩展名的文件。而BAM就是SAM的二进制文件(B取自binary)。SAM(Sequence Alignment/Map)格式是一种通用的比对格式，用来存储reads到参考序列的比对信息。

　　Bowtie2是现下最流行的短序列比对软件，SAM(Sequence Alignment/Map)格式是一种通用的比对格式，用来存储reads到参考序列的比对信息。

注释：以@开头的行，。行：除注释外，每一行是一个read。。列：

第一列：序列名字，read name，read的名字通常包括测序平台等信息；；eg.ILLUMINA-379DBF:1:1:3445:946#0/1

第二列：标记数字总和，sum of flags，为flag的总和（整数）,flag取值见备注(3)；；eg.16

　　1? 序列是一对序列中的一个

　　2? 比对结果是一个pair-end比对的末端

　　4? 没有找到位点

　　8? 这个序列是pair中的一个但是没有找到位点

　　16? 在这个比对上的位点，序列与参考序列反向互补

　　32? 这个序列在pair-end中的的mate序列与参考序列反响互补

　　64 序列是 mate 1

　　128 序列是 mate 2

　　假如说标记为以上列举出的数目，就可以直接推断出匹配的情况。假如说标记不是以上列举出的数字，比如说83=（64+16+2+1），就是这几种情况值和。

第三列：参考序列名字，reference sequence name，实际上就是比对到参考序列上的染色体号。若是无法比对，则是*；；eg.chr1

第四列：在参考序列上的位置，position，read比对到参考序列上，第一个碱基所在的位置。若是无法比对，则是0；；eg.36576599

第五列：map得分，Mapping quality，比对的质量分数，越高说明该read比对到参考基因组上的位置越唯一。；；eg.42。。bowtie2有时并不能完全确定一个短的序列来自与参考序列的那个位置，特别是对于那些比较简单的序列。但是bowtie2会给出一个值来显示出 这个段序列来自某个位点的概率值，这个值就是mapping qulity。Mapping qulity的计算方法是：Q=-10log10p，Q是一个非负值，p是这个序列不来自这个位点的估计值。

　　    假如说一条序列在某个参考序列上找到了两个位点，但是其中一个位点的Q明显大于另一个位点的Q值，这条序列来源于前一个位点的可能性就比较大。Q值的差距越大，这独特性越高。

第六列：CIGAR值，碱基匹配上的碱基数。match/mismatch、insertion、deletion 对应字母 M、I、D；；eg.36M 表示36个碱基在比对时完全匹配；比如3S6M1P1I4M，前三个碱基被剪切去除了，然后6个比对上了，然后打开了一个缺口，有一个碱基插入，最后是4个比对上了，是按照顺序的；

注：第七列到第九列是mate(备注1)的信息，若是单末端测序这几列均无意义。

第七列：MRNM(chr)，mate的参考序列名字(reference sequence name)，实际上就是mate比对到的染色体号，若是没有mate，则是*；；eg.*

第八列：mate position，mate比对到参考序列上的第一个碱基位置，若无mate,则为0；；eg.0

第九列：估计出的片段的长度，当mate 序列位于本序列上游时该值为负值。ISIZE，Inferred fragment size.详见Illumina中paired end sequencing 和 mate pair sequencing，是负数，推测应该是两条read之间的间隔(待查证)，若无mate则为0；；eg.0

第十列：Sequence，就是read的碱基序列，如果是比对到互补链上则对read进行了reverse completed；eg. CGTTTCTGTGGGTGATGGGCCTGAGGGGCGTTCTCN

第十一列：ASCII，read质量的ASCII编码。；；eg.PY[[YY_______________QQQQbILKIGEFGKB

例：GEPKO:00021:00039 4 * 0 0 * * 0 0TGTGACTGCTGTACCAAGATGTCAGGAACTACTGGTGAAAACACCGCAGCATGTCAAGATCACAGATTTTGGGCCGGCCAAACTGCCGGTGGCGGA >??CCCCCC???DDADA?@A@?CCC>???DCCAC=@;;;;,666066;;@A;;;<;6;;;;;;;;;;;;/;C9C>C>A7<<1<<</+/8<<+--)- ZP:B:f,0.00399726,0.00412064,0.000124653 ZA:i:202 ZG:i:339 ZB:i:30 ZC:B:i,339,338,1,0  

第十二列之后：Optional fields，以tab建分割。，，eg.AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:35T0 YT:Z:UU

• AS:i? 匹配的得分
• XS:i? 第二好的匹配的得分
• YS:i? mate 序列匹配的得分
• XN:i? 在参考序列上模糊碱基的个数
• XM:i? 错配的个数
• XO:i? gap open的个数
• XG:i? gap 延伸的个数
• NM:i? 经过编辑的序列
• YF:i? 说明为什么这个序列被过滤的字符串
• YT:Z
• MD:Z? 代表序列和参考序列错配的字符串

 

 gff格式

　　gff格式是Sanger研究所定义，是一种简单的、方便的对于DNA、RNA以及蛋白质序列的特征进行描述的一种数据格式，已经成为序列注释的通用格式，比如基因组的基因预测，许多软件都支持输入或者输出gff格式。存文本文件，由tab键隔开的9列组成，以下是各列的说明：

Column 1: “seqid”：序列的编号，编号的有效字符[a-zA-Z0-9.:^*$@!+_?-|]

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代替。键是区分大小写的，以大写字母开头的键是预先定义好的，在后面可能被其他注释信息所调用。

例如：

ctg123 . mRNA 1050 9000 . + . ID=mRNA00002;Parent=gene00001

ctg123 . mRNA 1300 9000 . + . ID=mRNA00003;Parent=gene00001

ctg123 . exon 1300 1500 . + . Parent=mRNA00003

 

GFF和GTF是大量用于存储注释信息的数据格式。经常可以看到交替使用这两种格式。然而，GFF（一般特征格式）实际上意味着可用于任何基因组特征，而GTF（基因转移格式）被严格用于基因。

而GTF格式: GTF文件的第9列同GFF文件不同，虽然同样是标签与值配对的情况，但标签与值之间以空格分开，且每个特征之后都要有分号；（包括最后一个特征），其内容必须包括gene_id和transcript_id。

目前两种文件可以方便的相互转化，比如:使用Cufflinks软件的 的gffread。

BED文件，

    BED 文件格式提供了一种灵活的方式来定义的数据行，以用来描述注释的信息。BED行有3个必须的列和9个额外可选的列。 每行的数据格式要求一致。

必须包含的3列：

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

2.chromStart 染色体或scaffold的起始位置，染色体第一个碱基的位置是0

3.chromEn 染色体或scaffold的结束位置，染色体的末端位置没有包含到显示信息里面。例如，首先得100个碱基的染色体定义为chromStart =0 . chromEnd=100, 碱基的数目是0-99

还有9 个额外的可选列:

4. name 指定BED行的名字，这个名字标签会展示在基因组浏览器中的bed行的左侧。

5. score 0到1000的分值，如果在注释数据的设定中将原始基线设置为１，那么这个分值会决定现示灰度水平（数字越大，灰度越高），下面的这个表格显示Genome Browser

6. strand 定义链的方向，''+” 或者”-”

7. thickStart 起始位置（The starting position at which the feature is drawn thickly）(例如，基因起始编码位置）

8. thickEnd 终止位置（The ending position at which the feature is drawn thickly）（例如：基因终止编码位置）

9. itemRGB 是一个RGB值的形式, R, G, B (eg. 255, 0, 0), 如果itemRgb设置为'On”, 这个RBG值将决定数据的显示的颜色。

10. blockCount BED行中的block数目，也就是外显子数目

11. blockSize 用逗号分割的外显子的大小, 这个item的数目对应于BlockCount的数目

12. blockStarts- 用逗号分割的列表, 所有外显子的起始位置，数目也与blockCount数目对应.

 例如：

chr7 127471196 127472363 Pos1 0 + 127471196 127472363 255,0,0 

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 

.vcf文件

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample

chr1 43814979 COSM27286 G A 78.2854 PASS

CHROM-(chromosome): 表示变异位点是在哪个contig 里call出来的，如果是人类全基因组的话那就是chr1…chr22，chrX,Y,M。

POS - position ：变异位点相对于参考基因组所在的位置，如果是INDEL（插入缺失），位置是INDEL的第一个碱基位置

ID - identifier:  variant的ID。如果call出来的SNP存在于dbSNP数据库里，就会显示相应的dbSNP里的rs编号；若没有，则用’.'表示其为一个novel variant。

REF - reference base(s): 参考碱基，染色体上面的碱基，必须是ATCGN中的一个，N表示不确定碱基，在这个变异位点处，参考基因组中所对应的碱基和研究对象基因组中所对应的碱基。

ALT - alternate base(s): 与参考序列比较发生突变的碱基

QUAL - quality: Phred格式(Phred_scaled)的质量值，表 示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为variant存在的概率; 。

FILTER - _filter status:  理想情况下，QUAL这个值应该是用所有的错误模型算出来的，这个值就可以代表正确的变异位点了，但是事实是做不到的。因此，还需要对原始变异位点做进一步的过滤。无论你用什么方法对变异位点进行过滤，过滤完了之后，在FILTER一栏都会留下过滤记录，如果是通过了过滤标准，那么这些通过标准的好的变异位点的FILTER一栏就会注释一个PASS，如果没有通过过滤，就会在FILTER这一栏提示除了PASS的其他信息。如果这一栏是一个“.”的话，就说明没有进行过任何过滤。。

INFO - additional information:  这一行是variant的详细信息，具体如下：

FORMAT 和 BC1-1-base.sorted.bam：这两行合起来提供了’ BC1-1-base′这个sample的基因型的信息。’ BC1-1-base′代表这该名称的样品，是由BAM文件中的@RG下的 SM 标签决定的。

 

GT：表示这个样本的基因型，对于一个二倍体生物，GT值表示的是这个样本在这个位点所携带的两个等位基因。0表示跟REF一样；1表示表示跟ALT一样；2表示第二个ALT。当只有一个ALT 等位基因的时候，0/0表示纯和且跟REF一致；0/1表示杂合，两个allele一个是ALT一个是REF；1/1表示纯和且都为ALT；样品的基因型（genotype）。或：两个数字中间用’/'分 开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele； 1 表示样品中variant的allele； 2表示有第二个variant的allele。因此： 0/0 表示sample中该位点为纯合的，和ref一致； 0/1 表示sample中该位点为杂合的，有ref和variant两个基因型； 1/1 表示sample中该位点为纯合的，和variant一致

GQ：表示GT的质量值。表示的意义同QUAL。

DP：覆盖到这个位点的总的reads数量，相当于这个位点的深度（并不是多有的reads数量，而是达到一定质量值要求的reads数）。

AD：对应两个以逗号隔开的值，这两个值分别表示覆盖到REF和ALT碱基的reads数，相当于支持REF和支持ALT的测序深度。

PL:对应3个以逗号隔开的值，这三个值分别表示该位点基因型是0/0，0/1，1/1的没经过先验的标准化Phred-scaled似然值（L）。如果转换成支持该基因型概率（P）的话，由于L=-10lgP，那么P=10^（-L/10），因此，当L值为0时，P=10^0=1。因此，这个值越小，支持概率就越大，也就是说是这个基因型的可能性越大。

举个例子说明一下：

chr1 899282 rs28548431  C  T  [CLIPPED]  GT:AD:DP:GQ:PL  0/1:1,3:4:25.92:103,0,26

在这个位点，GT=0/1，也就是说这个位点的基因型是C/T；GQ=25.92，质量值并不算太高，可能是因为cover到这个位点的reads数太少，DP=4，也就是说只有4条reads支持这个地方的变异；AD=1,3，也就是说支持REF的read有一条，支持ALT的有3条；在PL里，这个位点基因型的不确定性就表现的更突出了，0/1的PL值为0，虽然支持0/1的概率很高；但是1/1的PL值只有26，有10^(-2.6)=0.25%的可能性是1/1；但几乎不可能是0/0，因为支持0/0的概率只有10^(-10.3)=5*10^-11。

chr1  873762  .  T  G  [CLIPPED]  GT:AD:DP:GQ:PL  0/1:173,141:282:99:255,0,255

其中最后面两列是相对应的，每一个tag对应一个或者一组值，如：

chr1：873762，GT对应0/1；AD对应173,141；DP对应282；GQ对应99；PL对应255,0,255。

 

 

在SAM输出的结果中每一行都包括十二项通过Tab分隔，从左到右分别是：

1 序列的名字

2 概括出一个合适的标记，各个数字分别代表

• 1? 序列是一对序列中的一个
• 2? 比对结果是一个pair-end比对的末端
• 4? 没有找到位点
• 8? 这个序列是pair中的一个但是没有找到位点
• 16? 在这个比对上的位点，序列与参考序列反向互补
• 32? 这个序列在pair-end中的的mate序列与参考序列反响互补
• 64 序列是 mate 1
• 128 序列是 mate 2

假如说标记为以上列举出的数目，就可以直接推断出匹配的情况。假如说标记不是以上列举出的数字，比如说83=（64+16+2+1），就是这几种情况值和。

3? 参考序列的名字

4 在参考序列上的位置

5? mapping qulity?? 越高则位点越独特

bowtie2有时并不能完全确定一个短的序列来自与参考序列的那个位置，特别是对于那些比较简单的序列。但是bowtie2会给出一个值来显示出 这个段序列来自某个位点的概率值，这个值就是mapping qulity。Mapping qulity的计算方法是：Q=-10log10p，Q是一个非负值，p是这个序列不来自这个位点的估计值。

假如说一条序列在某个参考序列上找到了两个位点，但是其中一个位点的Q明显大于另一个位点的Q值，这条序列来源于前一个位点的可能性就比较大。Q值的差距越大，这独特性越高。

Q值的计算方法来自与SAM标准格式，请查看SAM总结。

6 代表比对结果的CIGAR字符串，如37M1D2M1I，这段字符的意思是37个匹配，1个参考序列上的删除，2个匹配，1个参考序列上的插入。M代表的是alignment match(可以是错配)

7? mate 序列所在参考序列的名称

8 mate 序列在参考序列上的位置

9? 估计出的片段的长度，当mate 序列位于本序列上游时该值为负值。

10 read的序列

11 ASCII码格式的序列质量

12 可选的区域

• AS:i? 匹配的得分
• XS:i? 第二好的匹配的得分
• YS:i? mate 序列匹配的得分
• XN:i? 在参考序列上模糊碱基的个数
• XM:i? 错配的个数
• XO:i? gap open的个数
• XG:i? gap 延伸的个数
• NM:i? 经过编辑的序列
• YF:i? 说明为什么这个序列被过滤的字符串
• YT:Z
• MD:Z? 代表序列和参考序列错配的字符串