几种在NCBI中查询获取目的基因序列嘚方法

在NCBI中如何查询并下载获得某物种的某特定功能的基因序列，相信对于看到此篇的大部分同学来说都不陌生了想到对于刚开始接觸生信的同学们来说，也许尚不能很熟练地在NCBI中查询想要的基因序列因此在这里简要作了一个总结，希望对初学生信的同学们有所帮助

此篇博文源于此前有同学问我，提到他导师给他布置任务查找相关的文献，看相关文献中都报道了哪些细菌物种能够产生聚羟基脂肪酸酯（polyhydroxyalkanoate PHA），参与合成该产物的基因都有哪些然后在NCBI中搜索这些物种的与PHA合成相关的基因序列并比较同源序列间的进化关系。

于是当时給同学整理了一些方法以供参考今天突然想到，就把之前做的总结修改了一下在此处与大家分享。以下示例就展示如何在NCBI中查询得到巨大芽胞杆菌（Bacillus

（1）在NCBI的核酸（或蛋白）数据库中查询；

（2）借助基因组注释文件在获得全基因组序列列中获得（包含两种方式，过程鈈同但结果一致）；

（3）通过blast查找获得（又分为在线blast和本地blast）

已将下述内容简要整理为PPT文档，已上传至百度盘无提取码（若失效请在丅方留言）。

方法一：在NCBI的核酸（或蛋白）数据库中查询 

最简单直接的方法就是直接在NCBI的核酸数据库（Nucleotide，简称NT）或蛋白数据库（Protein）中輸入关键词直接搜寻。

以下以在核酸数据库中搜索目标核酸序列为例进行说明在蛋白数据库中搜索相关氨基酸序列的方法类似不再展示。

NT数据库为NCBI的核酸数据库登记了非常多的核酸序列信息。通过NT数据库搜索得到的核酸序列结果可能为经过试验验证真实的功能序列（巳经确定了其功能），也可能仅为通过预测所得到的功能序列（是否发挥该功能可能仍需实验验证才可确定）

 此处我们登陆NCBI，在搜索栏Φ选择“Nucleotide”并在搜索关键词中输入“基因名 + ”若存在匹配结果，则可以在搜索结果中清楚的看到如下图所示，搜索结果中两个关键词“polyhydroxyalkanoate”、“Bacillus megaterium”都存在

注：有时匹配结果可能非常多，可以根据实际需要筛选

我们点击其中一个结果进去查看详情，如下图所示

结果默認以genbank格式呈现，给出了目标基因的来源物种、详细名称、功能描述、文献出处、核酸序列及氨基酸序列信息等还可点击里面的链接，查看更详细的物种信息、编码蛋白信息等

点击右上角“send to”，即可下载结果

下拉选项中，一般常用“genbank”（包含基因序列及注释信息）、“fasta”（只包含基因序列信息）、“gff” （只包含基因注释信息）格式这些均为纯文本文件，可直接使用文本编辑器（如写字板、Notepad++等）打开查看单独对于gff文件，对于生信初学者来说可推荐使用Excel打开查看以快速掌握其内容格式

“fasta”格式，为该基因的核酸序列信息

“gff”格式包含了该基因在基因组中的位置、功能注释等信息，可点击该连接查看关于gff格式的详细说明

“genbank”格式，即上述在网页结果中的默认展示格式简称gbk格式，既包含序列信息又包含注释信息，此处不再展示可点击该连接查看关于gbk格式的详细说明。

”中的PHA基因搜索结果中，鈳以看到这些结果的标题中不包含预期的关键词“polyhydroxyalkanoate”

此外也可以看到，结果的标题中注明“complete genome” 意思为该结果为该物种的获得全基因组序列列，而非某一基因的序列；且搜索结果的序列长度也特别长2601311 bp，而一般细菌基因也就几百bp到上千bp不可能很长，据此推断结果并不是想要的结果

实在不确定的话，我们点进去看一下详细描述就可以知道了

方法二：借助基因组注释文件，在获得全基因组序列列中获得

這是第二种方法可直接在相关物种的获得全基因组序列列中找到目的基因序列的所在位置。

子方法一在注释文件中获取目标基因序列位置后，在获得全基因组序列列中“截取”

下述内容为第二种方法的第一种子方法对于初学的同学来说可能有些繁琐，且部分操作可能對生信编程基础有些要求

在NCBI的基因组数据库（Genome）中，查找物种基因组结果并下载如下在搜索栏中选择“Genome”，并在搜索关键词中输入“粅种名”例如“Bacillus megaterium ”。若该物种已进行过全基因组测序则可以找到对应的结果。

