340 likes | 582 Views
青枯雷尔氏菌的深度测序、Scaffold/Contig的组装和基因注释. 唐唯其 2010年9月 福建农科院农业生物资源研究所微生物研究中心 http://fjat.yzxz.com/. 双击添加标题文字. 1. 深度测序的简要介绍. 青枯雷尔氏菌RS98和RS91序列拼接和注释. 2. 3. 比较基因组学的分析(部分). 青枯雷尔氏菌生物信息学的实验内容. 序列拼接 Assembly 对Scaffold和Contig进行概述性分析 原核生物基因预测 及功能注释
E N D
青枯雷尔氏菌的深度测序、Scaffold/Contig的组装和基因注释青枯雷尔氏菌的深度测序、Scaffold/Contig的组装和基因注释 唐唯其 2010年9月 福建农科院农业生物资源研究所微生物研究中心 http://fjat.yzxz.com/
双击添加标题文字 1 • 深度测序的简要介绍 • 青枯雷尔氏菌RS98和RS91序列拼接和注释 2 3 • 比较基因组学的分析(部分)
青枯雷尔氏菌生物信息学的实验内容 • 序列拼接 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.
新一代测序技术的优缺点 • 优点: • 体外构建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、CeleraAssembler(wgs) • CAP3 and PCAP • Newbler(454) 图片来源于http://www.illumina.com/Documents/products/technotes/technote_denovo_assembly.pdf
将qseq数据分别转换成velvet和SOAPdenovo所要求的FASTA文件将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)
拼接结果的概述 Source from "Genomes of three tomato pathogens within theRalstonia solanacearum species complex reveal significant evolutionary divergence"
拼接结果的概述2 Corverge计算公式:
无题 面临的问题? 同一套数据,采用不同的拼接软件,得到的结果却有较大的差异。 需要进一步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
基因预测的结果 ● 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(术语)
青枯雷尔氏菌比较基因组学的研究 • 数据源 • 精细图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鉴别直系同源基因,没有匹配的是特异基因。
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) • 其他的文章: • ?尚缺