support for gene_type

[nerilbl]

I'm trying to figure out the closest ncbi-datasets equivalent of this edirect commands:

esearch -db nuccore -query "ribosomal RNA[title] AND srcdb_refseq[PROP] AND 151:7000[SLEN]" | efetch -format fasta > "rrna_genes_refseq.fasta"

So filtering genes by sequence length (151 to 7000), and ribosomal rna" in title.
In current --include with:

      --include string(,string)    Specify the data files to include (comma-separated).
                                     * gene:           gene sequence
                                     * rna:            transcript
                                     * protein:        amino acid sequences
                                     * cds:            nucleotide coding sequences
                                     * 5p-utr:         5'-UTR
                                     * 3p-utr:         3'-UTR
                                     * product-report: gene transcript and protein locations and metadata

Would it be possible to add a rrna/trna/non-coding to that list?
Thanks for maintaining and developing this tool!

[olearyna]

Hi nerilbl,

Thanks for the feedback! At this time, we do not offer a service that retrieves sequences directly from the Nucleotide or Protein databases. However, this is a feature we are actively working on and plan to support in the future. In the meantime, you can access RefSeq ribosomal genes through our Gene service. Only the API supports filtering by gene type; this filter has not yet been implemented in the CLI.
You can retrieve rRNA gene IDs for a taxon using the following API call (note that paging is required, with a maximum page size of 1000):

curl -sX GET "https://api.ncbi.nlm.nih.gov/datasets/v2/gene/taxon/1returned_content=IDS_ONLY&types=rRNA&tax_search_subtree=true" \
     -H 'accept: application/json' | jq .

Example output

{
  "reports": [
    { "gene": { "gene_id": "948466" }, "query": ["511145"] },
    …
    { "gene": { "gene_id": "947771" }, "query": ["511145"] }
  ],
  "total_count": 1258783,
  "next_page_token": "eNrj4os2szAy0VGwNDE3NzeMBQAYrgM0"
}

The Gene service returns all genes associated with GeneIDs. If you are specifically looking for the curated RefSeq marker gene set, it is available via FTP. You can find more details in the project description here: https://www.ncbi.nlm.nih.gov/refseq/targetedloci/

Please let me know if you have any questions.

Nuala

[nerilbl]

Thank @olearyna !
The direct API call and targetedloci suggestions are very helpful, I might try them soon.
Partially related, is there a way to copy all of Genbank genomes? (not just refseq)
Something like https://ftp.ncbi.nlm.nih.gov/genomes/all/ would be great, but the approaches I tried seem pretty slow.
Specifically I tried using the genbank_assembly_summary.txt to fetch specific files using the ftp_path attribute, using requests and aria2c, and I've tried using lftp:

lftp ftp://ftp.ncbi.nlm.nih.gov -e "
       --include-glob='*/md5checksums.txt' \
       /genomes/genbank/ \
       ./local_genbank/
bye
"

which I suspect was very slow because of the glob/wildcard in path.
FWIW, my attempts are sort of documented in trrna_stats.ipynb. I eventually (or at least for now) decided to fetch the BLAST databases of the rRNA sequences, and convert them to fasta using blastcmd. This is pretty messy but requires much less download requests...
So, Is there a better way to get the data?
I'm mostly interested in genbank+refseq and in getting the "genomic.fna" and gff3 files. Once locally I could use gff2parquet or similar to get the nucleic regions of the different gene_types using the fna + gff files (I also want to extract e.g. only tRNAs, which I don't think I can find a blastdb for at this moment =\