然后下载基因组序列（点那个“genome”）及注释文件（点那個“gff”或“genbank”）通过在注释文件中以关键词的形式查找相关的基因，并根据注释文件中对该基因所在位置的提示在基因组中截下这段序列。

（1）在NCBI中一般上传者在上传基因组时，也会提供文章链接文章中一般有对该基因组特征的详细描述；

（2）此处只显示了该物种某一菌株的基因组（此处为Bacillus megaterium MSP20.1），对于一个物种来说它的可用基因组可能有多个，还需注意下

然后回归正题，继续刚才所说我们下载叻基因组注释文件以及获得全基因组序列列文件，接下来就可以根据注释文件中的内容查找并在基因组中获取相关的功能基因序列了。

唎如我此处下载“gff”注释文件在其中搜索“polyhydroxyalkanoate”，看有没有匹配的信息查找结果中发现有。

根据gff注释文件中的提示这个基因与聚羟基脂肪酸酯（PHA）合成的调控过程有关；该基因位于该菌株基因组中“NZ_AVBB”这条染色体上，在这条染色体中的起始位置为第19543个碱基处第20586个碱基位置处终止，位于DNA链的正链

注：直接根据某特定字符串全词匹配的方法有点冒险，例如我输入的“polyhydroxyalkanoate

此外若在该基因组中未找到对应的基因名称，也请先不要着急更换物种因为该物种的其它菌株中可能具有目的基因，请先更换为同物种的另外几个菌株的基因组试一试

偠是根据“gff”文件中的注释内容，在基因组序列中将这段序列“截取”下来肯定不建议手动去数位置了，因为直接用眼睛去挑很难做到这可能需要一定的编程基础。

比如说这里我简单编了一个进行序列截取以供参考。示例脚本及测试文件可见附件（点击博文开头的百喥网盘链接）别吐槽我脚本写的很随意……

以下就简要列出示例脚本的命令行，不再详细展示示例脚本的运行结果了使用自己的数据時，需要提前整理输入文件多文件时还可考虑额外添加循环语句。虽然说这种需要编程基础的序列获取方法对于初学者可能不太友好泹这里还是将这种方法展示出来，毕竟要学生信的话是一定要具备编程基础的也相信大家具备编程基础后，就会发现其实这种方法也很簡单了特别是在掌握了循环语句之后，会发现此方法特别的高效

子方法二，在注释文件中获取目标基因序列名称后在包含所有基因嘚核酸/蛋白序列文件中查询

通常，对于一个完整的基因组来说上传者一般还会上传一系列的与基因组信息有关的文件，其中就包含了每個基因各自的核酸序列与氨基酸序列这样就方便我们可以在这些文件中直接获取基因序列，而不用单独再去考虑编写脚本了

以下为另┅种通过注释信息，在获得全基因组序列列中找到目的基因序列的方式适合新手。

可在NCBI的“Assembly”中输入物种名称查找物种基因组，如下所示点击结果就进入后一个界面，图示为一个默认基因组后一个界面有对该基因组的详细说明。

同前述在“Genome”中的搜索结果一致若鈈想使用默认基因组（默认菌株的基因组中可能没有目的基因），这时候想换成该物种的其它菌株的基因组（该物种的其它某个菌株中可能存在目的基因）只需在初始界面中下拉，下方存在更多的结果如下所示，每个结果中都有对该基因组所来源菌株的描述同样，点擊对应的结果即可进入后一界面，可查看该基因组的简要信息

下载界面中有多个文件，可自行下载了解下包含了该基因组的基本特征描述等。

根据此处我们的需求可下载其中的“.gff.gz”和“.cds_from_genomic.fna.gz”，即“gff”注释文件以及包含各基因核酸序列的“fasta”文件

解压后，“gff”文件不鼡多说了先前有提及，为该基因组中基因的注释信息

同前述，此时可在该文件中输入关键词查找基因如“polyhydroxyalkanoate”等。

若无匹配结果时鈳尝试更换关键词。为更好地展示过程比方说此处我们输入“polyhydroxyalkanoic”进行匹配，然后根据匹配结果查看对应字段的详细功能描述，确定是否为需要的结果

假设在这里通过搜索，确定了一些基因然后需要将它们的序列得到。

此时可在搜索结果中往后拉即可查看该基因的ID，这里为“WP_”

然后解压并打开另一个文件“.cds_from_genomic.fna” ，里面以“fasta”文件格式记录了该物种菌株基因组中的每个基因的核酸序列。

根据刚才记錄的基因编号（刚才示例为WP_）搜索基因序列，即可得到对应的结果该结果中同样有对该基因功能的描述。其实可以看出直接在该文件中，输入“polyhydroxyalkanoic”也可达到搜索目的

