Wetawrap kraken2运行报错KeyError“2304686”解决方案

本文的bug来自输入经过初步处理的宏基因组测序序列，在进行kraken2模块时，使用taxtonomy分类数据库时，如下图出现KeyError“2304686”

我们根据提示路径打开kraken2_translate.py文件:

#!/usr/bin/env python
import sys
import os


def get_full_name(taxid, names_map, ranks_map):
    """
    Generate full taxonomic lineage from a taxid
    Args:
        taxid (str): taxid
        names_map (dict[str:str]): the taxonomic names for each taxid node
        ranks_map (dict[str:str]): the parent node for each taxid
    Returns:
        taxonomy (str): full taxonomic lineage string

    """
    taxid_lineage = list()
    while True:
        taxid_lineage.append(taxid)
        if taxid not in ranks_map:
            break
        new_taxid = ranks_map[taxid]
        if taxid == new_taxid:
            break
        else:
            taxid = new_taxid

    names_lineage = list()
    for taxid in taxid_lineage:
          name = names_map[taxid]
          names_lineage.append(name)

    taxonomy = ";".join(reversed(names_lineage))
    return taxonomy


def load_kraken_db_metadata(kraken2_db):
    """
    Load NCBI taxonomic name mappings to be able to convert taxids into taxonomic strings
    Args:
        kraken2_db (str): path to kraken2 standard database location
    Returns:
        names_map (dict[str:str]): the taxonomic names for each taxid node
        ranks_map (dict[str:str]): the parent node for each taxid

    """
    print("Loading NCBI node names")
    names_path = os.path.join(kraken2_db, "taxonomy", "names.dmp")
    names_map = dict()
    with open(names_path) as input:
        for line in input:
            cut = line.rstrip().split("\t")
            taxid = cut[0]
            name = cut[2]
            type = cut[6]
            if type == "scientific name":
                names_map[taxid] = name

    print("Loading NCBI taxonomic ranks")
    ranks_path = os.path.join(kraken2_db, "taxonomy", "nodes.dmp")
    ranks_map = dict()
    with open(ranks_path) as input:
        for line in input.readlines():
            cut = line.rstrip().split("\t")
            taxid = cut[0]
            parent_taxid = cut[2]
            rank = cut[4]
            ranks_map[taxid] = parent_taxid
    return names_map, ranks_map


def translate_kraken2_annotations(annotation_file=None, kraken2_db=None, output=None):
    """
    Translate the kraken2 annotations from taxids to taxonomy strings
    Args:
        annotation_file (str): kraken2 output file
        kraken2_db (str): path to kraken2 standard database
        output (str): path to final file with translated annotations
    Returns: None

    """
    print("Translating kraken2 annotations")
    if os.path.isfile(output):
        print("Looks like the translated kraken2 output already exists. Skipping...")
        return None

    names_map, ranks_map = load_kraken_db_metadata(kraken2_db)
    print("Writing translated taxonomy names to %s" % output)
    with open(output, "w") as output:
        with open(annotation_file) as input:
            for line in input:
                cut = line.rstrip().split("\t")
                contig = cut[1]
                if cut[0] == "U":
                    taxonomy = ""
                else:
                    taxid = cut[2].split()[-1][:-1]
                    taxonomy = get_full_name(taxid, names_map, ranks_map)
                output.write("%s\t%s\n" % (contig, taxonomy))


def main():
    """
    Load in arguments and start analysis
    Returns: None

    """

    database_location = sys.argv[1]
    kraken_file = sys.argv[2]
    output_file = sys.argv[3]
    print("Translating kraken2 annotations from %s, using metadata from the kraken2 database in %s; saving to %s" \
          % (kraken_file, database_location, output_file))
    translate_kraken2_annotations(annotation_file=kraken_file, kraken2_db=database_location, output=output_file)

if __name__ == '__main__':
    """ Launch script
    """
    main()

首先进行代码的阅读和分析，def为自定义的函数，def和（）之间为定义函数名。‘’‘    ‘’‘三引号之间为注释内容，代码不运行。自上而下

我们可以看到有4个def，其中第一个和第二个def为自定义的函数，第三个def调用了第一个和第二个def函数，第四个def函数调用了第三个def函数，同时第4个def是main（）函数，也叫做主函数，代码运行的基本框架在main函数中定义，代码最后一段的if语句表明了函数的执行从main函数开始。搞清代码的执行顺序后，我们再看开始的bug提示，30、98、114、120行代码运行错误，而根据具体每一行代码执行的顺序，我可以发现第一次出现错误在第30行，第30行代码为name = names_map[taxid]，问题就出现在这个赋值语句，将右边的值赋值到左边，names_map是一个字典数据类型，里面的taxid为字典的属性，也就是说将出现在names_map中某一个具体taxid属性对应的属性值赋值给name，如果names_map中没有这个taxid，那么这个字典的赋值运算就会失败，从而给出报错，并停止运行。那么为什么这个taxid会不存在呢？这就是因为两个数据库（kraken2的hash.k2d数据库和taxonomy数据库）存在一定信息偏差，也就是数学上的非充分必要包含关系。所以代码运行错误。知道了代码运行错误的背后逻辑，我们就可以删除或者不运行没有出现在name_map中的taxid，让其余正常的taxid继续运行，直至代码完全执行完毕。

第一步删除已经生成的final_assembly.kraken2文件，将下列代码中的*****替换成上面的KeyError中的数值。然后在python中运行：

f = open("./KRAKEN/final_assembly.krak2","a")
m = open("./KRAKEN/final_assembly.krak2_pure","a")
for each in f.readlines():
  b = each.find("taxid *****")
  if b == -1:
    m.write(each)

f.close()
m.close()

执行上述操作后，会生成一个新文件final_assembly.krak2_pure，此文件会跳过（taxid *****）的一行内容，删除旧的final_assembly.krak2文件，将final_assembly.krak2_pure文件名改为final_assembly.krak2。接着在linux下执行：

metawrap kraken2 -o KRAKEN -t 10 -s 100000 ASSEMBLY/final_assembly.fasta

 如此循环往复，便可将KeyError的数值所表示的contig 一一忽略。当最后一个KeyError的数值跳过后，translate的结果也运行完毕，将得到你想要的物种分类图。

     本推文旨在引导该方法出现bug解读和解决方案，以及删除关键字所在行的代码执行方案，对于解决本文，最好的办法还是从该软件的源代码进行矫正。

      只需要在kraken2_translate.py原来的代码基础上写一个if语句，也就4行内容，不但可以做到只运行一次就可得到最终的分类图，还可以输出KeyError的taxid到一个新的文件中，新文件可用于后续的结果检查。

       如果想要修复bug后的kraken2_translate.py，请关注下方微信公众号，并留言！我会在第一时间处理！本公众号还会继续对metawrap其他模块的bug问题进行一一解决，关注我，让新手metawrap一路畅行，谢谢！

        如果喜欢请分享，谢谢！