生信中常见的数据文件格式

DoubleHelix

发布于 2020-07-09 15:33:08

2.5K0

发布于 2020-07-09 15:33:08

前面我们介绍了各种测序技术的原理：illumina、Sanger、第三代和第四代测序技术原理，我们测序得到的是带有质量值的碱基序列fastq格式，参考基因组是fasta格式。⽤⽐对⼯具把fastq格式的序列回帖到对应的fasta格式的参考基因组序列，就可以产⽣sam格式的⽐对⽂件。把sam格式的⽂本⽂件压缩成⼆进制bam⽂件可以节省空间。如果是记录某些位点或者区域碱基的变化，就是VCF⽂件格式。如果对参考基因组上⾯的各个区段标记它们的性质，⽐如哪些区域是外显⼦，内含⼦， UTR等等，这就是gtf/gff格式。如果只是为了单纯描述某个基因组区域，就是bed格式⽂件，记录染⾊体号以及起始终⽌坐标，正负链即可。

1.fastq文件

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

FASTQ文件中每个序列通常有四行：

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

第二行是序列

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

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

例如：

质量评分指的是一个碱基的错误概率的对数值。其最初在Phred拼接软件中定义与使用，对于每个碱基的质量编码标示，不同的软件采用不同的方案，目前有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+ ，Phred quality score。

举个犊子，假设质量信息是5，这个5看成字符，不要看成数字，那么他对应的十进制数就是53，如果使用的是Sanger测序，那么质量评分采用的就是Phred+33的形式，那么Q = 53-33=20=-10lg(errorP)，所以errorP = 0.01。也就计算出错误率啦，就便于我们进行质控。每一个碱基都有一个质量评分，所以第2行和第4行的位数是相同的。

2.fasta文件

FASTA格式是一种用于表示核苷酸序列或多肽序列的文本格式。其中碱基对或氨基酸用单个字母来表示，且允许在序列前添加序列名及注释。该格式已成为生物信息学领域的一项标准。

FASTA文件各行记录信息如下：

第一行是由大于号">"开头的任意文字说明，用于序列标记，为了保证后续分析软件能够区分每条序列，单个序列的标识必须是唯一的。

从第二行开始为序列本身，只允许使用既定的核苷酸或氨基酸编码符号。通常核苷酸符号大小写均可，而氨基酸常用大写字母。注意有些程序对大小写有明确要求。一般每行60～80个字母。

核苷酸序列：

氨基酸序列：

fasta格式还是比较常见的，比如我们在NCBI查看基因的的时候通常就有fasta格式genebank格式。下面就是fasta格式的案例：

3.SAM/BAM

当我们测序得到的fastq数据map到基因组之后，会得到一个以sam或bam为扩展名的文件。这里，SAM的全称是sequence alignment/map format。而BAM就是SAM的二进制文件，也就是压缩格式的sam文件。

SAM格式文件包括头部注释部分和比对结果部分，头部分为’’可选部分’’。头部分位于比对部分之前，以“@”开头。比对部分有11列是固定的，其他多列可选。

（1）Header （标头注释部分）

@HD VN:1.0 SO:coordinate 
@SQ SN:chr1 LN:249250621 
@SQ SN:chr10 LN:135534747 
@SQ SN:chr11 LN:135006516
 ... 
@SQ SN:chrY LN:59373566 
@PG ID:TopHat VN:2.0.8b CL:/home/hpages/tophat-2.0.8b.Linux_x86_64/tophat --mate-inner-dist 150 --solexa-quals --max-multihits 5 --no-discordant --no-mixed --coverage-search --microexon-search --library-type fr-unstranded --num-threads 2 --output-dir tophat2_out/ERR127306 /home/hpages/bowtie2-2.1.0/indexes/hg19 fastq/ERR127306_1.fastq fastq/ERR127306_2.fastq

@HD：说明VN的版本以及比对有无排列顺序，这个例子没有排序。
@SQ：参考序列目录。SN：参考序列名字。LN：参考序列长度。
@PG：使用的比对程序名。

（2）比对结果部分

例如这样的：

