关键词:测序技术、序列格式、数据质控
阅读用时:全文共6小节,约4800字,约9分钟
关键词:测序技术、序列格式、数据质控
作为一名生信老鸟,深知一入生信深似海。
俗话说,不了解数据格式的程序猿不是好科研工作者,今天我们就聊聊随着高通量测序技术的发展,不同测序原理、不同测序仪的的原始下机文件格式都有哪些,走过了怎样的发展演化史?
2005年454公司的 Genome Sequencer 20 System,2006年illumina公司的Solexa和2007年ABI公司的SOLiD标志着第二代测序技术的诞生。之后陆续在06年推出454 FLX系统、10年推出illumina HiSeq系统和同年推出的基于半导体芯片测序技术的Ion PGM已经陆续成为当前基因测序的中坚力量。
从一代测序到二代测序,这段路走了近三十年,然而从二代迈向三代,仅用了不到10年。层出不穷的三代小鲜肉中,表现最为亮眼的莫过于10年推出的PacBio RS系列和12年推出的Nanopore MinION。15年PacBio又推出了万众瞩目的Sequel,Nanopore也在17年推出了GridION,为更方便、快捷、准确、经济的测序提供了无限可能。
NGS先行者:454及其序列数据
说到二代测序,454公司可谓是不折不扣的吃螃蟹的人。早在2005年,454就推出了第一个基于焦磷酸测序原理的高通量基因组测序系统——Genome Sequencer 20 System,这也是NGS领域的第一台测序仪。随后,454被罗氏收购,并在2006年推出了名动测序界的GS FLX测序系统。
454的技术核心在于乳液PCR和焦磷酸测序法,该技术能获得较长的序列,平均读长可达400 bp;相对于illumina平台其主要缺点是无法准确测定单碱基重复序列的长度。
454 GS-FLX测序仪产生的序列文件格式为SFF (Standard Flowgram Format),是一种二进制文件。虽然目前有很多软件可以直接读取SFF格式文件进行分析,但更多的软件还是需要把SFF格式转换为Fna-Qual格式,这种格式把序列和与序列对应的碱基质量分别放在以Fna和Qual为后缀的文件里(文本文件)。
①.fna数据文件保存碱基序列
>HKSD5CR01D6P3Ilength=70 xy=1599_2828 region=1 run=R_2012_03_15_01_23_26_GGAGTAGCATGCGTGACGAATCGTAGTTCCGACCATAACGATGCCGACCTTTGACCA
②.qual数据文件保存碱基质量值
>HKSD5CR01D6P3Ilength=70 xy=1599_2828 region=1 run=R_2012_03_15_01_23_26_40 40 40 40 40 40 4040 40 40 40 40 40 40 40 40 40 40 40 40 39 39 39 40 40 40 34 34 34 34 40 30 3030 40 39 39 39 38 38 37 40 40 38 38 32 24 17 17 20 20 26 30 30 36 36 37
我们可以采用多种方法将SFF文件转为Fna-Qual格式。在这列举三种方法:
方法一:python中的Biopython模块
from Bio import SeqIO
SeqIO.convert("raw.sff","sff-trim", "trimmed.fasta", "fasta")
SeqIO.convert("raw.sff","sff-trim", "trimmed.qual", "qual")
方法二:QIIME软件
process_sff.py -i raw.sff
方法三:Mothur软件
sffinfo(sfftxt=raw.sff.txt. fasta=True, flow=True)
随着测序成本的不断降低,尽管读长偏长,但成本过于昂贵的454测序技术逐渐退出了科研界。然而,虽然454的生命科学测序业务早已关停,454带来的影响却仍不容磨灭。来自众多大型科研项目(例如人类微生物组计划)中的海量454序列信息,仍是当今研究中不可或缺的重要参考数据和引用来源。
渐行渐远:SOLiD及其序列数据
再来谈谈几乎快被遗忘的ABI公司。事实上,在454和illumina走入市场之前,ABI一直处理测序市场上傲视群雄的霸主地位。在454和illumina初亮相后,由于SOLiD系统采用的是连接法而非PCR反应,因此对于高GC含量的样本具有非常大的优势。然而,后续SOLiD系统通量难以提升,且读长短、成本高,也已慢慢退出历史舞台。
SOLiD技术,其主要核心是油包水PCR和基于“双碱基编码原理”连接测序法。该技术主要优势在于每个位置的碱基均被检测两次,测序准确性高达99.9999%,算是第二代测序技术中准确性最高的;但是该技术的读长仅为2X50 bp,后续序列拼接等比较复杂,且由于是双碱基对应一个荧光信号,一旦发生错误,容易在荧光解码阶段产生连锁错误。
SOLiD技术的原始数据格式是csfasta及与其匹配的qual,可以利用http://edison.cremag.org/resources/seq-analysis/tools/solid2fastq/solid2fastq.py脚本或者http://www.bbmriwiki.nl/svn/bwa_45_patched/solid2fastq.pl脚本将其转为csfastq(color space fastq),与fastq格式一致,四行记录一碱基序列,区别在于第二行不是碱基序列,而是color的编码,如下:
@SRR2967009.1 100_1000_1168_F3
T10011023211201220121202030102221012302121010131001
+
2@@@@>@?@@@@<@@//;@@/@9?@8@=@@@6;6@66;<@6@67?2?;/@
colorspace编码规则为:
csfastq文件就可作为fastqc的输入文件进行质控啦,并可利用多款软件进行后续分析,如SHRiMP、Sequel、BFAST和Bowtie等等。
NGS领跑者:illumina及其序列数据
illumina在NGS乃至整个测序市场中的霸主地位已毋庸置疑。引用illumina官方说法:“世界上90%以上的测序数据都由Illumina仪器产生”,不较真的话,这句话确实在某种程度上反应了illumina雄踞NGS市场的现状。因此,我们用较大篇幅为大家分享illumina下机数据的格式和处理方法。
illumina技术核心在于桥式PCR和边合成边测序(SBS)。illumina测序系统的碱基读取也是基于化学发光法来的,给每一个碱基加入荧光基团,通过拍照捕捉发光的碱基,如此就得到了DNA的原始序列信息。
也就是说,illumina原始读取的是图像数据文件,解析后才形成碱基序列文件。对应到序列文件,其测序源文件为BCL格式文件(per-cycle BCL basecall file),BCL是一种包含碱基信号和图块质量信息的二进制文件,在进入下游分析前,BCL文件会经由Casava碱基识别(Base Calling)转化为原始测序序列(Sequenced Reads),转化得到的序列我们称之为Raw Data或Raw Reads,格式为Fastq。
下面我们就来详细解读下Fastq文件:
Fastq是一种文件文件,如果Fastq文件大小在M级,一般可以直接用文本编辑器(如Notepad++)打开。打开之后,就可以看到类似下面表格中第二列的内容,每四行记录一段序列,这个序列也是有名字的,叫reads。如果采用的是illumina MiSeq PE250,则第二行碱基和第四行碱基质量值均为250个,即reads长度为250 bp。
Illumina Seqence identifier | @HWUSI-EAS100R:6:73:941:1973#0/1 |
Read bases | GATTTGGGGTTCAAAGCAGTATCGATCAAATAG。。。 |
+ | + |
Phred quality scores(ASCII) | !''*((((***+))%%%++)(%%%%).1***-+。。。 |
首行“illumina sequence identifier”记录了哪些信息?
illumina sequence identifier中包含了此次测序的信息,根据平台不同包含的信息有所区别,一般来说可能包含以下的部分元素:
@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos>:<UMI> <read>:<is filtered>:<control number>:<index>
我们再回头看上面的例子,就可以对应查看这条序列文件中的测序信息了:
@ | Each sequence identifier line starts with @ |
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) |
注:Run、Flowcell、Lane和tile的详细辨析,可关注我们的喜马拉雅FM栏目《一分钟听懂NGS基础概念》第七期--“Run、Flowcell与Lane”,关注我们的公号并在后台回复“run”即可获得文字版名词辨析。
碱基质量值的读取和评判
illumina平台采用了Phred quality score的方式(即Phred=-10lgP,其中P是碱基被测错的概率)来作为碱基测序准确性的指标。Fastq文件不是直接存储Phred quality score的,而是利用ASCII码来存储,它们之间的转化方式为:若一个碱基被测错的概率为0.001,则Phred=30,Phred+33=63,查ASCII表,63对应的字符是“?”,则该碱基对应的质量代表值就是“?”。其实Phred值与ASCII码之间的转换关系经过了多次改变,如下图:
碱基质量解析必备神器--Fastqc软件
一般来说,生信工作者们会结合Fastqc软件识别和解析拿到的Fastq文件的碱基质量记录类型。
Fastqc软件下载地址为:
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
帮助文档在:
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/
Fastqc是java语言编写的,在安装有java的服务器上直接就可使用。使用方式也很简单,命令行输入“fastqcseqfile1.fq seqfile2.fq”即可。输出结果包括12模块,如下表:
1 | Basic Statistics | 碱基质量使用的版本就在该部分的Encoding说明 |
2 | Per base sequence quality | 统计每个cycle的碱基Phred值的平均值 |
3 | Per tile sequence quality | 统计每个tile中碱基Pherd值的平均 |
4 | Per sequence quality scores | 统计每条reads碱基Phred值的均值 |
5 | Per base sequence content | 统计每个cycle中A/T/C/G四种碱基的比例 |
6 | Per sequence GC content | 统计每条reads GC含量的平均值 |
7 | Per base N content | 统计每个cycle中N碱基(测序仪不能辨别某位置碱基是什么时,用N来表示)的比例 |
8 | Sequence Length Distribution | 统计reads长度 |
9 | Sequence Duplication Levels | 统计序列完全一致的reads的频数 |
10 | Overrepresented sequences | 统计序列完全一致的reads的比例,若超过0.1%,则认为该序列over-represented |
11 | Adapter Content | 统计reads中Adapter的比例 |
12 | Kmer Content | K默认值为5,统计是否存在连续的5个碱基在reads中大量出现 |
具体示例可见:
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/good_sequence_short_fastqc.html
具体每个模块的说明见:http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/:
Fastqc是个很贴心的小软件,除了illumina数据外,还可以将454的数据和三代PacBio数据转换成Fastq格式后进行质控。
① Fastqc用于454数据质控,示例见:http://www.bioinformatics.babraham.ac.uk/projects/fastqc/454_SRR073599_fastqc.html。
②Fastqc用于PacBio数据质控,示例见:
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/pacbio_srr075104_fastqc.html。
后起新秀:Ion Torrent及其序列数据
Ion Torrent平台最初由Jonathan Rothberg于2007年开发,几经转手,2013年后归属Thermo Fisher,是他们家的掌中宝。目前有2010年发布的PGM,2012年发布Proton,2015年发布了Ion S5/XL。主要利用半导体芯片(无需荧光标记)技术,具有经济、快速、简单、规模可扩展的特点,适用于临床诊断和医学等使用场景,目前国内已有达安OEM版本DA8600。
Ion Torrent的核心是半导体芯片测序技术,该技术采集的信号为pH值的变化,其他第二代测序技术主要是通过拍照采集荧光信号。该技术的主要优点是快速,主要缺点与Roche公司的454技术相同,在判读单碱基重复序列时存在困难。
该技术的序列文件是以Bam文件格式存储的,Bam也是一种二进制格式。计算时,Bam文件可以利用SAMtools软件将其转为文本格式的SAM文件:
samtools view -h raw.bam >raw.sam。
Bam/Sam文件有11列必填内容,第一列对应Fastq文件中的Seqence identifier,第10列是碱基序列,第11列是序列相对应的碱基质量值。Fastqc软件可以输入Bam/Sam文件,对测序序列进行质控。如何将Bam/Sam文件转成Fastq格式呢?大家当然会想到直接提取第一列、第11列、第12列信息就可以了。但是通常Bam/Sam文件是碱基序列比对后的存储文件,因此可能一条reads会对应多条比对信息,双端测序比对结果文件中可能已丢失配对序列相关信息,因此建议大家使用BEDTools、Picard等已有软件将Bam文件转成Fastq格式。
方法一:BEDTools软件
samtools sort -naln.bam aln.sort
bedtoolsbam2fastq -i aln.sort.bam -fq Read1.fq.gz -fq2 Read2.fq.gz
方法二:Picard软件
java -jar SamToFastq.jar INPUT=aln.sam FASTQ=Read1.fq.gz SECOND_END_FASTQ=Read2.fq.gz INCLUDE_NON_PF_READS=True VALIDATION_STRINGENCY=SILENT
三代明星:PacBio及其序列数据
再来看看江湖呼声渐涨的三代测序技术。目前三代测序市场上,表现最为抢眼的莫过于以PacBio公司的SMRT和Oxford Nanopore Technologies为代表的纳米孔单分子测序技术。与前两代相比,三代测序最为核心的特点就是单分子测序,测序过程无需进行PCR扩增。
PacBio SMRT技术其实也应用了边合成边测序的思想,并以SMRT芯片为测序载体。测序时,不需要对目标DNA进行PCR扩增,而是直接在目标片段两端加上两个发卡结构的接头,形成一个连续的环状结构。也因此,PacBio系统在读长上显示了极大的优势。目前比较受市场热捧的三代测序是PacBio的RSⅡ和2015年推出的Sequel。
PacBio下机产生的序列文件以HDF5格式存储。可以采用h5dump命令来查看H5文件内容。
1. 查看碱基序列:
h5dump –d /PulseData/BaseCalls/Basecall raw.h5 > Basecall.info ,文件内容如下:
DATA {(0): 67, 71,67, 67, 65, 71, 67, 71, 65, 65, 84, 71, 71, 67, 84, 71, 67, (17): 71, 71, 71,71, 65, 65, 71, 67, 65, 71, 65, 65, 65, 84, 84, 65, 84, (34): 67, 67, 71, 84,65, 65, 65, 67, 84, 71, 84, 84, 71, 67, 84, 71, 67,
该文件采用的ASCII码的编码方式存储的碱基序列:A=> 65, C=>67, G=>71, T=>84。
2. 查看碱基质量值:
h5dump -d /PulseData/BaseCalls/QualityValue raw.h5 > Basecall.quality,
文件内容如下,其碱基质量值采用与illumina技术一致:
DATA {(0): 51, 44,42, 44, 24, 24, 51, 51, 51, 51, 50, 20, 20, 20, 50, 51, 51, (17): 48, 48, 48,47, 9, 9, 9, 51, 51, 46, 31, 31, 31, 31, 44, 51, 51, 30, (35): 30, 51, 51, 7,7, 7, 7, 51, 51, 44, 44, 44, 51, 51, 50, 27, 27, 26,
长到天际:Nanopore及其序列数据
Oxford Nanopore 公司2005年在英国牛津成立,其运用的纳米孔测序技术使得DNA链在一个单通道中就能够被解码和识别,而不需要将长链打断成小短链。由于实现了DNA聚合酶内在自身的延续性和反应速度,Nanopore读长更长速度更快;同时由于能直接检测每个碱基的特征性电流,因而能对修饰碱基进行测序,对于表观遗传学研究具有极高的价值;因此,这款长到天际的测序仪,非常有潜力横扫当前测序格局。
2014年春天推出U盘大小的便携式MinION测序仪,仪器售价仅需$1000,据官网报道最长Reads可长达960 Kb,2014年10月推出平板大小的台式测序仪PromethION,有48个flow cell,可以单独运行也可以并行,2017年推出桌面式GridION X5测序仪。
Nanopore目前还主要在测试和生产阶段,尚未大规模应用,其应用主要体现在微生物等小基因组生物上。推出至今,其最亮眼的表现莫过于2014年西非埃博拉病毒爆发,MinION以最快的速度破译病毒序列,名噪一时。随着独特的纳米孔技术的成熟和完善,未来在即时检测、太空应用、大众检测等方面会有很大的想象空间。
Nanopore测序得到的序列文件的格式基础也是HDF5(https://support.hdfgroup.org/HDF5/),下机产生后缀为Fast5的序列文档。Fast5文件可经由Poretools软件(http://poretools.readthedocs.io/en/latest/)转换为Fastq文件或Fasta,然后进行后续数据分析。
① 应用Poretools将fast5转换为fastq,示例见:
http://poretools.readthedocs.io/en/latest/content/examples.html#poretools-fastq
② 应用Poretools将fast5转换为fasta,示例见:
http://poretools.readthedocs.io/en/latest/content/examples.html#poretools-fasta
总结一下,在测序市场中,一代测序因其准确度高,仍作为突变检测、单菌鉴定等的金标准而存在;以illumina HiSeq和MiSeq为代表的二代测序势头强劲,主打低成本和高通量,2017新机型NovaSeq更宣称已将测序成本降至百美金;科研市场上三代测序最常见的莫过于PacBio,辅以冉冉上升的新星Nanopore等,主打长读长策略,直击二代测序碎片化序列的软肋,在基因组de novo上表现不俗,错误率较高,但可被矫正。
回到我们今天的主题---数据格式上,一代测序主要是读取峰图文件后转化为Fasta格式;二代测序中illumina原始读取数据为BCL,下游分析中转化为Fastq格式;454下机序列为SFF格式,后续分析中转化为Fna-Qual格式使用;Ion Torrent下机序列为WELLS格式,下游分析中转化为Bam格式;三代测序的两大主流系统PacBio和Nanopore,其下机数据都以HDF5格式为基础,后续转化为Fastq格式进行下游分析。
不管一代、二代还是三代的数据分析中,原始下机数据都以二进制文件为主,原因无他,相比于文本文件,二进制文件在存储上更为经济集约。二进制文件本身是难于阅读的,并且很难改动,所以,我们可以乐观的认为,二进制文件造假的可能性是很低的。二进制的数据拿到之后,我们想要把数据转换成正常人能看懂的格式,这时身为文本文件的Fastq就应运而生了,Fastq文件会被用于质控及比对等后续分析。
总之,Fastq是当前最为主流认可的序列数据存储格式,不管哪一代测序技术,什么样的原始数据,都免不了要打上Fastq格式的烙印,Fastq文件的格式及使用已经成为高通量测序学习中当仁不让的第一站。
弄清下机数据格式,万里长征才算走完了第一步。
下一期,枫玲会接继续给大家分享有关人类基因组测序的数据分析干货。下期话题:序列比对,敬请期待~
【完】
- 本文固定链接: https://maimengkong.com/zu/1223.html
- 转载请注明: : 萌小白 2022年10月3日 于 卖萌控的博客 发表
- 百度已收录