亚洲激情专区-91九色丨porny丨老师-久久久久久久女国产乱让韩-国产精品午夜小视频观看

溫馨提示×

溫馨提示×

您好,登錄后才能下訂單哦!

密碼登錄×
登錄注冊×
其他方式登錄
點擊 登錄注冊 即表示同意《億速云用戶服務條款》

R語言 NCBI蛋白數據庫分庫的方法

發布時間:2022-03-21 10:05:10 來源:億速云 閱讀:257 作者:iii 欄目:開發技術

今天小編給大家分享一下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,分類信息對所有的蛋白序列nr庫進行分類

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蛋白數據庫分庫的方法”這篇文章的所有內容,感謝各位的閱讀!相信大家閱讀完這篇文章都有很大的收獲,小編每天都會為大家更新不同的知識,如果還想學習更多的知識,請關注億速云行業資訊頻道。

向AI問一下細節

免責聲明:本站發布的內容(圖片、視頻和文字)以原創、轉載和分享為主,文章觀點不代表本網站立場,如果涉及侵權請聯系站長郵箱:is@yisu.com進行舉報,并提供相關證據,一經查實,將立刻刪除涉嫌侵權內容。

AI

虹口区| 昔阳县| 西畴县| 沙雅县| 龙州县| 湖州市| 秀山| 惠来县| 讷河市| 霍邱县| 海阳市| 嘉祥县| 西乌| 宁远县| 南京市| 定日县| 巴里| 仁怀市| 玉田县| 石门县| 哈巴河县| 盐山县| 金山区| 磐石市| 平阴县| 宿松县| 马公市| 正镶白旗| 仁寿县| 米林县| 潼南县| 出国| 子洲县| 镇安县| 都江堰市| 临夏市| 崇文区| 富锦市| 靖州| 云浮市| 阿克苏市|