常⽤⽣物信息学格式介绍(fasta、fastq、gff2、
gtf(gff2.5)、gff3。。。
前⾔
在各个⾏业都是有⾏业标准的,这样才能统⼀规范⽽⽅便后⾯的分析,在⽣物信息学领域中主要是各种⼤量序列数据、注释数据等,这些都是有特定的格式去表⽰,下⾯列举⼏种常见的格式。了解这些是进⾏后续⽣物信息学分析的必备知识,有些⼈虽说是在做⽣物信息学分析,但是到现在可能还不知道什么是GFF3格式等。
fasta
>gi|129295|sp|P01013|OVAX_CHICK GENE X PROTEIN (OVALBUMIN-RELATED) QIKDLLVSSSTDLDTTLVLVNAIYFKGMWKTAFNAEDTREMPFHVTKQESKPVQMMCMNNSFNVATLPAE KMKILELPFASGDLSMLVLLPDEVSDLERIEKTINFEKLTEWTNPNTMEKRRVKVYLPQMKIEEKYNLTS VLMALGMTDLFIPSANLTGISSAESLKISQAVHGAFMELSEDGIEMAGSTGVIEDIKHSPESEQFRADHP
FLFLIKHNPTNTIVYFGRYWSP
其中主要分为两个部分
第⼀部分是序列的定义⾏(单⾏),该⾏的开头是>符号,紧跟着后⾯的就是该条序列的名称(具有唯⼀性,即不能和其它序列同名称),即>号和后⾯的名称的第⼀字符间是没有任何空⽩的。⼀般第⼀个空格后⾯的内容即为可选的描述信息。如上⾯,
gi|129295|sp|P01013|OVAX_CHICK为序列名称, ⽽GENE X PROTEIN (OVALBUMIN-RELATED)则为描述信息。注意:有点软件是把⼀整⾏当做名称的,所以在出现错误的时候可以查看下格式是否正确。
第⼆部分就是序列,所有的序列碱基或者氨基酸可以都放在⼀⾏存储,也可以多⾏存储,但是建议⼤家多⾏存储且单⾏长度不超过80个字符,因为这样容易阅读。且序列的多⾏之间不能有空⾏,序列信息描述的第⼀⾏与序列数据的第⼀⾏之间不能有空⾏。其中序列数据主要是按照密码表来表⽰的,*表⽰是蛋⽩质翻译的结束。
多⾏序列举例如下:
>SEQUENCE_1
MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK IPQFASR
烟草网络KQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL
>SEQUENCE_2
SATVSEINSETDFVAKNDQFIALTKDTTAHIQSNSLQSVEELHSSTINGVKFEEYLKSQI ATIGENLVVRRFATLKAGANGVVNGYIHTNGRVGVVIAAACDSAEVASKSRDLLRQICMH
fastq
fastq格式⽂件中⼀个完整的单元分为四⾏,每⾏的含义如下:
第⼀⾏: 以@开头,内容同fasta的描述⾏类似
第⼆⾏:具体的碱基序列
第三⾏:以+开头,后⾯的内容可以和第⼀⾏类似,也什么都没有只留+
第四⾏:以ASCII字符集(分数)编码来表⽰对应碱基的测序质量
⽐如下⾯的这个例⼦:
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65
下⾯以Illumina和NCBI SRA两个测序数据来源来讲讲它们之间的区别:
小和尚念经歇后语
通常我们获取测序数据有两种途径,⼀种是⾃⼰通过仪器测定,⼀种是在公共数据库中(⽐如之前说到的NCBI中的SRA数据库)获取,这两种⽅式主要是在序列名称的命名上和测序质量表⽰⽅式上有所不同。
Illumina 序列名称:
@HWUSI-EAS100R:6:73:941:1973#0/1
上述以:隔开的每个字段的含义如下:
| HWUSI-EAS100R | the unique instrument name |
| 6 | flowcell lane |
| 73 | tile number within the flowcell lane |
| 941 | ‘x’-coordinate of the cluster within the tile |
| 1973 | ‘y’-coordinate of the cluster within the tile |
| #0 | index number for a multiplexed sample (0 for no indexing) |
| /1 | the member of a pair, /1 or /2 (paired-end or mate-pair reads only) |
NCBI SRA数据库:
将测序数据提交到NCBI的SRA数据库时,SRA数据库会为每⼀个样本提供⼀个编号,⼀般是SRRxxxxx,所以从SRA数据库上下载公共的测试数据(原始格式为
.sra, 需特定⼯具转换为fastq),其fastq格式⽂件中每个单元的名称是以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
需要注意的是:当把测序数据上传到SRA数据库时,它通常会将表⽰质量的分数 **转换为标准的Sanger格式 **。
质量分数表⽰法:
晋祠简介由于测序仪器的不同等因素所以对碱基测序质量的表⽰⽅式也不相同,在Fastq格式⽂件中,⽤ASCII码表来表⽰每个碱基的测序质量,下⾯介绍⼏种不同的⽅案:
image.png
其中有五种表⽰⽅法,Sanger的码表范围为!⾄I,其对应的数值为33-73,如果减去33(即Phred+33表⽰法)这个基数则范围转换为0-40,即如果某⼀个碱基的测序质量为!则对应的测序质量分数为0,表⽰测序质量低。其它⼏种表⽰法类似(X,I,J,L)。这⾥介绍测序质量的表⽰⽅法是因为后⾯有的软件是要指定测序数据的质量表⽰⽅法。
gff2
GFF(General Feature Format)是⼀种⽤于描述基因或者其它序列元素的⽂件格式,GFF有⼏个版本,早期的第Version 2和现在的Version 3. Version 2 是由Sanger机构所制定的,⽽Version 3是由Sequence Ontology Project制定。正是由于有统⼀的格式来表⽰基因等元素,使得GFF格式的⽂件被⼴泛的使⽤与mapping与基因组数据可视化⽅⾯。
GFF2⽂件格式是由tab隔开的九列值,每⼀⾏的九个字段的含义如下:
Chr1 curated CDS 365647 365963 . + 1 Transcript "R119.7"
第⼀列: reference quence, 该列表⽰的是特征元素所在的染⾊体(或者scaffold,或者contig),也就是在基因组中的坐标系统,后续⼀切的注释信息都是基于此列。
第⼆列:source,该列表⽰改⾏注释信息的来源,⽐如上述的⼀⾏表⽰该⾏的CDS注释信息来⾃名为“curated”的注释。
第三列:feature,或者说是method,type, 表⽰的是该注释的类型,⽐如上述表⽰改⾏注释为CDS信息,可以将source和feature结合起来描述的更加详细。
第四列:start position,在reference quence上的开始位置(坐标),通常是从1为起点⽽不是0。
第五列:end position, 在reference quence上的结束位置(坐标),⼀般是⼤于start position的。
第六列:score, 表⽰该⾏feature的分数,⽐如序列相似性等,如果没有对应的分数可以⽤.代替。
第七列:strand,feature所在链,+表⽰正链,-表⽰负链,.表⽰不确定或者与链⽆关。
第⼋列:pha,与蛋⽩质编码相关,⼀般是⽤于CDS,值的范围为0-2,表⽰编码时阅读框的移动相位。
下⾯这段描述很详细:
‘0’ indicates that the specified region is in frame, i.e. that its first ba corresponds to the first ba of
a codon. ‘1’indicates that there is one extra ba, i.e. that the cond ba of the region corresponds to the first ba of a codon, and ‘2’ means that the third ba of the region is the first ba of a codon. If the strand is ‘-‘, then the first ba of the region is value of <end>, becau the corresponding coding region will run from <end> to <start> on the rever strand.
第九列:group,或者称为attributes,是⽤于对改⾏注释更多的描述,以键值对的形式,⽐如上⾯的例⼦表⽰该CDS是属于名为R119.7的transcript。该列中可以存在多个属性,属性之间是⽤;隔开的。
对于GFF格式的理解主要是集中在最后⼀列,有以下集中情况:
1. 对于单个feature
Chr3 giemsa heterochromatin 4500000 6000000 . . . Band 3q12.1
2. 对于属于同⼀集合的多个feature
IV curated exon 5506900 5506996 . + . Transcript B0273.1
IV curated exon 5506026 5506382 . + . Transcript B0273.1
IV curated exon 5506558 5506660 . + . Transcript B0273.1
IV curated exon 5506738 5506852 . + . Transcript B0273.1
⽐如上⾯这个例⼦就表⽰这四个exon都是属于同⼀个名为B0273.1的transcript,这是表⽰⼀个完整transcript结构的最基本要求。
GFF2还可⽤于序列⽐对结果表⽰等其它⽅⾯,这⾥不做介绍了。
gtf(gff2.5)
GTF(Gene Transfer Format)格式是借鉴于GFF2格式,也被称为GFF2.5,⼤部分字段的定义是和GFF2相同的,只是每⾏的第九列必须带有如下四个域,具体为gene_id value; transcript_id value; 这样的设计是为了适应⼀个基因的多个转录本这种情况。⽐如下⾯的这个例⼦:
AB000123 Twinscan CDS 193817 194022 . - 2 gene_id "AB000123.1"; transcript_id "AB00123.1.2";
AB000123 Twinscan CDS 199645 199752 . - 2 gene_id "AB000123.1"; transcript_id "AB00123.1.2";
AB000123 Twinscan CDS 200369 200508 . - 1 gene_id "AB000123.1"; transcript_id "AB00123.1.2";
陈坤志AB000123 Twinscan CDS 215991 216028 . - 0 gene_id "AB000123.1"; transcript_id "AB00123.1.2";
AB000123 Twinscan start_codon 216026 216028 . - . gene_id "AB000123.1"; transcript_id "AB00123.1.2";
AB000123 Twinscan stop_codon 193814 193816 . - . gene_id "AB000123.1"; transcript_id "AB00123.1.2";
gff3
GFF2格式早期⽤的⽐较多,但是现在⽤的多的是GFF3格式,这也是好多软件所⽀持的,⽐如Gbrow, Jbrow等基因组数据可视化⼯具。
先看下⾯这个简单的例⼦:
##gff-version 3
ctg123 . exon 1300 1500 . + . ID=exon00001
新生儿溶血病ctg123 . exon 1050 1500 . + . ID=exon00002
ctg123 . exon 3000 3902 . + . ID=exon00003
ctg123 . exon 5000 5500 . + . ID=exon00004
ctg123 . exon 7000 9000 . + . ID=exon00005
第⼀⾏的##gff-version 3通常是需要的,⽽且必须是在⽂件的第⼀⾏。
前⼋列和GFF2、GFF2.5类似,但是有⼏点是要特别注意的,主要是将GFF3注释数据⽤于基因组浏览器时,字段中的⼀些特殊字符⽐如空格,> %等都需要使⽤URL编码进⾏转换才能准确的在web中进⾏展⽰。
第九列同样是表⽰attributes,采⽤的同样是键值对的形式(tag=value),只是这⾥有⼏个特定的键,具体如下:
ID,feature在整个GFF3⽂件中唯⼀的标识符;
Name,feature的名字,不同于ID,Name不要求唯⼀,只是⽅便⽤户浏览;
Alias, 相当于feature的别名;
Parent,表明该feature所属的上⼀级feature 的ID,这种关系可⽤于exons-transcripts,transcripts-genes,可以看出⼀个feature可以拥有多个⼦feature;
Target, 主要是⽤于序列⽐对结果的展⽰,value的格式为target_id start end [strand], 其中如果target_id中含有空格则需转换
为%20;
后⾯还有些其它属性⽐如Note等,这⾥不再做详细描述。
下⾯再来看下典型的例⼦:
蛋⽩质编码基因结构
ctg123 example gene 1050 9000 . + . ID=EDEN;Name=EDEN;Note=protein kina
ctg123 example mRNA 1050 9000 . + . ID=EDEN.1;Parent=EDEN;Name=EDEN.1;Index=1
ctg123 example five_prime_UTR 1050 1200 . + . Parent=EDEN.1
ctg123 example CDS 1201 1500 . + 0 Parent=EDEN.1
ctg123 example CDS 3000 3902 . + 0 Parent=EDEN.1
ctg123 example CDS 5000 5500 . + 0 Parent=EDEN.1
ctg123 example CDS 7000 7608 . + 0 Parent=EDEN.1
ctg123 example three_prime_UTR 7609 9000 . + . Parent=EDEN.1
ctg123 example mRNA 1050 9000 . + . ID=EDEN.2;Parent=EDEN;Name=EDEN.2;Index=1
ctg123 example five_prime_UTR 1050 1200 . + . Parent=EDEN.2
ctg123 example CDS 1201 1500 . + 0 Parent=EDEN.2
ctg123 example CDS 5000 5500 . + 0 Parent=EDEN.2
ctg123 example CDS 7000 7608 . + 0 Parent=EDEN.2
ctg123 example three_prime_UTR 7609 9000 . + . Parent=EDEN.2
ctg123 example mRNA 1300 9000 . + . ID=EDEN.3;Parent=EDEN;Name=EDEN.3;Index=1
怪你过份美丽
ctg123 example five_prime_UTR 1300 1500 . + . Parent=EDEN.3
ctg123 example five_prime_UTR 3000 3300 . + . Parent=EDEN.3
ctg123 example CDS 3301 3902 . + 0 Parent=EDEN.3
ctg123 example CDS 5000 5500 . + 1 Parent=EDEN.3
ctg123 example CDS 7000 7600 . + 1 Parent=EDEN.3
ctg123 example three_prime_UTR 7601 9000 . + . Parent=EDEN.3
⼀个名为EDEN的基因拥有三个转录本,分别名为EDEN.1 EDEN.2 EDEN.3, 每个转录本⼜有UTR和CDS等信息。
序列⽐对
ctg123 est EST_match 1050 1500 . + . ID=Match1;Name=agt830.5;Target=agt830.5 1 451
ctg123 est EST_match 3000 3202 . + . ID=Match1;Name=agt830.5;Target=agt830.5 452 654
ctg123 est EST_match 5410 5500 . - . ID=Match2;Name=agt830.3;Target=agt830.3 505 595
ctg123 est EST_match 7000 7503 . - . ID=Match2;Name=agt830.3;Target=agt830.3 1 504
ctg123 est EST_match 1050 1500 . + . ID=Match3;Name=agt221.5;Target=agt221.5 1 451
ctg123 est EST_match 5000 5500 . + . ID=Match3;Name=agt221.5;Target=agt221.5 452 952
兔子怎么喂在最美的年华里遇见你是什么歌ctg123 est EST_match 7000 7300 . + . ID=Match3;Name=agt221.5;Target=agt221.5 953 1253