青枯雷尔氏菌的深度测序、Scaffold/Contig的组装和基因注释 唐唯其 2010年9月 福建农科院农业生物资源研究所微生物研究中心 http://fjat.yzxz.com/

双击添加标题文字 深度测序的简要介绍 青枯雷尔氏菌RS98和RS91序列拼接和注释 2 3 比较基因组学的分析（部分） 1

青枯雷尔氏菌生物信息学的实验内容 序列拼接 Assembly 对Scaffold和Contig进行概述性分析 原核生物基因预测及功能注释 青枯雷尔氏菌精细图GMI1000、PSI07、CFBP2957、CMR15，草图MolK2、IPO1609、UW551，和RS98、RS91比较基因组学分析 共有基因和特异基因 同源性和共线性分析 基因组尺度的进化分歧 个例分析：识别并分析致病基因 这是我们的实验内容。这篇文章的结构，我们已经构建好了，只要逐步把实验完成就好了。 但是，我还没有全部做完。从刚接手到现在，中间也已经转了几套方案。 这次主讲-序列拼接，以及对上游的测序技术进行介绍。

深度测序、Deep Sequencing、高通量测序、第二代测序、新一代测序、Next-Generation Sequencing Roche 454 焦磷酸测序 http://www.454.com/ illumina-Solexa http://www.illumina.com/ ABI SOLiD（supported oligo ligation detetion） http://www.appliedbiosystems.com.cn/ Helicos 单分子测序 http://www.helicosbio.com/ 纳米孔测序 nanopore sequencing（第三代） 深度测序介绍：http://blog.sina.com.cn/s/blog_4dd865a70100g9n9.html

测序设备

科技的发展 过去： 现在： 将来： 13年，30亿，3Gb，HGP人类基因组计划（vs 阿波罗登月计划、曼哈顿核弹计划） 炎黄计划，1000 genome计划（已完成三个先导项目） 四周，48000美元，Heliscope单分子测序仪 HiSeq 2000，x30，1万美元，100bp，50万美元左右>10万美元。 4400美元，Complete Genomics纳米阵列技术（第三代测序技术），2009年底。定价2万美元。 将来： 2013年、1000美元、15分钟、微波炉，个人基因组测序仪 2009年Heliscope测序仪：Single-molecule sequencing of an individual human genome，《Nature》，2009 Complete Genomics纳米阵列技术：Human genome sequencing using unchained base reads on self-assembling DNA nanoarrays，《Science》2009 Nanopore DNA sequencing with MspA，2010，《PNAS》 测序费用逐年减少 仪器小型化，设备价格在减少。

传统的Sanger法测序原理和循环芯片测序法（Cyclic-Arrary Sequencing）原理图示 生命奥秘新一代测序技术专题：http://www.lifeomics.com/?p=21006 传统的Sanger法测序和循环芯片测序法（cyclic-array sequencing）的原理： a是Sanger法测序的图示，b是循环芯片测序法的图示。 将毛细管（凝胶电泳）比作电子管，芯片比作晶体管集成电路。 Jay Shendure & Hanlee Ji. (2008) Next-generation DNA sequencing. Nature Biotechnology, 26(10):1135-1145.

454焦磷酸测序法原理 焦磷酸测序的原理是通过检测碱基掺入时发出的光来进行测序的 微乳PCR和桥式PCR

Solexa测序原理 可逆性末端终止反应

SOLiD测序原理

Helicos单分子测序原理 HeliScope测序仪最大的特点是无需对测序模板进行扩增，它使用了一种高灵敏度的荧光探测仪直接对单链DNA模板进行合成法测序。首先，将基因组DNA切割成随机的小片段DNA分子，并且在每个片段末端加上poly-A尾。然后通过poly-A尾和固定在芯片上的poly-T杂交，将待测模板固定到芯片上，制成测序芯片。最后借助聚合酶将荧光标记的单核苷酸掺入到引物上（图6d）。采集荧光信号，切除荧光标记基团，进行下一轮测序反应，如此反复，最终获得完整的序列信息。根据最近的报道，经过数百轮这种单碱基延伸可以获得25bp或更长的测序长度。

Complete Genomics 纳米阵列技术

