lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

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

已有 20295 次阅读 2019-1-20 00:18 |个人分类:经验分享|系统分类:科研笔记| _信息检索, _数据库

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

 

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

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

于是当时给同学整理了一些方法以供参考,今天突然想到,就把之前做的总结修改了一下,在此处与大家分享。以下示例就展示如何在NCBI中查询得到巨大芽胞杆菌(Bacillus megaterium)这个物种中与聚羟基脂肪酸酯(polyhydroxyalkanoate PHA)合成相关的基因。

此处共介绍3种方法:

1)在NCBI的核酸(或蛋白)数据库中查询;

2)借助基因组注释文件,在全基因组序列中获得(包含两种方式,过程不同但结果一致);

3)通过blast查找获得(又分为在线blast和本地blast)。


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

https://pan.baidu.com/s/1NC674WW_s9WH4ZbOKycoYg



方法一:在NCBI的核酸(或蛋白)数据库中查询 


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

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


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

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

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

1.png

 

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

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

2.png

 

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

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

3.png


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

4.png

 

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

https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md

5.png

 

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

https://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html

 

有时在搜索时匹配不到结果,如下所示,查找了物种“Corynebacterium hydrocarboxydans ”中的PHA基因,搜索结果中,可以看到这些结果的标题中不包含预期的关键词“polyhydroxyalkanoate”

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

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

6.png

 

若未找到想要的结果时,请先不要着急,这时候可以考虑换其它的关键词搜索试试,如使用“pha”“phaA”“hydroxybutyrate”“polyhydroxyalkanoic”“poly-beta-hydroxybutyrate”等,或许可以搜索到。

 


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


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


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


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

 

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

7.png

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

 

补充说明:

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

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


8.png

 

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

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

10.png

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

注:直接根据某特定字符串全词匹配的方法有点冒险,例如我输入的“polyhydroxyalkanoate ”,但有时候“gff”文件中的注释的结果可能为“polyhydroxyalkanoic”之类的,会导致匹配不到这种结果。同前述,这时候可以考虑换其它的关键词搜索试试,如使用“pha”“phaA”“PHA”“hydroxybutyrate”“polyhydroxyalkanoic”“poly-beta-hydroxybutyrate”等。

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

 

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

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

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

示例命令行:

python3 seq_select.py –i test.fasta –l test.list.txt –o result.fasta

以下为“python3 seq_select.py -h”输出。

11.png

 

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


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

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

 

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

12.png

 

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

13.png

 

然后点击右方的“Download the RefSeq assembly”,可进入下载界面。

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

14.png

 

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

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

15.png

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

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

16.png

 

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

此时可在搜索结果中往后拉,即可查看该基因的ID,这里为“WP_013055934.1”。

17.png

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

18.png

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

19.png

 

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

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

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

 

补充说明:

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

 


方法三:通过blast查找获得



在线blast


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

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

 

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

首先进入NCBI在线blast界面:https://blast.ncbi.nlm.nih.gov/Blast.cgi

若输入序列为核酸序列,就用“Nucleotide BLAST”;若为氨基酸序列,就用 Protein BLAST”。

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

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

 

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

20.png

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

21.png

 

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

22.png

 

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

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

23.png

 

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

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

24.png

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

25.png

 

本地blast


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

本地blast的大致思路:

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

 

补充说明:

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

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

 

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

注:此处需要安装formatdb和blast工具Linux系统),以下示例工作路径在当前路径下进行

 

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

26.png

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

formatdb -i test.fasta -p F -o T

索引建立完成后,使用blast,示例命令行:

blastall -p blastn -e 1e-5 -m 8 -i genome.fasta -d test.fasta -o blast.result.txt

此处为核酸-核酸比对模式使用blastn,比对e值为1e-5,输入基因组“fasta”文件(可以将多个物种基因组合并至一个“fasta”文件中),指定blast数据集,输出结果为“blast.result.txt”,输出格式为“m8”格式(http://www.drive5.com/usearch/manual/blast6out.html)。

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

 




http://blog.sciencenet.cn/blog-3406804-1158033.html

上一篇:关于在ggplot2中绘制截断坐标轴的方法
下一篇:R包randomForest的随机森林方法鉴定组间代表物种并进行样本分类

2 陈洋 孙雪

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 18:53

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部