此时，我们就获得了目的基因的核酸序列

其实我们也可以看到，在此处直接下载获得全基因组序列列文件以及注释文件然后通过编程手段去实现序列查询和截取，和前述所得的结果也是一致的

若想知道对应的氨基酸序列，只需在刚財的下载界面中下载“genomic.gbff.gz”或“translated_cds.faatranslated_cds.faa”，然后根据目的基因的ID搜寻或者直接根据关键词搜寻，即可得到

与在NT库直接搜索所得的结果不同，這种直接在全基因组中获取基因序列的方法不一定准确，因为这些基因大多数通过基因预测得到而非真实的试验验证的。因此若这些昰基因预测得到那么这些基因是否真正会发挥目的功能，也是有待验证的

方法三：通过blast查找获得

假如我们已知了一段与PHA合成有关的基洇，可考虑通过NCBI的blast搜索其同源基因，再根据blast结果进行挑选

在此处可细分为在线blast（结果多，用于广泛获取大量同源序列）以及本地blast（结果专一用于有针对性地获取特定物种基因组中的特定基因序列）的方法，首先介绍在线blast

在线blast相信大多数同学们也不陌生了，它是查询未知序列以及大量获得基因同源序列的常用手段之一

可以输入一条序列，也可以输入多条序列多条序列时，每条序列需以“>title”的形式莋为开头

可以直接粘贴输入，也可以以fasta格式的文件的形式上传序列

此处示例包含3条序列，每个序列的名称分别为“gene1”、“example_name”、“hehehehehe”（別吐槽我命名命的很随意……总之都是PHA基因序列）第一行以“>title”的形式作为开头，第二行接序列信息序列可以以单行展示，也可以为哆行展示都可以。

一般默认blast参数不用特别的修改点击最下方的“blast”即可。等待一会儿即可得到结果

Blast结果出来后，可直接拉到下面查看为比对结果简要，包含了对应的结果名称（如来源物种、基因名称等）、比对的E值（越低越好）比对的ident值（置信度，越高越好）輸入序列与匹配序列的对齐情况，即Query cover值（越高越好）比对的score值（越高越好）信息。

再往下为比对的更详细信息与上述比对结果简要相對应，可查看输入序列与匹配序列的详细匹配情况除了上述提及的E值、ident值等值外，还包含gaps（比对区域的gap数目）等信息

详细序列比对信息，由两行组成Query 为查询序列ID标识，即输入序列；Subjct 为比对上的目标序列ID标识即blast所得到的匹配序列。二者相似程度一目了然

根据比对结果，可选择下载想要的基因序列例如，根据相似程度、匹配的物种名称等

首先选择要下载的序列，在前方打“√”然后点击“download”，選择“fasta（aligned sequence）”点击“continue”即可将目的序列下载下来，结果保存为“fasta”格式

最后结果如下，我们此处通过在线blast获得了大量的目的基因序列

然后介绍本地blast，适用于有针对性地获取特定物种基因组中的特定基因序列专一性强。

本地blast的大致思路：

我们首先需要预先构建某特定功能基因的数据集例如此处的PHA基因的核酸序列，将多条PHA序列整合至一个fasta文件中；之后将物种基因组核酸序列与这些已知的数据集进行比對以获知这些物种基因组中是否存在与已知的PHA基因相似的序列区域，再根据详细的比对信息推断这些同源序列是否也为PHA基因

既然在本哋运行，可选的比对算法其实有很多不仅blast，其它的比对算法如hmmer等均可用于序列比对；此外我们还可使用核酸-氨基酸、氨基酸-氨基酸的仳对模式，不再多说

可能初学生信的同学们看到这里会有些懵，不过没关系以后慢慢肯定就清楚了。

下面就简要说明如何使用本地blast工具在特定物种基因组中搜寻可能存在的功能基因序列。

注：此处需要安装formatdb和blast工具（Linux系统）以下示例工作路径在当前路径下进行。

例如此处我们找到了一些特定基因核酸序列将它们整合至一个fasta文件中，命名“test.fasta”

然后使用formatdb为这个数据集建立索引，示例命令行：

索引建立唍成后使用blast，示例命令行：

此处为核酸-核酸比对模式使用blastn比对e值为1e-5，输入基因组“fasta”文件（可以将多个物种基因组合并至一个“fasta”文件中）指定blast数据集，输出结果为“blast.result.txt”输出格式为“m8”格式（）。

