番茄青枯雷尔氏菌致病力分化的全基因组分析Word文件下载.docx

番茄青枯雷尔氏菌致病力分化的全基因组分析Word文件下载.docx

《番茄青枯雷尔氏菌致病力分化的全基因组分析Word文件下载.docx》由会员分享，可在线阅读，更多相关《番茄青枯雷尔氏菌致病力分化的全基因组分析Word文件下载.docx（21页珍藏版）》请在冰豆网上搜索。

个特异获得基因、？

个特异丢失基因，在强致病力菌株FJAT-91中，共有？

特异获得基因，？

个特异丢失基因。

ANI进化分歧

组成成分的差异

关键词：

青枯雷尔氏菌，全基因组Solexa测序，denovo拼接，致病力分化

1前言

青枯雷尔氏菌（Ralstoniasolanacearum）引起的细菌性青枯病（Bacterialwilt）是危害最大、分布最广且造成损失最为严重的植物病害之一（Tsuchiyaetal.,2005）。

青枯雷尔氏菌主要分布在热带、亚热带和温带地区，环境适应能力强；

其寄主中有木本植物和草本植物，一年生植物和多年生植物，双子叶植物和单子叶植物（Janse,2005），对茄科植物如番茄、辣椒、茄子、烟草、马铃薯的危害尤为严重，成为生产的主要障碍之一（Yangetal.,2006），该菌亦可危害花生、香蕉、甘薯、生姜、桉树等植物（Sood,2005），寄主范围达到44个科的300多种植物（Alvarez,2005）。

因其侵染性强，寄主范围大，地域分布广，在水中和土壤环境中的生存能力强，青枯雷尔氏菌被认为是世界上危害最大植物病原菌之一（Genin,2010），美国在“911”之后一度将其列在十大农业恐怖微生物之首，多年来一直是相关学科的研究热点。

青枯雷尔氏菌表现出高度的种下分化的多态性，可按生理小种（Carmeilleetal.,2006）、生化型（Stefanietal.,2005）、血清型、基因型、致病型等指标进行种下分类（Thurston,2001）。

青枯雷尔氏菌的致病力分化极为显著的，自然环境中就存在强致病力、无致病力以及介于其中的过渡菌株。

强致病力菌株侵入寄主植物并导致寄主发病，无致病力菌株可侵染寄主植物却不引起寄主发病（Balabeletal.,2005）。

现有研究表明，无致病力菌株无法自然恢复为强致病力菌株（Robertsonetal.,2004），而强致病力菌株却会自然丧失致病性，弱化形成无致病力菌株，致病力丧失的原因可能是营养条件和极端气候引起的基因突变（Kandaetal.,2003）。

青枯雷尔氏菌存在的复杂多态性是其内部基因组的复杂性的体现，青枯雷尔氏菌能自然发生基因重组（Salanoubatetal.,2002），也存在从环境中获得遗传物质的水平基因转移（Guidotetal.,2009），模式菌株GMI1000中16%的基因可能来自于水平转移（Nakamuraetal.,2004），这些都是青枯雷尔氏菌的菌种复杂性的原因（Coupatetal.,2008）。

青枯雷尔氏菌的种内差异甚至让人怀疑这些菌株是否还应该属于一个物种。

2004年，强致病力的标准菌株GMI1000的基因组全序列被测序并发表（Nakamuraetal.,2004），迄今共有6个公开的强致病力菌株基因组序列（Remenantetal.,2010），其中包含4个完成图，2个草图。

通过对这些基因组的分析表明，青枯雷尔氏菌基因组呈现镶嵌结构（Genin,2004），一些持家基因及毒力基因构成大致的骨干，但不同的生态表型或寄主特异性则可能具有不同的基因组模块（Poueymiroetal.,2009）。

然而，目前尚无一个完整的无致病力菌株基因组完成图数据，数据的缺乏使得青枯雷尔氏菌复杂的致病力分化机制难以分析。

2材料与方法

2.1菌株的获得

通过常规组织分离方法，在番茄青枯病发生田块，通过大量采集与分离获得大量青枯雷尔氏菌菌株，结合TTC培养基鉴定和致病性检测，获得亲缘性关系极高的番茄无致病力菌株FJAT-1458和强致病力菌株FJAT-91。

