块茎作物马铃薯的基因组序列和分析

　　库是由DM基因组DNA构建的，并在Sanger，Illumina Genome Analyzer 2（GA2）和Roche 454平台上进行了测序（请参阅补充文本）。使用Sanger平台对BAC库和三个Fosmid库进行了测序。对于Illumina GA2平台，我们产生了70.6 GB的37-73 bp配对末端读数，其中16个库的插件长度为200-811 bp（补充表21和22）。我们还产生了18.7 GB的Illumina Mate Pair库（2、5和10 KB插入尺寸）。总共生成了7.2 GB的454个单端数据，并应用于间隙填充以改善组件，其中4.7 GB（12,594,513读）被纳入了最终组件中。对于8和20 KB 454配对读数，分别代表0.7和1.0 GB的原始数据，分别为90.7 MB（511,254读）和211 MB（1,525,992读），分别纳入了最终组装中。 　　我们使用简短的读取软件Soapdenovo4（版本1014）生成了高质量的马铃薯基因组。我们首先将69.4 GB的Ga2配对短读读数组装成重叠群，它们是序列组件，而没有由重叠读数组成的差距。为了提高组装精度，仅考虑了具有高质量的读取的78.3％。然后将重叠群进一步链接到脚手架中，通过配对 - 末端关系（300至550 bp插入尺寸），伴侣读数（2至大约10 kb），fosmid末端（40 kb，90,407对末端序列）和BAC末端（100 kb，71,375最终序列）。然后，我们使用Illumina GA2读取生成的整个短阅读数据填补了空白。主要的重叠群N50大小（重叠长度使得使用相等或更长的重叠群产生一半的组装基因组碱基的一半）为697 bp，间隙填充后增加到1,318 kb（补充表23和24）。当在组装过程中仅使用配对端的关系时，N50支架的大小为22.4 kb。添加配偶对的读数为2、5和10 kb插入尺寸，N50脚手架尺寸分别增加到67、173和389 kb。当与较大插入尺寸的其他文库（例如Fosmid和BAC末端序列）集成时，N50达到1,318 kb。最终组装尺寸为727 MB，其中93.87％是非环绕序列。我们进一步填补了454个数据的6.74倍覆盖范围，这将N50重叠群的大小提高到31,429 bp，填补了15.4％。 　　组装的单基准确性是通过差异读数的深度和比例来估计的。对于DM v3.0组件，使用最佳参数，可以将8.8亿可用读数中的95.45％映射到组装基因组（参考文献34）中。计算了每个基因组位置的读取深度和整个基因组的峰值深度，CDS区域分别为100和105。大约96％的组装序列具有超过20倍的覆盖率（补充图1）。马铃薯基因组的总GC含量约为34.8％，GC含量和测序深度之间存在正相关（数据未显示）。DM马铃薯应该很少有杂合地点，而93.04％的站点可以由至少90％的读取支持，这表明基本质量和准确性很高。 　　基因型RH的全基因组测序使用各种片段大小在Illumina GA2平台上进行，并读取长度，从而产生144 GB的原始数据（补充表25）。使用自定义C程序过滤这些数据，并使用Soapdenovo 1.03（参考文献4）组装。此外，在Roche 454 Titanium Sequencer上测序了四个20-kb的序列库，总计为581 MB的原始数据（补充表26）。使用自定义Python脚本过滤所得的序列进行重复。 　　使用Sanger和454测序在各种覆盖范围（补充表9–11）中对RH BAC进行测序。使用与先前所述的类似方法（补充文本），使用自定义Python和C脚本校正了BAC序列中的共识碱基调用错误。使用BLAST 2.2.21（参考文献36）中的Megablast鉴定了同一物理瓷砖路径内BAC之间的序列重叠，并与Megamerger合并了Megamerger 6.1.0 package37。使用同一条管道，通过对齐初步的RH全基因组组件来封闭几个千个酶大小的间隙。通过使用Soapalign 2.20（参考文献34）映射RH全基因组Illumina和454个配合序列，通过映射RH全基因组光明和454个配合序列，然后使用自定义Python脚本处理这些映射结果，从而对所产生的非冗余重叠群进行脚打脚打。然后根据FPC图的平铺路径中的BAC顺序将脚手架订购为超交对管。该过程除去了25 MB的冗余序列，将序列片段的数量从17,228减少到3,768，并将N50序列长度从24 kb增加到144 kb（补充表9和10）。 　　为了沿着染色体锚定和完全定向物理重叠群，使用序列标记的位点（STS）标记，包括简单序列重复序列（SSR），SNP和多样性阵列技术（DART），从头开发了遗传图。SSR和SNP标记是直接从组装序列支架设计的，而多态性的飞镖标记序列被搜索到用于高质量独特匹配的支架上。总共有4,836个STS标记物，包括2,174个飞镖，2,304个SNP和358个SSR，分析了来自回击种群的180个后代克隆（（（DM×DI）×Di）在DM和DI之间开发的，DM和DI之间的CIP（CIPNo。703825），均为heterozygous diploID（cipno。703825）。Stenotomum SSP。使用JoinMap 4（参考文献38）分析了来自2,603个多态性STS标记的数据，其中包括1,881个飞镖，393个SNP和329个SSR等位基因，并产生了预期的12个马铃薯链接组。补充图3代表用7号染色体以7为例的马铃薯基因组的映射和锚定。 　　锚定DM基因组是使用直接和间接方法完成的。直接方法采用了（（DM×DI）×DI）链接图，其中2,603个STS标记中有2,037个由1,402个飞镖，376个SNP和259个SSR组成，可以独特地固定在DM SuperSpersCaffolds上。这种方法锚定了52％（394 MB）的组装，该组件排列成334个超交对车（补充表27和补充图3）。 　　RH是使用RHPOTKEY BAC Library 39的超高密度（UHD）链接MAP28的映射种群的男性父母。使用RH遗传和物理MAP28,40以及SGN（http://solgenomics.net/）的番茄遗传图数据，在硅锚固中利用的间接映射方法。来自RH遗传图的扩增片段长度多态性标记通过BLAST ARIGNMENT 36与DM序列支架相关联，该基因组序列序列TAGS41在RH物理图中从锚固的种子BAC中获得，或直接将序列的RH RHED SEED BAC与DM序列的直接比对。将组合标记的比对处理成可靠的锚点。使用SSAHA2将来自遗传图的番茄序列标记对准DM组装（参考文献42）。手动检查并纠正模棱两可的超级贩子的位置。这种方法锚定了另外32％的组件（229 MB）。在294个情况下，两种独立的方法提供了直接的支持，将相同的脚手架固定在两个地图上的同一位置上。 　　总体而言，这两种策略将649个超交对行为固定在覆盖623 MB长度的马铃薯遗传图上的近似位置。623 MB（86％）锚定的基因组包括39,031个预测基因中的90％。在未锚定的超夫妇中，在N90中发现了84个（622个脚手架大于0.25 MB），构成了17 MB的整体组装或组装基因组的2％。最长的SuperScaffold是7 MB（染色体1），最长的未经常规超级助手为2.5 MB。 　　在DNA和蛋白质水平上鉴定出马铃薯基因组组装中的转座元素（TES）。使用repbase43应用repotMasker29在DNA水平上进行TE识别。在蛋白质水平上，在针对TE蛋白数据库的WublastX36搜索中使用了重复蛋白毫米29,44，以进一步识别TES。对属于同一重复类的重叠TE进行了整理，如果序列重叠> 80％并且属于不同的重复类，则将其删除。 　　为了预测基因，我们对重复掩盖的基因组进行了从头算预测，然后将结果与蛋白质和转录本的剪接比对进行了使用，并使用GLEAN30进行了基因组序列。通过超过500 bp的重复序列掩盖了马铃薯基因组，除了通常在基因或内含子内部附近发现的微型倒置重复转座元件外，45。Augustus46和Genscan47软件用于从头开始预测，该参数训练了塔利亚纳。为了基于相似性的基因预测，我们将四个测序植物（A. thaliana，Carica Papaya，V。vinifera和Oryza sativa）的蛋白质序列对齐在马铃薯基因组上使用TBLASTN，具有1×10-5的E-VALUE临界值，然后使用相似的基因组序列，然后使用相似的基因组序列，用于使用匹配的蛋白质。在基于EST的预测中，使用Blat（身份≥0.95，覆盖范围≥0.90）将11个溶剂物种的EST序列对齐，以产生剪接的比对。所有这些资源和预测方法都是由GLEAN30组合在一起的，以建立共识基因集。为了最终确定基因集，我们将32个库中的RNA-seq对准了RNA-seq，其中8个库是单端和配对末端读数测序的，使用TopHat31将基因组对齐，然后使用默认参数将对齐方式用作CuffLinkS32的输入。过滤基因，转录本和肽组，以去除跨测序间隙，TE-编码基因以及其他不正确注释的小基因，基因。最终基因集包含39,031个基​​因，其中有56,218个蛋白质编码转录本，其中52,925种非相同蛋白被保留以进行分析。 　　从代表发育，非生物应激和生物应力条件的DM和RH的许多组织中分离RNA（补充表4和补充文本）。构建了cDNA库（Illumina），并在单个和/或配对模式下在Illumina GA2上进行测序。为了表示每个基因的表达，我们通过从每个基因中选择最长的CD从每个基因模型中选择了一个代表性转录本。对齐的读取数据是由TopHat31生成的，所选的转录本用作CuffLinkS32的输入，这是一种简短读取的转录物组件，该组件计算每种映射读数（FPKM）的片段作为每个转录物的表达值。袖扣以默认设置运行，最大内含子长度为15,000。报告了FPKM值并为每个转录本列表（补充表19）。 　　使用Orthomcl49使用11种植物物种的预测蛋白质组鉴定了寄生虫和直系同源簇（补充表28）。在删除未在早期注释步骤中未过滤的1,602个与TE相关的基因之后，使用初始的Orthomcl聚类鉴定了Asterid特异性和马铃薯特异性基因，然后使用BLAST搜索（E-VALUE截止）（1×10-5），针对plantgdb Project的ests和153的eStssberies（1×10-5）（nonaster nondersid）;57 Asterid物种；补充图5和补充表29。使用Intercoscan搜索标识的Interpro（http://www.ebi.ac.uk/interpro）对蛋白质结构域进行分析。我们将Asterid特异性和马铃薯特异性集的PFAM结构域与至少与其他一个非碱基基因组或转录组共享的域进行了比较。Fisher的精确测试用于检测蛋白质组之间的PFAM表示存在显着差异。 　　删除自我和多次匹配后，使用MCSCAN9和I-Adphore 3.0（参考文献50）鉴定了同步块（≥5个基因），基于对齐蛋白基因对（补充表8）。对于自我对齐的结果，每个对齐块代表源自基因组重复的寄生段对，而对于种间比对结果，每个比对块代表源自共享祖先的直系同源对。我们从对齐块中计算了每个基因对的4DTV（第三个密码子的四倍变性同义位点），并给出4DTV值的分布，以估算进化历史记录中发生的物种或WGD事件。 　　使用HMMER v.3（http://hmmer.janelia.org/software）筛选了Tuberosum组Phureja组装V3的注释中预测的公开阅读框架（ORF）。HMM从PFAM主页（http://pfam.sanger.ac.uk/）下载。使用NBS域的原始HMM进行分析，导致351个候选者。从中，高质量的蛋白质集（<1 × 10−60) was aligned and used to construct a potato-specific NBS HMM using the module ‘hmmbuild’. Using this new potato-specific model, we identified 500 NBS-candidate proteins that were individually analysed. To detect TIR and LRR domains, Pfam HMM searches were used. The raw TIR HMM (PF01582) and LRR 1 HMM (PF00560) were downloaded and compared against the two sets of NBS-encoding amino acid sequences using HMMER V3. Both TIR and LRR domains were validated using NCBI conserved domains and multiple expectation maximization for motif elicitation (MEME)51. In the case of LRRs, MEME was also useful to detect the number of repeats of this particular domain in the protein. As previously reported52, Pfam analysis could not identify the CC motif in the N-terminal region. CC domains were thus analysed using the MARCOIL53 program with a threshold probability of 90 (ref. 52) and double-checked using paircoil2 (ref. 54) with a P-score cut-off of 0.025 (ref. 55). Selected genes (±1.5 kb) were searched using BLASTX against a reference R-gene set56 to find a well-characterized homologue. The reference set was used to select and annotate as pseudogenes those peptides that had large deletions, insertions, frameshift mutations, or premature stop codons. DNA and protein comparisons were used. 　　RH reads generated by the Illumina GA2 were mapped onto the DM genome assembly using SOAP2.20 (ref. 34) allowing at most four mismatches and SNPs were called using SOAPsnp. Q20 was used to filter the SNPs owing to sequencing errors. To exclude SNP calling errors caused by incorrect alignments, we excluded adjacent SNPs separated by <5 bp. SOAPindel was used to detect the indels between DM and RH. Only indels supported by more than three uniquely mapped reads were retained. Owing to the heterozygosity of RH, the SNPs and indels were classified into heterozygous and homozygous SNPs or indels. 　　On the basis of the annotated genes in the DM genome assembly, we extracted the SNPs located at coding regions and stop codons. If a homozygous SNP in RH within a coding region induced a premature stop codon, we defined the gene harbouring this SNP as a homozygous premature stop gene in RH. If the SNP inducing a premature stop codon was heterozygous, the gene harbouring this SNP was considered a heterozygous premature stop codon gene in RH. In addition, both categories can be further divided into premature stop codons shared with DM or not shared with DM. As a result, the numbers of premature stop codons are 606 homozygous PS genes in RH, 1,760 heterozygous PS genes in RH but not shared with DM, 288 PS in DM only, and 652 heterozygous premature stop codons in RH and shared by DM. 　　To identify genes with frameshift mutations in RH, we identified all the genes containing indels of which the length could not be divided by 3. We found 80 genes with frameshift mutations, of which 31 were heterozygous and 49 were homozygous. 　　To identify DM-specific genes, we mapped all the RH Illumina GA2 reads to the DM genome assembly. If the gene was not mapped to any RH read, it was considered a DM-specific gene. We identified 35 DM-specific genes, 11 of which are supported by similarity to entries in the KEGG database57. To identify RH-specific genes, we assembled the RH Illumina GA2 reads that did not map to the DM genome into RH-specific scaffolds. Then, these scaffolds were annotated using the same strategy as for DM. To exclude contamination, we aligned the CDS sequences against the protein set of bacteria with the E-value cut-off of 1 × 10−5 using Blastx. CDS sequences with >90% identity and >90％的覆盖范围被认为是污染物，被排除在外。另外，将所有DM RNA-SEQ读数映射到CDS序列上，并且排除具有同源读数的CDS序列，因为这些基因可能是由于组装不正确所致。我们总共预测了246个RH特异性基因，其中34个由基因本体学注释支持17。