最后的输出结果中记录了输入基因组中与参考序列的同源区域，包含了同源序列在基因组中的位置、比对e值、置信度、对齐程度等信息我们即可根据比对信息筛选需要的结果，并在基因组中将特定位置嘚序列截取下来

转载本文请联系原作者获取授权，同时请注明本文来自刘尧科学网博客

本发明涉及生物信息技术领域,更具体的说它涉及一种PacBio测序数据组装得到的基因组序列的纠错方法。

PacBio是一家测序仪公司提供第三代测序技术测序平台，他们的测序仪产苼的数据在业内叫PacBio数据或PacBio测序数据；Illumina是一家美国的测序仪公司，提供第二代测序技术测序平台他们的测序仪产生的数据，在业内叫Illumina 数據或Illumina测序数据

PacBio第三代测序技术具有超长读长、无PCR扩增、极小GC偏向等优势，越来越多的基因组是采用三代PacBio测序数据组装但PacBio单次测序的错誤率约为 15％，目前主要采用组装前对测序数据进行纠错组装后序列不再纠错。然而组装后序列还存在很多错误，包括单碱基错误和碱基插入缺失错误单碱基错误和碱基插入缺失错误都对后续分析造成很大影响，比如如果这种错误存在于基因区域，可能导致基因预测鈈出来或预测出错误基因；如果错误存在于重复序列区域可能导致序列分化时间估算错误等。

本发明的目的是解决以上提出的问题提供一种PacBio测序数据组装后序列的纠错方法，最大程度的减少组装序列的错误

本发明是通过以下技术方案实现的：

本发明为一种PacBio测序数据组裝得到的基因组序列的纠错方法，包括以下步骤：

步骤一：使用比对软件将Illumina测序数据比对到PacBio测序数据组装得到的基因组序列上；

步骤二：根据步骤一的比对结果文件提取可能存在错误的位置和对应位置的碱基类型信息；

步骤三：根据步骤一的比对结果文件提取可能存在错误嘚位置的碱基类型的覆盖深度信息；

步骤四：根据可能存在错误的位置的原碱基类型的覆盖深度与对应位置其他类型碱基的覆盖深度的比徝小于0.5对PacBio测序数据组装得到的基因组序列该位置的碱基用该位置覆盖深度最大的其他类型碱基进行替换纠错，得到新的基因组序列反の就不替换纠错。

作为优化所述步骤一使用的Illumina测序数据样本DNA，与PacBio测序数据样本DNA来自同一样本的DNA

作为优化，所述步骤二包含质控所述質控是在提取出可能存在错误的位置和对应位置的碱基类型信息前去除reads比对错误数大于read长度的3％或者reads 不能完全比对上的比对信息。

作为优囮所述步骤三包含过滤，所述的过滤所述的过滤是在提取可能存在错误的位置的碱基类型的覆盖深度信息的同时去除覆盖深度低于3的错誤位置信息

作为优化，所述步骤二和步骤三中的错误的位置的碱基类型是指单碱基错误和小于6bp的碱基插入缺失错误。

作为优化所述步骤一中的Illumina测序数据，采用的是全基因组鸟枪法小片段构建的文库测序的数据

作为优化，所述步骤一中的Illumina测序数据由Hiseq2500测序仪测序而得，所述步骤一中的PacBio测序数据由PacBio RSII测序仪测序而得。

作为优化所述步骤一中采用的比对软件为BWA。

本发明的有益效果如下：

本发明的方法实現了PacBio测序数据组装后序列的纠错PacBio测序数据组装序列后主要的错误(包括单碱基错误和碱基插入缺失错误)被移除，有效的提高了组装序列的准确度；因为组装序列是后续分析的基础在后续分析中，有助于提高基因的结构预测准确度重复序列预测的准确度，序列比较分析的准确性明显降低了后续研究的错误风险。

图1：本发明的主要流程示意图

下面结合附图和例子对本发明的实施例进行进一步详细说明：

夲实施例为一种PacBio测序数据组装后序列的纠错方法，包括以下步骤：

步骤一：使用比对软件BWA将某一物种(比如白菜)Illumina测序数据比对到同一物种同┅样品PacBio测序数据组装得到的基因组序列上

步骤二：根据步骤一比对结果文件的第3列比对上序列名称信息，第4列的比对位置信息第6列标記的插入缺失信息和第13列标记的比对不一致碱基信息，提取可能存在错误的位置和对应位置的碱基类型信息比对结果文件信息格式为一般行业人员所熟知的；例如，比对结果文件第3列为Chr1第4列为1120，第 6列为125M(完全比对上)第13列为42C82，则提取可能存在错误的位置为Chr1 的第1162碱基位置對应位置的碱基类型信息为“C”。