青枯雷尔氏菌无致病力菌株和强致病力菌株的菌落形态以及致病力分化的检测结果见图a。

通过青枯雷尔氏菌其菌落形态来判定其致病性强弱，其中强致病性的青枯雷尔氏菌菌落形态表现为：

流动性强，中间为粉红色，并且白边比较宽，表面比较湿润。

弱致病性或无致病性的青枯雷尔氏菌菌落形态表现为：

无流动性，中间为暗红色，并且白边比较窄，表面比较干燥。

在相同接菌量的情况下，强致病性青枯雷尔氏菌采取剪叶方式在第4天就开始发病；

并且在第6天已经完全发病（死亡率达到100%）。

而弱致病性青枯雷尔氏菌在这段期间剪叶均不发病。

2.2高通量测序

本文采用Solexa-illumina新一代测序技术分别对青枯雷尔氏菌FJAT-91和FJAT-1458的全基因组进行shotgun测序，其中前者是强致病力菌株，而后者是无致病力菌株，测序平台为GAIIx。

Solexa测序读段（shortreads）的读长为54bp，分别构建了两种双末端测序文库（Library），插入序列长度（Insertlength）分别为300bp和2500bp。

其中300bp的插入文库采用Paired-End测序，两端读段的方向相向（inward）；

2500bp文库采用Mate-Pair测序，两端读段的方向相反（outward）。

图1.生物信息学分析的技术路线

图2.青枯雷尔氏菌致病力分化

注：

青枯雷尔氏菌强致病力菌株的菌落形态（a：

普通相机；

b：

体视显微镜），青枯雷尔氏菌无致病力菌株的菌落形态（d：

e：

体视显微镜），致病力分化的检测结果（c：

强致病力菌株；

f：

无致病力菌株）

2.3denovo拼接

2.3.1拼接前的质量控制

在对高通量读段进行拼接之前，需要进行一系列质量控制，包括：

格式转换，质量评估，质量修剪和长度过滤。

Solexa高通量测序生成原始数据为一套qseq文件，需要根据标签信息转换为fastq格式，并根据illumina的自动筛选结果输出高质量的读段。

qseq转换为fastq文件的perl脚本采用了多线程编程技术，可根据标签信息，并行式地从各个qseq文件中提取目标样本的读段，并统计目标样本测序的总数据量。

读段的质量值采用了illumina1.3+标准，即PHRED值0至62分别由第64位ASCII码字符（@）至第126位字符（~）表示。

本文首先利用SolexaQA（Coxetal.,2010）对全部读段的测序质量进行评估，评测结果可为质量修剪提供依据。

本文在质量修剪过程中，首先统一剪除3'

末端最后一个碱基，使所有读段只为53bp；

继而采用Modified-Motttrimming算法（为CLCGenomicsWorkbench所采用质量修剪算法）进行Perl编程，对所有读段进行质量修剪和过滤，过程如下：

首先通过公式2，将PHRED质量值换算成单碱基测序错误的概率（Pe），接着对每个碱基都计算了一个新值Val，令Val=Limit-Pe，本文将Limit设为0.05。

再设置一个新值sum，初始为0；

从读段的末端开始扫描序列，每扫描到一个碱基，就把碱基对应的Val值累加，赋予sum，即sum=sum+Val；

若碰到sum<

0，则令sum=0；

从第一次出现sum>

0的碱基开始到出现sum最大值的碱基结束，这段序列就是质量修剪后（trimmed）保留下来的高质量（Highquality）读段，而低质量的序列末端已被切除。

若扫描整条序列后，sum仍为0，则舍弃整个读段。

经过质量修剪之后，对含有超过2个N（BaseCalling模糊不清的碱基）的读段进行修剪，只取最长的一段少于等于2个N的读段。

最后根据长度阀值进行过滤，本文只保留读段长度大于阀值35bp的读段。

若PE配对的两条读段中仅其中一条被过滤，则将另一条被保留的读段另存至一个SE（SingleEnd）文件中，使之不影响PE读段的对齐。

公式1

公式2

