Merge branch 'release/0.5.4'

Author: simonvh

CHANGELOG.md

@@ -6,11 +6,18 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/).
 
 ## [Unreleased]
 
-## [0.5.3] - 2019-03-19
+## [0.5.4] - 2019-03-19
+
+### Added
+
+- Downloading annotation from NCBI now implemented.
+- Genbank assemblies at NCBI can be searched and downloaded
 
 ### Fixed
 
-- Fixed Ensembl downloads
+- Fixed #23.
+- Fixed #26.
+- Fixed Ensembl downloads (#30)
 - Fixed FTP tests for CI
 
 ## [0.5.2] - 2018-09-11

README.md

@@ -39,6 +39,7 @@ utilities and make sure they are in your PATH:
 * `genePredToGtf`
 * `bedToGenePred`
 * `gtfToGenePred`
+* `gff3ToGenePred`
 
 You can find the binaries [here](http://hgdownload.cse.ucsc.edu/admin/exe/).
 

genomepy/__about__.py

@@ -1,3 +1,3 @@
 """Metadata"""
-__version__ = '0.5.3'
+__version__ = '0.5.4'
 __author__ = "Simon van Heeringen"

genomepy/provider.py

@@ -617,17 +617,23 @@ def _get_genomes(self):
 
         names = [
                 "assembly_summary_refseq.txt", 
+                "assembly_summary_genbank.txt", 
                 "assembly_summary_refseq_historical.txt",
                 ]
 
+        sys.stderr.write("Downloading assembly summaries from NCBI, this will take a while...\n")
+        seen = {}
         for fname in names:
             urlcleanup() 
             response = urlopen(self.assembly_url + "/" + fname)
             lines = response.read().decode('utf-8').splitlines()
             header = lines[1].strip("# ").split("\t")
             for line in lines[2:]:
                 vals = line.strip("# ").split("\t")
-                genomes.append(dict(zip(header, vals)))
+                # Don't repeat samples with the same BioSample ID
+                if vals[2] not in seen:   # BioSample ID
+                    genomes.append(dict(zip(header, vals)))
+                    seen[vals[2]] = 1
 
         return genomes
 
@@ -780,3 +786,67 @@ def _post_process_download(self, name, genome_dir, mask="soft"):
 
         # Rename tmp file to real genome file
         shutil.move(new_fa, fa)
+
+    def download_annotation(self, name, genome_dir, version=None):
+        """
+        Download annotation file to to a specific directory
+
+        Parameters
+        ----------
+        name : str
+            Genome / species name
+        
+        genome_dir : str
+            Directory to install annotation
+        """
+        sys.stderr.write("Downloading annotation...\n")
+        if not self.genomes:
+            self.genomes = self._get_genomes()
+        
+        for genome in self.genomes:
+            if genome["asm_name"] == name:
+                #ftp_path': 'ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/004/195/GCF_000004195.3_Xenopus_tropicalis_v9.1'
+                url = genome["ftp_path"]
+                url += "/" + url.split("/")[-1] + "_genomic.gff.gz"
+  
+        out_dir = os.path.join(genome_dir, name)
+        if not os.path.exists(out_dir):
+            os.mkdir(out_dir)
+        
+        # Download the file
+        try:
+            response = urlopen(url)
+            gff_file = out_dir + "/"+ name + ".annotation.gff.gz"
+            with open(gff_file, "wb") as f:
+                f.write(response.read())
+        except Exception:
+            sys.stderr.write("WARNING: Could not download annotation from NCBI, skipping.\n")
+            sys.stderr.write("URL: {}\n".format(url))
+
+            sys.stderr.write("If you think the annotation should be there, please file a bug report at:\n")
+            sys.stderr.write("https://github.com/simonvh/genomepy/issues\n")
+            return
+        
+        cmd = "gff3ToGenePred {0} /dev/stdout | wc -l"
+        out = sp.check_output(cmd.format(gff_file), shell=True)
+        if out.strip() == b"0":
+            sys.stderr.write("WARNING: annotation from NCBI contains no genes, skipping.\n")
+        else:
+            # Convert to BED file
+            bed_file = gff_file.replace("gff.gz", "bed")
+            cmd = "gff3ToGenePred -rnaNameAttr=gene {0} /dev/stdout | genePredToBed /dev/stdin {1} && gzip {1}"
+            ret = sp.check_call(cmd.format(gff_file, bed_file), shell=True) 
+
+            # Convert to GTF file
+            gtf_file = gff_file.replace("gff.gz", "gtf")
+            cmd = "gff3ToGenePred -geneNameAttr=gene {0} /dev/stdout | genePredToGtf file /dev/stdin {1} && gzip {1}"
+            ret = sp.check_call(cmd.format(gff_file, gtf_file), shell=True)
+        
+        readme = os.path.join(genome_dir, name, "README.txt")
+        with open(readme, "a") as f:
+            f.write("annotation url: {}\n".format(url))
+
+        return out_dir
+
+
+