R语言 NCBI蛋白数据库分库的方法

今天小编给大家分享一下R语言 NCBI蛋白数据库分库的方法的相关知识点，内容详细，逻辑清晰，相信大部分人都还太了解这方面的知识，所以分享这篇文章给大家参考一下，希望大家阅读完这篇文章后有所收获，下面我们一起来了解一下吧。

NCBI 蛋白数据库分库

1.下载数据：

#wget -c https://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz #wget -c https://ftp.ncbi.nlm.nih.gov/genbank/livelists/gi2acc_mapping/gi2acc_lmdb.db.gz #wget -c https://ftp.ncbi.nlm.nih.gov/genbank/livelists/gi2acc_mapping/gi2accession.py #wget -c https://ftp.ncbi.nih.gov/blast/db/FASTA/nr.gz #wget -c https://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz #快速下载方法 /share/work/biosoft/aspera/latest/cli/bin/ascp -v -k 1 -T -l 400m -i ~/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/genbank/livelists/gi2acc_mapping/gi2acc_lmdb.db.gz ./ /share/work/biosoft/aspera/latest/cli/bin/ascp -v -k 1 -T -l 400m -i ~/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/pub/taxonomy/accession2taxid/prot.accession2taxid.gz ./ /share/work/biosoft/aspera/latest/cli/bin/ascp -v -k 1 -T -l 400m -i ~/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nr.gz ./

2.写脚本根据 prot.accession2taxid.gz　文件中的蛋白ID，分类信息对所有的蛋白序列ｎｒ库进行分类

0       |       BCT     |       Bacteria        |               | 1       |       INV     |       Invertebrates   |               | 2       |       MAM     |       Mammals |               | 3       |       PHG     |       Phages  |               | 4       |       PLN     |       Plants and Fungi        |               | 5       |       PRI     |       Primates        |               | 6       |       ROD     |       Rodents |               | 7       |       SYN     |       Synthetic and Chimeric  |               | 8       |       UNA     |       Unassigned      |       No species nodes should inherit this division assignment        | 9       |       VRL     |       Viruses |               | 10      |       VRT     |       Vertebrates     |               | 11      |       ENV     |       Environmental samples   |       Anonymous sequences cloned directly from the environment        |

代码如下:

perl /share/work/huangls/piplines/01.script/split_taxid_ncbiv2.pl division.dmp nodes.dmp prot.accession2taxid.gz nr.gz ./ #ls *fa|while read a;do b=${a%%.fa};echo mv "$a ${b}_nr.fa";done /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in INV_nr.fa -dbtype prot -title INV_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in PLN_nr.fa -dbtype prot -title PLN_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in MAM_nr.fa -dbtype prot -title MAM_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in PHG_nr.fa -dbtype prot -title PHG_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in PRI_nr.fa -dbtype prot -title PRI_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in ROD_nr.fa -dbtype prot -title ROD_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in SYN_nr.fa -dbtype prot -title SYN_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in UNA_nr.fa -dbtype prot -title UNA_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in VRL_nr.fa -dbtype prot -title VRL_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in VRT_nr.fa -dbtype prot -title VRT_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in ENV_nr.fa -dbtype prot -title ENV_nr -parse_seqids /share/work/biosoft/blast/ncbi-blast-2.6.0+/bin/makeblastdb -in BCT_nr.fa -dbtype prot -title BCT_nr -parse_seqids

perl代码, 此代码需要320G以上的内存,运行时间约3天

split_taxid_ncbiv2.pl