2.3.2ABySS和Velvet拼接

在拼接之前，还需要进行以下处理：

Mate-Pair文库（MPLibrary），即本文的2.5kb文库，双末端配对读段的方向为outward，需要将forward和reverse读段都进行反向互补，使两读段方向转为inward。

velvet和ABySS都只识别inward的PE读段，转换之后会使拼接质量提高。

本文先分别利用ABySS（Simpsonetal.,2009）和Velvet（Zerbino,Birney,2008）单独各自对青枯雷尔氏菌FJAT-91和FJAT-1458进行拼接。

对Kmer值从21至上限35（读段长度过滤的阀值）之间的每个奇数都进行调试，确认达到最佳拼接效果时的Kmer值为35bp；

另外，Velvet的参数min_pair_count和ABySS的参数n意义是一样的，都设为20，表示构成一个框架序列Scaffold至少需要20对读段；

ABySS不需要再设其他参数，Velvet需要参数exp_cov，表示期望的覆盖度（Coverage），可设为auto，还需要指定插入文库的长度及标准差SD，SD的设定只需要取一个尽量大的值，不需要调试，本文设为20%，即300bp的SD为60bp，2500bp的SD为500bp。

本文进一步将ABySS第一次拼接（Kmer=35）得到的框架序列（Scaffold）拆分为重叠群（Contig），Contig长度长于Velvet长序列拼接的长度限制（32700bp）时，将ABySS的Contig拆分为三分之二互有重叠的两段序列。

最后将ABySS的Contig序列作为长序列（longsequence）和所有短序列（shortread）一起通过Velvet进行二次拼接，Velvet的其他参数不变（Kmer=35）。

2.3.3读段定位

将读段和拼接后的框架图（Scaffold）进行比对，读段的定位（Mapping）信息可用来计算读段的覆盖率（Coverage）、SNP频率和测序错误率（Error）。

本文利用大规模高通量读段快速联配作图工具Bowtie（Langmeadetal.,2009）分别将青枯雷尔氏菌FJAT-91和FJAT-1458质量修剪后的读段，与Velvet-ABySS拼接得到的框架图进行比对。

首先，选择序列长度大于300bp的Scaffold，并按长度排序以及重命令，如RS91va_0001。

然后，通过bowtie-build建立框架图的ebwt索引。

最后，选择Bowtie比对模式为“-v2”，即半全局联配模式（全局的读段-vs-局部的染色体），允许有2个错配，bowtie尚不支持Gapped比对，比对结果中不会含有InDel；

其他参数有--phred64-quals，即本文的读段质量标准为illumina1.3+；

比对结果保存为SAM格式（参数-S）；

此外，采用Paired-End模式（参数-q-1m1-2m2），不能达到理想的比对效果，因而本文未采用该模式。

2.3.4计算覆盖率（深度和广度）

Bowtie的比对结果保存为SAM格式，可利用Perl语言读取SAM文件计算定位的读段（Mappedread）总数量和碱基总数。

利用SAMTools（Lietal.,2009）处理SAM文件（文本文件），过滤未定位的读段，并转换为BAM格式（压缩的二进制SAM文件），对BAM文件排序并建立索引，最后将多个基于同一个框架图的BAM文件合而为一。

本文利用BEDTools（Quinlan,Hall,2010）的genomeCoverageBed程序计算per-basecoverage，即单碱基覆盖率，需要参数“-d”。

根据perbasecoverage通过Perl编程计算所有定位读段对整个框架图（含Gap）的总覆盖率，包括测序深度（Depth）和广度（Breadth）。

测序深度为所有定位读段覆盖框架图的倍数，广度为所有定位读段覆盖框架图的完整程度，在本文研究中，框架图中未覆盖读段的部分即为Gap区域。

通过扫描perbasecoverage，计算各个框架序列（Scaffold）的覆盖度（深度Depth和广度Breadth），同时计算各个重叠群（Contig，不含Gap）的测序深度。

可以根据覆盖率的高低评估框架图的拼接质量，滤去覆盖率过低（本文设为30×

）的Scaffold；

而覆盖率数倍于平均覆盖率的序列，则可能是重复序列。

2.3.5检测SNP和测序错误

