Transcript 青枯雷尔氏菌的深度测序
青枯雷尔氏菌的深度测序、Scaffold/Contig的组 装和基因注释 唐唯其 2010年9月 福建农科院农业生物资源研究所微生物研究中心 http://fjat.yzxz.com/ •深度测序的简要介绍 •青枯雷尔氏菌RS98和RS91序列拼接和注 释 •比较基因组学的分析（部分） 青枯雷尔氏菌生物信息学的实验内容 序列拼接 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（第三代） 测序设备 科技的发展 过去： – 13年，30亿，3Gb，HGP人类基因组计划（vs 阿波罗登月计划、曼哈顿核弹计划） 现在： – 炎黄计划，1000 genome计划（已完成三个先导项目） – 四周，48000美元，Heliscope单分子测序仪 – HiSeq 2000，x30，1万美元，100bp，50万美元左右>10万美元。 – 4400美元，Complete Genomics纳米阵列技术（第三代测序技术），2009年底。定价 2万美元。 将来： – 2013年、1000美元、15分钟、微波炉，个人基因组测序仪 传统的Sanger法测序原理和循环芯片测序法（Cyclic-Arrary Sequencing）原理图示 . Jay Shendure & Hanlee Ji. (2008) Next-generation DNA sequencing. Nature Biotechnology, 26(10):1135-1145 454焦磷酸测序法原理 Solexa测序原理 SOLiD测序原理 Helicos单分子测序原理 Complete Genomics 纳米阵列技术 第三代测序技术-纳米孔技术原理 新一代测序技术的优缺点 优点： – 体外构建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 – 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/techno te_denovo_assembly.pdf de novo Assembly算法 将qseq数据分别转换成velvet和SOAPdenovo所要求的FASTA文件 velvet和SOAPdenovo接受的输入文件主要有两种格式： – FASTA • >seq1 • 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 [LIB] 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 – 打开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 Gap Num Length Length Scf Velvet SOAPdenovo 971 10255 Total Num Avg. Max 5.9Mb 3661 1623bp 58kb 402kb 923kb Con 5.5Mb 4332 1279bp 53kb Scf 8.5Mb 13380 638bp 15kb Con 7.6Mb 23635 322bp 15kb Corverge计算公式： Coverage GC Cotent x49.9 65.9% x84.1 59.7% (length coverage) length 无题 面临的问题？ 同一套数据，采用不同的拼接软件，得到的结果却有较大的差异。 需要进一步Check！策略： 1，整合多种拼接结果（取共有的结果）；2，Velvet的结果再拼接； 3，设定一个简单的标准：覆盖度x10以上，长度200bp以上。 4，基于Mapping的拼接（GMI1000等精细图作参考） 在青枯雷尔氏菌基因组研究上，做别人没有做过的，基于Solexa测序。 – 检测SNP（单核苷酸多态性），分析等位基因和菌株杂合率。 – 检查Scaffold中的重复序列，评估repeat对Short Reads Assemble的影响 原核生物基因预测的方法和工具 工具： – Glimmer3 • 在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预测 6.2Mb 68.0% 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） 方法 – 双向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) 其他的文章： – ？尚缺 感谢您的关注！