搜索
查看: 3992|回复: 1

【学习笔记】从Nr库中任意提取各分类水平的子库

[复制链接]

58

主题

103

帖子

756

积分

版主

Rank: 7Rank: 7Rank: 7

积分
756
QQ
发表于 2016-11-29 21:50:51 | 显示全部楼层 |阅读模式
本帖最后由 Panda姐 于 2016-11-30 20:07 编辑


此文来自我的博客:从Nr库中任意提取各分类水平的子库

写在前面
[size=15.5556px]月初的时候我就在研究这件事,在我坚持不懈的寻找下,终于在github上找到了跟我有一样想法的现成代码啦!但是拥有跟前辈一样的思考思路,我还是挺开心的。从Nr库中任意提取各分类水平的子库,比如蓝藻库啊,真核藻库啊,病毒库啊,细菌库啊~而且还可以进行完整基因组的下载等等[size=15.5556px]~
[size=15.5556px]
[size=15.5556px]下面是我运行过程中使用的代码记录:

#########################################################################
安装
[Shell] 纯文本查看 复制代码
wget -c [url=ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zip]ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zip[/url]
wget -c [url=ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz]ftp://ftp.ncbi.nih.gov/pub/taxon ... .accession2taxid.gz[/url]
wget -c [url=ftp://ftp.ncbi.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_refseq.txt]ftp://ftp.ncbi.nih.gov/genomes/A ... _summary_refseq.txt[/url]
unzip taxdmp.zip
gzip -d prot.accession2taxid.gz
  • 建关联库:ncbi_taxonomy.db
    耗时较长,耐心等待~~

[Shell] 纯文本查看 复制代码
create table acc_taxid (accession text, accession_version text, taxid integer, gi integer);
.tables
.mode list
.separator \t
.import prot.accession2taxid acc_taxid # 导入prot.accession2taxid数据
CREATE UNIQUE INDEX accvers_idx_on_accvers_taxid ON acc_taxid(accession_version);
# 唯一索引 使用唯一索引不仅是为了性能,同时也为了数据的完整性。唯一索引不允许任何重复的值插入到表中。
CREATE INDEX taxid_idx_on_accvers_taxid ON acc_taxid(taxid); # 单列索引
# CREATE INDEX语句创建索引,它允许命名索引,指定表及要索引的一列或多列,并指示索引是升序排列还是降序排列。
.exit

  • 检查关联库的完整性

[Shell] 纯文本查看 复制代码
sqlite3 ncbi_taxonomy.db
SELECT taxid FROM tree WHERE name = "Acidobacterium capsulatum";
.exit
[size=15.5556px]检验成功啦~~
[Shell] 纯文本查看 复制代码
## 将taxomia链接到python安装目录下,可以让python找到taxomia
[shenmy@localhost taxomias]
ln -s /home/shenmy/biosoft/taxomias /home/shenmy/anaconda2/lib/python2.7
chmod +x setup_taxomias.py taxomias.py
测试
  • 从NR库中分出微囊藻属的非冗余蛋白acc号列表

[Shell] 纯文本查看 复制代码
## 下载Nr库
wget [url=ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz]ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz[/url]
gzip -d nr.gz
[Python] 纯文本查看 复制代码
# 打开python运行以下代码:
#载入所需要的包
import sqlite3,taxomias
# 看看微囊藻的分类号(taxid)是多少?答案应该是1125。
print taxomias.TaxidByName("Microcystis")
# 把所有属于微囊藻的accession version 号存入一个set,set比list会更快些,因为集合是使用哈希列表。
wanted= set(taxomias.AllAccByTaxid(1125))
acc=open('nr_Microcystis_acc_list.txt', 'a')
acc.write('\n'.join(wanted))
  • 从fasta文件中提取id列表中对应的序列
    可以自己编写脚本,我采用现成的工具(可能不太对,还是要自己写脚本~~):

[Shell] 纯文本查看 复制代码
wget -c [url=http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faSomeRecords]http://hgdownload.cse.ucsc.edu/a ... 86_64/faSomeRecords[/url]
chmod +x faSomeRecords
faSomeRecords /home/shenmy/Metavirome/RefSeq/nr nr_Microcystis_acc_list.txt nr_Microcystis.fa

Tips:
从安装到测试完成花了我6个小时左右(解决中间各类报错,读懂运行的所有python代码,修补命令);提取的微囊藻的nr蛋白子库的正确性还没有验证,但是比我之前发布的方法靠谱,之前的文章,回复“2016-11-02”即可获得。



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:线粒体研究小结
下一篇:使用DESeq2寻找差异表达的基因中result函数
回复

使用道具 举报

1

主题

11

帖子

132

积分

注册会员

Rank: 2

积分
132
发表于 2017-4-11 15:49:21 | 显示全部楼层
2016-11-02
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-21 05:35 , Processed in 0.033750 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.