#The gi_taxid_nucl.dmp is about 160 MB and contains  two columns: the nucleotide's gi and taxid. #The gi_taxid_prot.dmp is about 17 MB and contains two columns:  the protein's gi  and taxid. #Divisions file (division.dmp): #       division id                             -- taxonomy database division id #       division cde                            -- GenBank division code (three characters) #       division name                           -- e.g. BCT, PLN, VRT, MAM, PRI... #       comments #0       |       BCT     |       Bacteria        |               | #1       |       INV     |       Invertebrates   |               | #2       |       MAM     |       Mammals |               | #3       |       PHG     |       Phages  |               | #4       |       PLN     |       Plants and Fungi        |               | #5       |       PRI     |       Primates        |               | #6       |       ROD     |       Rodents |               | #7       |       SYN     |       Synthetic and Chimeric  |               | #8       |       UNA     |       Unassigned      |       No species nodes should inherit this division assignment        | #9       |       VRL     |       Viruses |               | #10      |       VRT     |       Vertebrates     |               | #11      |       ENV     |       Environmental samples   |       Anonymous sequences cloned directly from the environment        | #nodes.dmp file consists of taxonomy nodes. The description for each node includes the following #fields: #        tax_id                                  -- node id in GenBank taxonomy database #        parent tax_id                           -- parent node id in GenBank taxonomy database #        rank                                    -- rank of this node (superkingdom, kingdom, ...) #        embl code                               -- locus-name prefix; not unique #        division id                             -- see division.dmp file #        inherited div flag  (1 or 0)            -- 1 if node inherits division from parent #        genetic code id                         -- see gencode.dmp file #        inherited GC  flag  (1 or 0)            -- 1 if node inherits genetic code from parent #        mitochondrial genetic code id           -- see gencode.dmp file #        inherited MGC flag  (1 or 0)            -- 1 if node inherits mitochondrial gencode from parent #        GenBank hidden flag (1 or 0)            -- 1 if name is suppressed in GenBank entry lineage #        hidden subtree root flag (1 or 0)       -- 1 if this subtree has no sequence data yet #        comments                                -- free-text comments and citations die "perl $0 <division.dmp> <nodes.dmp> <gi_taxid_n/p.dmp> <nr.fa/nt.fa> <od>" unless(@ARGV==5); use Math::BigFloat; use Bio::SeqIO; use Bio::Seq; use Data::Dumper; use PerlIO::gzip; use FileHandle; use Cwd qw(abs_path getcwd); if ($ARGV[3]=~/gz$/){         open $Fa, "<:gzip", "$ARGV[3]" or die "$!"; }else{         open $Fa, "$ARGV[3]" or die "$!"; } my$od=$ARGV[4]; $od||=getcwd; $od=abs_path($od); unless(-d $od){ mkdir $od;} my $in  = Bio::SeqIO->new(-fh => $Fa ,                                -format => 'Fasta'); my %division2name=(); if ($ARGV[0]=~/gz$/){         open IN, "<:gzip", "$ARGV[0]" or die "$!"; }else{         open IN,"$ARGV[0]" or die "$!"; } while(<IN>){         chomp;         my($taxid,$name)=split(/\s+\|\s+/);         $division2name{$taxid}=$name; } close(IN); print Dumper(\%division2name); #die; my%fout=(); for my$i (keys %division2name){         my$FO=FileHandle->new(">$od/$division2name{$i}.fa");         my $out = Bio::SeqIO->new(-fh => $FO ,  -format => 'Fasta');         $fout{$i}=$out; } print "$ARGV[0] readed\n"; #print Dumper(\%fout); #print Dumper(\%division2name); if ($ARGV[1]=~/gz$/){         open IN, "<:gzip", "$ARGV[1]" or die "$!"; }else{         open IN,"$ARGV[1]" or die "$!"; } my%taxid2division=(); while(<IN>){         chomp;         my@tmp=split(/\|/);         #for($i=0;$i<@tmp;$i++){         #       $tmp[$i]=~s/^\s+|\s+$//g;         #}         $tmp[0]=~s/^\s+|\s+$//g;         $tmp[4]=~s/^\s+|\s+$//g;         $taxid2division{$tmp[0]}=$tmp[4];         #last if $.>100; } close(IN); print "$ARGV[1] readed\n"; my %gi2taxid=(); my %prot2gi=(); if ($ARGV[2]=~/gz$/){         open IN, "<:gzip", "$ARGV[2]" or die "$!"; }else{         open IN,"$ARGV[2]" or die "$!"; } while(<IN>){         chomp;         my@tmp=split(/\s+/);         $gi2taxid{$tmp[3]}=$tmp[2];         $prot2gi{$tmp[1]}=$tmp[3]         #last if $.>100; } close(IN); print "$ARGV[2] readed\n"; #print Dumper(\%gi2taxid); $out="nr.gi.fa.gz"; open my $GZ ,"| gzip >$out" or die $!; my$out_nr = Bio::SeqIO->new(-fh => $GZ , -format => 'fasta'); while( my $seq = $in->next_seq()){         my $id=$seq->id;         my $sequence=$seq->seq;         my $desc=$seq->desc;         #my ($gi)=($id=~/gi\|(\d+)\|ref\|/);         if( exists $prot2gi{$id}){                 my $s=Bio::Seq->new(-seq=>$sequence,                     -id=>"gi|$prot2gi{$id}|ref|$id|",                     -desc=>$desc                 );                 $out_nr->write_seq($s);                 my$gi=$prot2gi{$id};                 if(exists($gi2taxid{$gi}) and exists($taxid2division{$gi2taxid{$gi}})){                         $fout{$taxid2division{$gi2taxid{$gi}}}->write_seq($s);                 }else{                         print "unknown tax for gi: $gi\n";                 }         }else{                 print "unknown prot id :$id\n";         } } $out_nr->close(); $in->close(); for my $i (keys %fout){         $fout{$i}->close(); }

以上就是“R语言 NCBI蛋白数据库分库的方法”这篇文章的所有内容，感谢各位的阅读！相信大家阅读完这篇文章都有很大的收获，小编每天都会为大家更新不同的知识，如果还想学习更多的知识，请关注亿速云行业资讯频道。