蛇、熊猫和NCBI:我如何将熊猫数据与远程NCBI搜索结合起来？

Question

我对Snakemake还是很陌生的，我一直在尝试写一条规则时遇到了麻烦。

我一直试图将使用snakemake.remote.NCBI与访问熊猫数据和使用通配符结合起来，但我似乎无法做到这一点。

我有一个名为genomes.tsv的tsv文件，它有几个列，其中每一行都是一个物种。其中一个是" id“，并拥有该物种基因组的genbank id。另一种“物种”对每个物种都有一个独特的短字符串。在我的Snakefile中，genomes.tsv作为基因组被导入，只有id和物种列，然后将物种作为基因组指数并从基因组中删除。

我希望在我的工作流程中使用“物种”中的值作为通配符{物种}的值，我希望我的规则使用snakemake.remote.NCBI以fasta格式下载每个物种的基因组序列，然后将其输出到一个文件“{物种}_gen.fa”中。

from snakemake.remote.NCBI import RemoteProvider as NCBIRemoteProvider
import pandas as pd

configfile: "config.yaml"

NCBI = NCBIRemoteProvider(email=config["email"]) # email required by NCBI to prevent abuse

genomes = pd.read_table(config["genomes"], usecols=["species","id"]).set_index("species")

SPECIES = genomes.index.values.tolist()

rule all:
    input: expand("{species}_gen.fasta",species=SPECIES)

rule download_and_count:
    input:
        lambda wildcards: NCBI.remote(str(genomes[str(wildcards.species)]) + ".fasta", db="nuccore")
    output:
        "{species}_gen.fasta"
    shell:
        "{input} > {output}"

目前，试图运行我的代码会导致一个关键错误，但它说关键是物种的值，所以它应该能够从基因组中获得相应的genbank id。

编辑:这是错误

InputFunctionException in line 18 of /home/sjenkins/work/olflo/Snakefile:
KeyError: 'cappil'
Wildcards:
species=cappil

cappil是{物种}的一个有效值，我认为它应该可以用作一个索引。以下是前几行基因组，供参考：

species id  accession   name    assembly
cappil  8252558 GCA_004027915.1 Capromys_pilorides_(Desmarest's_hutia)  CapPil_v1_BIUU
cavape  1067048 GCA_000688575.1 Cavia_aperea_(Brazilian_guinea_pig) CavAp1.0
cavpor  175118  GCA_000151735.1 Cavia_porcellus_(domestic_guinea_pig)   Cavpor3.0

更新：

我尝试将输入行更改为：

lambda wildcards: NCBI.remote(str(genomes[genomes['species'] == wildcards.species].iloc[0]['id']) + ".fasta", db="nuccore")

但这给了我错误信息：

回溯(最近一次调用)：文件"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/init.py"，第547行，在snakemake export_cwl=export_cwl中)文件"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/workflow.py"，第421行，在execute dag.init() File "/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/dag.py"，中在init "/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/dag.py"，=self.update(self.update，progress=progress)文件"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/dag.py"，第603行，在update progress=progress)文件"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/dag.py"，第666行，在update_ progress=progress)文件“/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/dag.py”，第603行，(在update "/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/dag.py"，中)文件"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/jobs.py"，行655，在update_ missing_input = job.missing_input File "/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/jobs.py"，第398行，在missing_input中为f在self.input文件“/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/jobs.py”，中第399行(如果不是"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/io.py"，，在self.subworkflow_input中不是f)文件f.exists第208行，在存在的返回self.exists_remote文件的第119行中，在包装器v=func中(self，*args，**kwargs)文件"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/io.py"，第258行，在exists_remote返回self.remote_object.exists()文件第72行，在exists_remote likely_request_options = self._ncbi.guess_db_options_for_extension(self.file_ext，db=self.db，rettype=self.rettype中，文件"/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/remote/NCBI.py"，第110行，在file_ext添加版本中，file_ext = self._ncbi.parse_accession_str(self.local_file()) File "/home/sjenkins/miniconda3/envs/olflo/lib/python3.7/site-packages/snakemake/remote/NCBI.py"，第366行，在parse_accession_str assert file_ext中，“必须定义file_ext：{}.{}..。可能的值包括：{}“.format(加入，版本"，”.join(列表(self.valid_extensions))“AssertionError: file_ext必须定义：.可能的值包括: est、ss档、gb.xml、docset、fasta.xml、fasta、fasta_cds_na、摘要、txt、gp、medline、chr、flt、同系物、对齐分数、gbwithparts、seqid、fasta_cds_aa、gpc、uilist、uilist.xml、rsr、xml、gb、gene_table、gss、ft、gp.xml、acc、asn1、gbc。

Answer 1

好的，所以我得到AssertionError: file_ext must be defined:的原因是，如果它指定的文件名没有有效的Genbank登录号，NCBIRemoteProvider就无法识别文件扩展名。我给了它带有genbank it的文件名，所以它返回了这个错误。

而且，似乎整个基因组序列没有返回所有序列的登录号。取而代之的是wgs报告的登录号，然后是每个脚手架的登录号。我决定尝试下载我需要的基因组手动，而不是尝试下载所有的脚手架，然后组合它们。

Answer 2

我觉得你应该：

genomes = pd.read_table(config["genomes"], usecols=["species","id"])

SPECIES = list(genomes['species'])

然后访问给定物种的ID：

lambda wildcards: str(genomes[genomes['species'] == wildcards.species].iloc[0]['id'])