第三代测序技术-纳米孔技术原理 专题：http://www.lifeomics.com/?p=21027 新型纳米孔测序法（nanopore sequencing）是采用电泳技术，借助电泳驱动单个分子逐一通过纳米孔来实现测序的。由于纳米孔的直径非常细小，仅允许单个核酸聚合物通过，因而可以在此基础上使用多种方法来进行高通量检测。此外，纳米级别的孔径保证了检测具有良好的持续性，所以测序的准确度非常高。对于长达1,000个碱基的单链DNA分子、RNA分子或者更短的核酸分子而言，根本无需进行扩增或标记就可以使用纳米孔测序法进行检测，这使得便宜、快速地进行DNA测序成为可能。如果对现有纳米孔测序法进行进一步发展和改进，那么它将有望成为第三代测序技术（也可称为下、下一代测序技术），从而帮助人们实现24小时内只花费1,000美元完成二倍体哺乳动物基因组测序这一目标。

新一代测序技术的优缺点 优点： 体外构建DNA文库，随后体外克隆扩增DNA片段。这就解决了传统Sanger测序技术中好几个限制平行测序规模的瓶颈问题，比如转化大肠杆菌以及阳性克隆挑选等问题。 基于芯片的测序方式能极大地提高平行测序的能力。 试剂消耗少，测序费用低。 缺点： 测序片段长度短。但随着技术的发展，单个测序反应获得的片段长度也越来越长。Sanger法一个反应一般有500～800bp。Solexa测序从一开始的20～30bp到现在的75bp，454焦磷酸测序从一开始的100bp到现在的400bp。 准确率不高。平均来说，新一代测序的准确率比传统测序技术低10倍。但可以通过高覆盖度来验证测序错误。

测序结果--本实验的原始数据 测序返回的结果有 询问之后，对方的说明： Solexa测序返回的short reads共有两套数据，s_1和s_7。 每套数据分别有360个文件构成，s_1分为s_1_1系列的120个qseq文件、s_1_2系列的120个qseq文件和s_1_3系列的120个qseq文件；s_7同理。 其中1和3系列是的short reads长度为54bp，而2系列的长度只有7bp。 每个qseq文件均包含20万条左右short read，其中s_1的有效序列占总数的86%左右，在s_7中，该平均值为87%。 已经拼接过的Scaffold序列：RS98，7.2Mb（有效6.2Mb），RS91，6.6Mb（有效5.5Mb）。 询问之后，对方的说明： s_1_1和s_1_3是对应的双末端序列，library长度分为两种，分别是300bp和2500bp。 s_1_2系列的7bp的短序列是测序中不同样本的编号。 s_1系列中，编号为CACTTGAA的是RS91的300bp插入文库 s_7系列中，编号为CACTTGAA的是RS91的2500bp插入文库 s_1系列中，编号为CGATCAGA的是RS1458的300bp插入文库 s_1系列中，编号为CCTTGTAA的是RS1458的2500bp插入文库 除了这些序列外，其他编号的序列是同次测序中的其他样品。

illumina测序结果的qseq数据格式说明 Illumina pipeline versions 1.3 and later说明  Machine name: unique identifier of the sequencer.  Run number: unique number to identify the run on the sequencer.  Lane number: positive integer (currently 1-8).  Tile number: positive integer.  X: x coordinate of the spot. Integer (can be negative).  Y: y coordinate of the spot. Integer (can be negative).  Index: positive integer. No indexing should have a value of 1.  Read Number: 1 for single reads; 1 or 2 for paired ends.  Sequence (BASES)  Quality: the calibrated quality string. (QUALITIES)  Filter: Did the read pass filtering? 0 - No, 1 - Yes.

重复太多，不正常，需要重新拼接

Short Reads Assembly de novo reference 参照基因组（Mapping） 传统的拼接软件 velvet（EBI），基于de Bruijn graph SOAPdenovo（BGI华大），基于图论 Euler、Euler-SR、ABySS、ALLPATHS（Broad Institute）、SHARCGS、SSAKE、VCAKE、Edena、Forge reference 参照基因组（Mapping） SHORE和genomemapper（1001genome.org），MAQ（Sanger Center）、SOAP2（BGI）、ELAND（illumina） BFAST、Bowtie、BWA（Maq） 传统的拼接软件 phred+phrap+consed AMOS、Celera Assembler(wgs) CAP3 and PCAP Newbler（454） 图片来源于http://www.illumina.com/Documents/products/technotes/technote_denovo_assembly.pdf

de novo Assembly算法