E00514:173:H3C3JCCXY:4:1124:12398:67234    337 Chr00   32904   0   150M    Chr09   33498107    0   TCAATTTCACTTGAAGCTTACTTGTAGTTTCAGGCTTGGTCAAGCGCGATACAAACCATGTAGTAGGAGTCCTCCAAGTCGCCAAGCTAGGGGATCTGCTGAAAGAGGTGACAGACAAGGTAAGCAATCAGAGCTCTAAGCAATCAGTCC  iieiiiii`eiiiiiiiiiiiiiiieiiiiiiiieiiiiiiiiiiiiiiiiiiiiieiiiiiiiiiiiiieiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiieiiiiiiiii`iiieeieieeieee``  AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:136C13 YT:Z:UU NH:i:8  CC:Z:Chr10  CP:i:18604313   HI:i:0  RG:Z:J36CK1
E00514:173:H3C3JCCXY:4:1124:12398:67234    369 Chr00   32904   0   150M    Chr16   2469225 0   TCAATTTCACTTGAAGCTTACTTGTAGTTTCAGGCTTGGTCAAGCGCGATACAAACCATGTAGTAGGAGTCCTCCAAGTCGCCAAGCTAGGGGATCTGCTGAAAGAGGTGACAGACAAGGTAAGCAATCAGAGCTCTAAGCAATCAGTCC  iieiiiii`eiiiiiiiiiiiiiiieiiiiiiiieiiiiiiiiiiiiiiiiiiiiieiiiiiiiiiiiiieiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiieiiiiiiiii`iiieeieieeieee``  AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:136C13 YT:Z:UU NH:i:8  CC:Z:Chr10  CP:i:18604313   HI:i:2  RG:Z:J36CK1
E00514:173:H3C3JCCXY:4:1124:12398:67234    369 Chr00   32904   0   150M    Chr16   29515410    0   TCAATTTCACTTGAAGCTTACTTGTAGTTTCAGGCTTGGTCAAGCGCGATACAAACCATGTAGTAGGAGTCCTCCAAGTCGCCAAGCTAGGGGATCTGCTGAAAGAGGTGACAGACAAGGTAAGCAATCAGAGCTCTAAGCAATCAGTCC  iieiiiii`eiiiiiiiiiiiiiiieiiiiiiiieiiiiiiiiiiiiiiiiiiiiieiiiiiiiiiiiiieiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiieiiiiiiiii`iiieeieieeieee``  AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:136C13 YT:Z:UU NH:i:8  CC:Z:Chr10  CP:i:18604313   HI:i:4  RG:Z:J36CK1
E00514:173:H3C3JCCXY:4:1124:12398:67234    369 Chr00   32904   0   150M    Chr17   31040767    0   TCAATTTCACTTGAAGCTTACTTGTAGTTTCAGGCTTGGTCAAGCGCGATACAAACCATGTAGTAGGAGTCCTCCAAGTCGCCAAGCTAGGGGATCTGCTGAAAGAGGTGACAGACAAGGTAAGCAATCAGAGCTCTAAGCAATCAGTCC  iieiiiii`eiiiiiiiiiiiiiiieiiiiiiiieiiiiiiiiiiiiiiiiiiiiieiiiiiiiiiiiiieiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiieiiiiiiiii`iiieeieieeieee``  AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:136C13 YT:Z:UU NH:i:8  CC:Z:Chr10  CP:i:18604313   HI:i:6  RG:Z:J36CK1
E00514:173:H3C3JCCXY:4:1212:19025:24532    409 Chr00   33538   0   150M    *   0   0   GATTCCAAGTGCTGACTGATTGCTCTCTTTCTCCTTGTCTTGCAGGTAAGAACAAGGCCAAAGGAAAAGACAGGGAAAAAACATGAAATGAGATACTCTTGCTTTTAACCCTGATGATATGAGATATTCTTGCTCTAGTATAGCTTGTTT  ii`e`ei[iiiiiiiiiiiiie[ieeieieiiiiieiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiieee``  AS:i:-12    XN:i:0  XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:52T33T63   YT:Z:UU NH:i:20 CC:Z:Chr01  CP:i:11331871   HI:i:0  RG:Z:J36CK1

字段之间也就是列之间由Tab隔开，每一字段具体含义参考下图：

qname：测序的reads的名称；
Flag：提供了一个比对文件的信息，比对文件可以设置或取消flag，flag整数是多个flag的同时表示。表示比对的结果，由数字表示，不同的数值含义不同：
- 1 ：代表这个序列采用的是PE双端测序
- 2：代表这个序列和参考序列完全匹配，没有错配和插入缺失
- 4：代表这个序列没有mapping到参考序列上
- 8：代表这个序列的另一端序列没有比对到参考序列上，比如这条序列是R1,它对应的R2端序列没有比对到参考序列上
- 16：代表这个序列比对到参考序列的负链上
- 32 ：代表这个序列对应的另一端序列比对到参考序列的负链上
- 64 ：代表这个序列是R1端序列， read1;
- 128 : 代表这个序列是R2端序列，read2；
- 256：代表这个序列不是主要的比对，一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置，其他都是次要的
- 512：代表这个序列在QC时失败了，被过滤不掉了（# 这个标签不常用）
- 1024: 代表这个序列是PCR重复序列（#这个标签不常用）
- 2048: 代表这个序列是补充的比对（#这个标签具体什么意思，没搞清楚，但是不常用。

具体还可以参考官网：https://www.samformat.info/sam-format-flag

rname：map到参考基因组后的染色体/序列/ contig的名称；
contig：拼接软件基于reads之间的overlap区，拼接获得的序列称为Contig（重叠群）
strand：比对的链；
pos：比对的最左边部分的坐标；
mapq：比对的映射质量；
CIGAR：CIGAR字符串，一个数字与字母交替构成的字符串，标记了这段reads不同位置的match情况。

qwidth：读段的长度；
RNEXT：双末端测序中下一个reads比对的参考系列的名称，如果没有则用 " * " 表示，如果和前一个reads比对到同一个参考序列则用" = "表示；
PNEXT：下一个reads比对到参考序列上的位置，如果没有则用0表示；
TLEN：序列模板的长度；
seq：比对的实际顺序；
qual：比对的质量字符串(fasta文件中的质量得分)；
cigar中会包含数字，代表了特定match持续了多少nt；以及不同的字符，代表了不同的match情况。如36M表示它没有插入或删除。

由于sam格式的文件通常都非常大，所以为了节省存储空间而将sam转换为二进制格式以便于存储，也就是bam文件。sam/bam文件可以由特定的一些软件（比如samtools）来处理的，包括格式互转、排序、建立索引等操作。

4.GTF和GFF文件

GTF全称是Gene transfer format，用于储存基因结构信息，总共有9列。关于gtf文件，我在文章：生信中各种ID转换中就已经有所应用。我之前在TCGA数据库差异分析的文章中，也是通过gtf文件进行ID转换的。

GFF全称为general feature format，这种格式主要是用来注释基因组。描述了基因组上各种特征的区间信息，包括染色体，基因，转录本等。GFF文件本质上是一个\t分隔的，和gtf文件差不多，共9列的纯文本文件。

NC_000010.11    BestRefSeq%2CGnomon    gene    35126830    35212958    .    +    .    ID=gene27850;Dbxref=GeneID:1390,HGNC:HGNC:2352,MIM:123812;Name=CREM;description=cAMP responsive element modulator;gbkey=Gene;gene=CREM;gene_biotype=protein_coding;gene_synonym=CREM-2,hCREM-2,ICER
NC_000010.11    BestRefSeq    mRNA    35126841    35179847    .    +    .    ID=rna82191;Parent=gene27850;Dbxref=GeneID:1390,Genbank:NM_001881.3,HGNC:HGNC:2352,MIM:123812;Name=NM_001881.3;gbkey=mRNA;gene=CREM;product=cAMP responsive element modulator%2C transcript variant 2;transcript_id=NM_001881.3
NC_000010.11    BestRefSeq    exon    35126841    35127193    .    +    .    ID=id995818;Parent=rna82191;Dbxref=GeneID:1390,Genbank:NM_001881.3,HGNC:HGNC:2352,MIM:123812;gbkey=mRNA;gene=CREM;product=cAMP responsive element modulator%2C transcript variant 2;transcript_id=NM_001881.3
NC_000010.11    BestRefSeq    exon    35148368    35148491    .    +    .    ID=id995819;Parent=rna82191;Dbxref=GeneID:1390,Genbank:NM_001881.3,HGNC:HGNC:2352,MIM:123812;gbkey=mRNA;gene=CREM;product=cAMP responsive element modulator%2C transcript variant 2;transcript_id=NM_001881.3
NC_000010.11    BestRefSeq    exon    35178889    35178986    .    +    .    ID=id995820;Parent=rna82191;Dbxref=GeneID:1390,Genbank:NM_001881.3,HGNC:HGNC:2352,MIM:123812;gbkey=mRNA;gene=CREM;product=cAMP responsive element modulator%2C transcript variant 2;transcript_id=NM_001881.3
NC_000010.11    BestRefSeq    exon    35179134    35179847    .    +    .    ID=id995821;Parent=rna82191;Dbxref=GeneID:1390,Genbank:NM_001881.3,HGNC:HGNC:2352,MIM:123812;gbkey=mRNA;gene=CREM;product=cAMP responsive element modulator%2C transcript variant 2;transcript_id=NM_001881.3
NC_000010.11    BestRefSeq    CDS    35148372    35148491    .    +    0    ID=cds57086;Parent=rna82191;Dbxref=CCDS:CCDS7184.1,GeneID:1390,Genbank:NP_001872.3,HGNC:HGNC:2352,MIM:123812;Name=NP_001872.3;Note=isoform 2 is encoded by transcript variant 2;gbkey=CDS;gene=CREM;product=cAMP-responsive element modulator isoform 2;protein_id=NP_001872.3
NC_000010.11    BestRefSeq    CDS    35178889    35178986    .    +    0    ID=cds57086;Parent=rna82191;Dbxref=CCDS:CCDS7184.1,GeneID:1390,Genbank:NP_001872.3,HGNC:HGNC:2352,MIM:123812;Name=NP_001872.3;Note=isoform 2 is encoded by transcript variant 2;gbkey=CDS;gene=CREM;product=cAMP-responsive element modulator isoform 2;protein_id=NP_001872.3
NC_000010.11    BestRefSeq    CDS    35179134    35179329    .    +    1    ID=cds57086;Parent=rna82191;Dbxref=CCDS:CCDS7184.1,GeneID:1390,Genbank:NP_001872.3,HGNC:HGNC:2352,MIM:123812;Name=NP_001872.3;Note=isoform 2 is encoded by transcript variant 2;gbkey=CDS;gene=CREM;product=cAMP-responsive element modulator isoform 2;protein_id=NP_001872.3

第一列是seqid, 代表序列ID, 通常是染色体的ID, 每条染色体拥有一个唯一的ID。
第二列是source, 代表基因结构的来源，可以是数据库的名称，比如来自genebank数据库，也可以是软件的名称，比如用GeneScan软件预测得到，当然，也可以为空，用.点号填充。
第三列是type, 代表区间对应的特征类型，比如gene, exon等。
第四列是start, 代表区间的起始位置。
第四列是end, 代表区间的终止位置。
第六列是score, 软件提供了统计值，如果没有，就用.填充。
第七列是strand, 代表正负链的信息, +表示正链，-表示负链，?表示不清楚正负链的信息，当正负链信息没有意义时，可以用.填充。
第八列是phase，当描述的是CDS区间信息时，需要指定翻译时开始的位置，取值范围包括0,1,2。
第九列是attributes, 表示属性，每种属性采用key=value 的形式，多个属性之间用;分号分隔。

gtf与gff的比较

5.BED文件

BED文件每行至少包括chrom，chromStart，chromEnd三列必选;另外还可以添加额外的9列可选，这些列的顺序是固定的。

chr3    124792319       124792562       ENSG00000276626 RF00100 -
chr1    92700819        92700934        ENSG00000201317 RNU4-59P        -
chr14   100951856       100951933       ENSG00000200823 SNORD114-2      +
chr22   45200954        45201019        ENSG00000221598 MIR1249 -
chr1    161699506       161699607       ENSG00000199595 RF00019 +

必选的三列：

chrom - 染色体的名称（例如chr3，chrY，chr2_random）或支架（例如scaffold10671）。
chromStart- 染色体或支架中特征的起始位置，染色体中的第一个碱基编号为0。
chromEnd- 染色体或支架中特征的结束位置。所述 chromEnd碱没有包括在特征的显示。

例如，染色体的前100个碱基定义为chromStart = 0，chromEnd = 100，并跨越编号为0-99的碱基。

9个可选的BED字段：

name - 定义BED行的名称。当轨道打开到完全显示模式时，此标签显示在Genome浏览器窗口中BED行的左侧，或者在打包模式下直接显示在项目的左侧。
score - 得分在0到1000之间。如果此注释数据集的轨迹线useScore属性设置为1，则得分值将确定显示此要素的灰度级别（较高的数字=较深的灰色）。
此表显示 Genome Browser将BED分数值转换为灰色阴影：
strand - 定义strand。要么“.” （=无绞线）或“+”或“ - ”。
thickStart- 绘制特征的起始位置（例如，基因显示中的起始密码子）。当没有厚部分时，thickStart和thickEnd通常设置为chromStart位置。
thickEnd - 绘制特征的结束位置（例如基因显示中的终止密码子）。
itemRgb- R，G，B形式的RGB值（例如255,0,0）。如果轨道行 itemRgb属性设置为“On”，则此RBG值将确定此BED行中包含的数据的显示颜色。注意：建议使用此属性的简单颜色方案（八种颜色或更少颜色），以避免压倒Genome浏览器和Internet浏览器的颜色资源。
blockCount- BED行中的块（外显子）数。
blockSizes- 块大小的逗号分隔列表。此列表中的项目数应与blockCount相对应。
blockStarts - 以逗号分隔的块开始列表。应该相对于chromStart计算所有 blockStart位置。此列表中的项目数应与blockCount相对应。

BED文件与GFF文件的区别与联系：