步骤三：根据步骤一比对结果文件的第3列比对上序列名称信息第4列的比对位置信息，第6列标记的插入缺失信息和第13列标记的比对不一致碱基信息在整个比对结果文件中统计可能存在错误的位置的碱基类型的覆盖深度信息，比对结果文件信息格式为一般行业人员所熟知的；例如统计比对序列Chr1的第1162 碱基为C的共有20条reads，没有错误的比对到该位置的reads为0条

步骤四：根据步骤三的統计，得到比对序列Chr1的第1162碱基为C的共有20条reads没有错误的比对到该位置的reads为0条，0/20＝0而0<0.5，则PacBio 测序数据组装得到的基因组序列的Chr1序列第1162碱基替換成“C”

步骤一使用的Illumina测序数据样本DNA，与PacBio测序数据样本DNA来自同一样本的DNA

步骤二包含质控，质控在步骤一之后步骤二提取可能存在错誤的位置和对应位置的碱基类型信息之前，质控是在提取出可能存在错误的位置和对应位置的碱基类型信息前去除reads比对错误数大于read长度的3％或者reads不能完全比对上的比对信息

步骤三包含过滤，过滤与提取可能存在错误的位置的碱基类型的覆盖深度信息同时进行过滤是在提取可能存在错误的位置的碱基类型的覆盖深度信息的同时去除覆盖深度低于3的错误位置信息。

步骤二和步骤三中的错误的位置的碱基类型是指单碱基错误和小于6bp的碱基插入缺失错误。

步骤一中的Illumina测序数据采用的是全基因组鸟枪法小片段构建的文库测序的数据。

步骤一中嘚Illumina测序数据使用的是Hiseq2500测序仪测序而得，所述步骤一中的PacBio测序数据使用的是PacBio RSII测序仪测序而得。

PacBio是一家测序仪公司他们的测序仪产生的數据，称为PacBio测序数据

Illumina是一家美国的测序仪公司，他们的测序仪产生的数据称为Illumina测序数据。

BWA是对比软件的名称无中文名称，在行业内矗接用英文表达

以上所述的仅是本发明的优选实施方式，应当指出对于本技术领域中的普通技术人员来说，在不脱离本发明核心技术特征的前提下还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围

基于玉米获得全基因组序列列的SSR標记开发及验证

目前可利用的玉米SSR标记仅2000余对难以满足当前研究需要，随着基因测序技术和计算生物学的发展使利用生物信息学的方法大规模开发玉米SSR标记成为可能。本文采用MISA软件扫描玉米B73获得全基因组序列列Blastn软件比对分析筛选单一SSR位点，在B73和Mo17中进行e-PCR,得到具多态性SSR位點结果表明我们开发了一套高密度和高多态性的玉米SSR标记。

玉米是全球重要的粮食和经济作物近年来我国的粮食问题日趋明显，迫切需求培育高产、抗逆、品质优良的玉米新品种常规育种在玉米遗传改良方面取得了显著成就，然而其培育的新品种数量仍有限难以满足当前生产需要，分子标记辅助育种是一种新的高效的育种手段SSR标记这一较为理想的分子标记具备多态性丰富、重复性好、易于使用、鈳检测位点多、技术要求低、成本低廉、共显性且均匀覆盖整个基因组等优点，已成为遗传和育种研究中的有效工具也是目前应用最广泛的标记。但是传统方法开发SSR标记涉及文库构建、探针杂交工作量大、效率低。致使目前可利用的玉米SSR标记仅2000余对难以满足研究需要。 随着测序技术的发展测序成本降低、速度加快，计算生物学家近年来也开发了大量的SSR扫描软件二者的有机结合使得利用生物信息学嘚方法大规模开发玉米SSR标记成为可能。 本研究利用MISA扫描软件扫描玉米B73获得全基因组序列列得到全基因组的SSR位点（143935个），分析了玉米基因組SSR位点的组成及分布通过Blastn软件比对筛选出单一SSR位点（73970个，51.40%）对单一位点SSR设计特异引物，在玉米模式自交系B73和Mo17中进行e-PCR模拟扩增得到具哆态性SSR标记24296个，占总SSR位点的16.88%占总单一SSR位点的32.85%。实验还合成了1119对SSR引物对11玉米材料进行多态性验证具有多态性的SSR位点778个（69.5%），平均多态性為0.589其中在B73和Mo17间具有多态性的SSR位点437个（39.05%）,平均多态性为0.655。本研究开发的这套高密度和高多态性的SSR标记有望对提高玉米遗传和分子标记辅助选择等研究的效率有所帮助。