将qseq数据分别转换成velvet和SOAPdenovo所要求的FASTA文件 atcgactg >seq2 atcgcgat FASTQ，多了一份测序质量的信息。 不同文库插入长度的双末端（Paired-End）序列要根据不同的要求转换成相应的fasta文件 velvet：300bp和2500bp各一个文件，Paired-End对应序列都放在同一个文件中，第1，2，3，4……，其中第1和第2是一个Paired-End的两端short reads，第3和第4是另一个Paired-End的双末端，以下依次排列。 SOAPdenovo：300bp和2500bp分开，文库长度300bp的分两个文件，f1和f2，f1里的第一个short read对应f2里的第一个short read，以下依次对应排列。2500bp亦同理。 qseq2fasta.ShortPaired2.pl（for velvet）=> fasta1-2.SOAPdenovo.pl（for SOAPdenovo）

Velvet 使用说明 执行命令： 参数说明： 输出结果（在保存目录下）主要有： velveth 保存目录 31 -fasta -shortPaired RS98.300.fasta -shortPaired2 RS98.2500.fasta velvetg 保存目录 -cov_cutoff auto -exp_cov auto -read_trkg yes -amos_file yes -unused_reads yes -ins_length 300 -ins_length2 2500 参数说明： 31是K-Mer的值（默认31）， -fasta是读取序列的格式。 通过-shortPaired 和 -shortPaired2分别指定两个fasta文件。 cov_cutoff和exp_cov 设定为auto，自动取相应的覆盖率。 read_trkg amos_file设定为yes，额外产生一个afg文件（显示short reads拼接状况）。 unused_reads设定为yes，列出所有未拼接至contig或scaffold的short reads。 ins_length 300和ins_length2 2500分别指定两种Insert Length。 输出结果（在保存目录下）主要有： contigs.fa（FASTA格式） velvet_asm.afg UnusedReads.fa stats.txt LastGraph Log

SOAPdenovo 使用说明 首先需自定义一个配置文件（可参考说明书附录的样文） 执行命令如下： 输出的主要结果有： max_rd_len=54 [LIB] avg_ins=300 reverse_seq=0 asm_flags=3 pair_num_cutoff=3 map_len=32 f1=RS98.seq1.300.1.fasta f2=RS98.seq1.300.2.fasta avg_ins=2500 reverse_seq=1 asm_flags=2 rank=2 pair_num_cutoff=5 map_len=35 f1=RS98.seq2.2500.1.fasta f2=RS98.seq2.2500.1.fasta 执行命令如下： soapdenovo all –s config_file –K 25 –R –o graph_prefix 说明：config_file即以上自定义的配置文件，graph_prefix指的是输出文件的前缀（比如“RS98”） 输出的主要结果有： RS98.scafSeq、RS98.contig、RS98.gapSeq、RS98.scaf

Visualization 可视化（略） AMOS Hawkeye Gbrowse UCSC Browser 打开afg文件，显示拼接结果 安装需要：QT环境，在win下还需要cygwin Gbrowse UCSC Browser Staden Tools (GAP4, TGAP) illumina GenomeStudio Mapview Xmatchview（Python） IGV（Integrative Genomics Viewer）

拼接结果的概述 Scaffold /contig base Contig Gap Gap Number GC Content 旧 270/1584 7,182,105 bp 6,210,616 bp 971,489 bp 1314个 68.7% Velvet 3361/4332 5,941,578 bp 5,539,152 bp 402,426 bp 971个 65.9% Source from "Genomes of three tomato pathogens within the Ralstonia solanacearum species complex reveal significant evolutionary divergence"

拼接结果的概述2 Corverge计算公式： Gap Length Coverage GC Cotent Num Total Avg. Max Velvet 971 402kb Scf 5.9Mb 3661 1623bp 58kb x49.9 65.9% Con 5.5Mb 4332 1279bp 53kb SOAPdenovo 10255 923kb 8.5Mb 13380 638bp 15kb x84.1 59.7% 7.6Mb 23635 322bp Corverge计算公式：

无题 面临的问题？ 同一套数据，采用不同的拼接软件，得到的结果却有较大的差异。 需要进一步Check！策略： 1，整合多种拼接结果（取共有的结果）；2，Velvet的结果再拼接； 3，设定一个简单的标准：覆盖度x10以上，长度200bp以上。 4，基于Mapping的拼接（GMI1000等精细图作参考） 在青枯雷尔氏菌基因组研究上，做别人没有做过的，基于Solexa测序。 检测SNP（单核苷酸多态性），分析等位基因和菌株杂合率。 检查Scaffold中的重复序列，评估repeat对Short Reads Assemble的影响

