2016-02-20 74 views
0

我試圖從ncbi下載與一個有機體相關的所有fasta文件。從ncbi下載多個fasta文件

我試圖wget -r -l3 -A "*.fna.gz" ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/Microcystis_aeruginosa/獲得從第三級倒在.fna.gz結尾的文件,但當時它只是拒絕一切與以下的輸出:

刪除「ftp.ncbi.nlm.nih.gov/基因組/以RefSeq /菌/ Microcystis_aeruginosa/latest_assembly_versions/.listing」。 拒絕「GCF_000010625.1_ASM1062v1」。 拒絕「GCF_000307995.1_ASM30799v2」。 拒絕「GCF_000312165.1_ASM31216v1」。 拒絕「GCF_000312185.1_ASM31218v1」。 拒絕「GCF_000312205.1_ASM31220v1」。 拒絕「GCF_000312225.1_ASM31222v1」。 拒絕「GCF_000312245.1_ASM31224v1」。 拒絕「GCF_000312265.1_ASM31226v1」。 拒絕「GCF_000312285.1_ASM31228v1」。 拒絕「GCF_000312725.1_ASM31272v1」。 拒絕「GCF_000330925.1_MicAerT1.0」。 拒絕「GCF_000332585.1_MicAerD1.0」。 拒絕「GCF_000412595.1_spc777-v1」。 拒絕「GCF_000599945.1_Mic70051.0」。 拒絕「GCF_000787675.1_ASM78767v1」。 拒絕「GCF_000981785.1_ASM98178v1」。

關於爲什麼拒絕這些目錄的任何想法?謝謝你的幫助。

+0

我在想你在他們的服務器上太頻繁地請求太多,所以他們把你踢出去了。你應該真的編寫一個shell腳本,在每個wget之間休眠,以免重載服務器。 – Seekheart

回答

0

不完全確定它爲什麼會拒絕你的請求,但是當我還在做這種事情時,我發現如果我不下載小批量查詢,NCBI服務器會將我的IP定時並阻止我的IP而之前我可以再次下載。這似乎不是你看到的同樣的問題,但也許這個腳本可能會完成相同的事情。讓我知道這是否有幫助。

#!/usr/bin/env python 

from Bio import Entrez 

search_term = raw_input("Organism name: ") 

Entrez.email = "[email protected]" # required by NCBI 
search_handle = Entrez.esearch(db="nucleotide", term=search_term, usehistory="y") 
search_results = Entrez.read(search_handle) 
search_handle.close() 

gi_list = search_results["IdList"] 
count = int(search_results["Count"]) 
webenv = search_results["WebEnv"] 
query_key = search_results["QueryKey"] 

batch_size = 5 # download sequences in batches so NCBI doesn't time you out 

with open("ALL_SEQ.fasta", "w") as out_handle: 
    for start in range(0, count, batch_size): 
     end = min(count, start+batch_size) 
     print "Going to download record %i to %i" % (start+1, end) 
     fetch_handle = Entrez.efetch(db="nucleotide", rettype="fasta", retmode="text",retstart=start, retmax=batch_size, webenv=webenv, query_key=query_key) 
     data = fetch_handle.read() 
     fetch_handle.close() 
     out_handle.write(data) 

print ("\nDownload completed") 
+0

嗨,感謝您的幫助。不幸的是,該代碼似乎仍然超載了他們的服務器。但我實際上試圖使用基因組db而不是核苷酸db來分離全基因組。我認爲這需要使用elink將基因組數據庫中的ID與核苷酸數據庫中的數據庫進行關聯,這是數據實際存儲的位置。 – michberr

0

我發現一個perl腳本讓我接近完成這項任務從here。不幸的是,這個腳本只是返回基因組的ID,而不是實際的序列。

例如,我的輸出的頭:

GI | 425458296 |裁判| NZ_CAIN00000000.1 | NZ_CAIN01000000銅綠微囊藻PCC 9808,全基因組鳥槍法測序項目

GI | 425448636 |參考| NZ_CAIK00000000.1 | NZ_CAIK01000000銅綠微囊藻PCC 7941,全基因組鳥槍測序項目

任何perl用戶都知道發生了什麼事?

use strict; 
use LWP::Simple; 
my ($name, $outname, $url, $xml, $out, $count, $query_key, $webenv, $ids); 
my @genomeId; 
my $base = 'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/'; 
my $limit = 'wgs[prop]+AND+srcdb+refseq[prop])'; 
my @species = ('Microcystis aeruginosa'); 

foreach my $s (@species) { 
    undef @genomeId; 
    $query_key = $webenv = ''; 
    $s =~ s/ /+/g; 

    # ESearch 
    $url = $base . "esearch.fcgi?db=genome&term=$s"; 
    $xml = get($url); 
    $count = $1 if ($xml =~ /<Count>(\d+)<\/Count>/); 

    if ($count > 30) { 
    $url = $base . "esearch.fcgi?db=genome&term=$s&retmax=$count"; 
    $xml = get($url); 
    } 

    while ($xml =~ /<Id>(\d+?)<\/Id>/gs) { 
    push(@genomeId, $1); 
    } 

    $ids = join(',', @genomeId); 

    # ELink 
    $url = $base . "elink.fcgidbfrom=genome&db=nuccore&cmd=neighbor_history&id=$ids&term=$limit"; 
    $xml = get($url); 
    $query_key = $1 if ($xml =~ /<QueryKey>(\d+)<\/QueryKey>/); 
    $webenv = $1 if ($xml =~ /<WebEnv>(\S+)<\/WebEnv>/); 

    # EFetch 
    $url = $base . "efetch.fcgidb=nuccore&query_key=$query_key&WebEnv=$webenv&rettype=fasta&retmode=text"; 
    $out = get($url); 

    open (OUT, ">$s.fna"); 
    close OUT; 
}