读段的Bowtie比对结果已转换为SAMTools处理之后的BAM格式，利用SAMTools的mpileup生成pileup文件，在此文件中包括了每个匹配位点的定位详情，pileup格式说明详见SAMTools网站上的说明：

通过perl语言编程来处理pileup文件。

首先从pileup文件中列出所有差异位点（读段与参照序列之间有错配）。

其次，对这些差异位点进行SNP过滤，SNP或者Error的区别只是频率不同。

某个位点发生多次测序错误的概率会随着错配次数的增多而急剧降低，错配次数越多，说明该位点存在SNP的概率越高。

比如说，参照序列为碱基“A”，读段错配碱基为“T”和“G”，A至T的次数和频率都很高，则有存在SNP的可能性就高，而A至G的次数仅一次，则完全不能排除是测序错误的可能性。

因此，进行SNP过滤需要设定一系列阀值，本文设定的阀值如下：

差异位点的perbasecoverage在30以上，500以下，错配次数（绝对数量）在5次以上（包括5次），SNP频率大于等于2%。

最后，对不符合SNP筛选条件的错配当作测序错误，错配数进行累积，从而计算测序错误率。

测序的错误率为所有测序错误碱基数除以所有定位读段的总碱基数。

另外，框架图中少数几个位点需要进行SNP校正。

本文进行读段定位的参照序列是新拼接的框架图序列，来自于这些定位读段，在计算差异位点的错配比率时，部分SNP频率高于参照碱基的频率，并且在所有SNP中频率最高，将该SNP和参照碱基交换，进行校正。

对框架图进行SNP校正之后，再计算SNP频率和错误率。

2.3.6框架图拼接后的质控统计

青枯雷尔氏菌FJAT-91和FJAT-1458基因组拼接之后，需要对其拼接结果（框架图）进行评估，从而判断结果是否理想。

首先，通过框架序列（Scaffold）的长度阀值和最低覆盖率对框架图进行质量控制，本文设定的长度阀值为300bp，最低覆盖率为30×

。

然后，对质控之后的框架图进行评估，计算以下几个数值：

Scaffold的数量，长度总和（基因组规模），Scaffold长度的最大值和平均值，Gap的比例，Scaffold的N50、N90序列数量和N50、N90长度值，Contig的数量，Contig的长度总和，Contig长度的最大值和平均值，Contig的N50、N90序列数量和N50、N90长度值。

其中N50、N90值的概念如下：

将各个序列按长度大小排序，从大至小逐一扫描各个序列的长度值，进行累积，当该累积值第一次超过所有序列总长的50%时，此时扫描到的序列，其长度值即为N50值，而此时已扫描的序列数量即为N50序列数量。

N90值亦同理。

2.4基因预测与注释

本文利用原核生物基因预测软件Glimmer3、GeneMarkS和Prodigal分别在青枯雷尔氏菌FJAT-91和FJAT-1458基因组框架图的重叠群序列（Contig）上预测编码基因序列。

因本文测得的框架图尚不完整，含有较多的Gap，基因预测需包括未达全长的基因片段（Fragment）。

基因功能注释，将FJAT-91和FJAT-1458的基因分别与Swiss-Prot和trEMBL数据库进行BLASTP同源比对，选择最相似的比对的结果，并且要求序列比对相似性在30%以上，比对双方中较短的那条序列的全长，联配长度需占其60%以上，期望值E-value在0.001以上。

通过Swiss-prot的GO关联（GOA，GOAssociation），分别对FJAT-91和FJAT-1458中具有GO注释的基因利用WEGO进行GO功能分类。

蛋白质结构域分析，利用Hmmer3的hmmscan程序将FJAT-91和FJAT-1458所有预测基因的蛋白质序列与本地化PFAM数据库进行比对，从而鉴定这些基因的蛋白质结构域。

非编码RNA的预测包括核糖体RNA（rRNA）和转运RNA（tRNA）。

rRNA可通过青枯雷尔氏菌GMI1000中rRNA序列进行同源预测，tRNA利用rRNAScan-SE和aragorn进行预测。