原核生物基因预测的方法和工具 工具： Glimmer3 GeneMark：GeneMarkS getorf (EMBOSS) 在Linux 64位下安装不成功，仅在NCBI Bacteria Tools Glimmer3网页上运算。 GeneMark：GeneMarkS 优点：可以得到基因片段（Fragment），应用了HMM算法。 getorf (EMBOSS) 仅仅用于搜索ORF perl的批量调用和改编输出结果 GeneMark需要对单个Contig/Scaffold进行计算，然后将结果整合在一起。 执行batch_run_gm.pl rRNA识别和tRNA预测 通过同源比对寻找rRNA（16S、5S、23S） aragorn tRNAscan-SE

基因预测的结果 基因组平均规模 GC含量平均值 基因数量 基因密度均值 基因的平均长度 六种已测序青枯菌 5.68Mb 66.6% 5330个 84.20% 921bp RS98 Glimmer3预测 6.2Mb 68.0% 6017个 66.27% 684bp RS98 GeneMarkS预测 7553个 83.35% 685bp ● Glimmer3和GeneMarkS的结果相比较，发现共有5404个同时预测到的基因，但这些基因的起始位置并不完全相同。Glimmer3尚有613个基因，未被GeneMarkS预测到，而GeneMarkS还有2149个基因是Glimmer3没有预测到的，其中有1593个是fragment。 ● 获得65个tRNA，同源比对得到rRNA一组（16S+5S+23S），但多了一个5S，疑为拼接错误所致。

基因功能的注释方法 和NCBI Bacteria WGS（Whole genome Sequence）数据库所有的蛋白序列进行比对，共有1162个物种(2203个refseq ID)。 数据载自NCBI FTP，下载all.faa.tar.gz文件。 筛选：取最佳匹配（Blast的Score分值最高），且匹配片段占其中较短的那条序列总长的60%以上，相似性在30%以上的结果，作为同源注释的依据。 仅和GMI1000进行比对，通过GMI1000的同源基因进行注释。 选择标准：Best Hit，60%，30%。（附：80%，40%） 和UniProt进行比对 Swiss-Prot和trEMBL 和GO进行比对 三个方面：Biological Process、Cellular Component、Molecular Function 和致病相关的GO Terms（术语）

RS98最佳匹配项所属物种 物种名 最佳匹配原始结果 筛选后结果 All 7553 6879 GMI1000 6170 6039 PSI07 326 318 CFBP2957 120 119 Ralstonia pickettii 69 65 Ralstonia eutropha 19 17 Burkholderia 124 106 Others 725 215

青枯雷尔氏菌比较基因组学的研究 数据源 下载数据 实验内容 精细图finished：GMI1000、PSI07、CMR15、CFBP2957 草图draft：MolK2、IPO1609、UW551 ours：RS98和RS91 下载数据 从NCBI上下载全基因组序列、所有基因的核酸以及蛋白序列、gb文件等。 GeneScope（论文的官方网站），cds、pep、genome、gff、ncrna、repeat 实验内容 区分各共有和特异基因 同源性、共线性分析：Megablast和Mummer，以及作图！ 基因组尺度的进化分歧：ANI Average nucleotide identity

BLAST双向最佳匹配（BBH） 方法 BBH鉴别直系同源基因，没有匹配的是特异基因。 CFBP2957 CMR15 PSI07 双向BLASTP最佳比对（通过设定参数max_target_seqs为1），并进行筛选。 筛选标准：选择最佳匹配（Score分值最大），并要求匹配片段长度占较短的比对序列的60%以上，相似性在30%以上。 BBH鉴别直系同源基因，没有匹配的是特异基因。 CFBP2957 CMR15 PSI07 GMI1000 基因数目 5310 5149 5247 5635 RS98-vs- 6280 6636 6106 7125 -vs-RS98 4154 4209 4219 4668 双向最佳匹配（BBH） 3681 3839 3695 4251

Virulence Factors致病基因的生物信息学分析 分析BMC Genomics 2010, 11:379这篇文章归纳的约130个致病基因 TTE(type III effectors) 该类基因的生物信息学预测方法(PLoS pathog 2009, 5: e1000376) EPS(the exopolysaccharide biosynthetic genes) CWDE(the cell wall-degrading enzyme) 其他的文章： ？尚缺

感谢您的关注！