重复序列注释，通过RepeatModeler（调用RECON和RepeatScout）分析FJAT-91和FJAT-1458的框架图序列，发现重复序列，构建物种特异的重复序列数据库，联同RepBase，利用RepeatMasker检测框架图中的重复序列，通过和已有的TE数据库进行比对，对重复序列进行注释。

2.5比较基因组学研究

青枯雷尔氏菌CMR15、PSI07、CFBP2957、GMI1000和Po82全基因组完成图（CompletedGenome），青枯雷尔氏菌IPO1609、UW551和MolK2基因组草图（GenomeDraft），以及青枯雷尔氏菌近缘物种Ralstoniapickettii12J全基因组完成图，均可从NCBI或者GenoscopeRalstoniaScope网站下载数据，后者的网址如下：

s.fr/agc/microscope/about/collabprojects.php?

P_id=67

将新测序的青枯雷尔氏菌强致病力菌株FJAT-91和无致病力菌株FJAT-1458联同以上9个菌株的基因组一起进行比较基因组学研究。

2.5.1双向最佳比对

在全基因组水平上，两个物种之间的蛋白质序列进行双向最佳比对（BBH，BidirectionalBestHit），可以鉴定得到这两个物种之间的直系同源基因（OrthologousGenes）。

利用BLAST+程序包中的BLASTP程序，将以上每个菌株全基因组中的所有蛋白质序列分别和其他10个菌株进行BLASTP，只选择最佳匹配项，并要求序列比对相似性在30%以上，比对双方中较短的那条序列的全长，联配长度需占其60%以上，期望值E-value在0.001以上。

通过Perl编程，各个菌株两两之间的最佳比对结果进行聚类……

2.5.2全基因组比对（megablast）和基因组岛鉴定

2.5.3计算AAI并构建物种进化树

通过蛋白质序列双向最佳比对（BBH）鉴定各个物种两两之间的直系同源基因（OrthologousGenes），从而计算这些直系同源基因的平均氨基酸一致性（AAI，AverageAminoacidIdentity）。

各菌株之间的分歧程度通过100%减去AAI值来表示，以此生成距离矩阵，并保存为Phylip所能读取的格式。

之后，利用phylip的neighbor程序，采用邻接法（NB）构建物种进化树。

3结果与分析

3.1框架图拼接结果分析

3.1.1数据量和测序深度

经过Solexa高通量测序，青枯雷尔氏菌强致病力菌株FJAT-91共产生了1.3Gb数据量（测得碱基总数），无致病力菌株FJAT-1458共产生了1.1Gb数据量。

经过illumina自动筛选之后，FJAT-91的有效读数，即高质量读段的数量，为2109万个，碱基总数为1.1Gb，占了原始数据的87.6%，测序深度为203×

（已知青枯雷尔氏菌基因组平均规模为5.6Mb）；

而FJAT-1458的有效读数为1751万个，碱基总数为0.9Gb，是原始数据的87.5%（详见表1），测序深度为169×

如表1所示，读段经过质量修剪之后，青枯雷尔氏菌FJAT-91和FJAT-1458实际用于拼接的测序深度分别为156×

和131×

，各自修剪了23.2%和22.5%的原始序列，读段的平均长度分别为49.8bp和50.0bp，原始读段长度为54bp。

在denovo拼接之后，将读段定位至框架图上，未能联配至框架图上的读段，可以认为这些读段没有在拼接中用到。

读段定位结果显示，实际用于拼接的读段覆盖率深度在FJAT-91和FJAT-1458中分别有129×

和103×

，分别下降了17%和21%，用于拼接的读段分别仅占原始数据的56%和53%。

表1.青枯雷尔氏菌FJAT-91和FJAT-1458读段数量和碱基总量

FJAT-91（virulent）

FJAT-1458（avirulent）

Rawdata

readlength

54bp

readsnumber

24072750

20001978

totalbases

1299928500bp

1080106812bp

genomesize

avg.5.6Mb

depth1）

232×

193×

Illuminafiltered

21094106

17507820

1139081724bp

945422280bp

depth

203×

169×

Qualitytrimmed

35~53bp,mean=49.8bp

35~53bp,mean=50bp

17501038

14740010

871237877bp

736276730bp