idseq基因测序流程中使用blastn_match_taxid（物种谱系）_20211213

import click
import pandas as pd
import os

@click.command()
@click.option('-r', '--reference', help='accession_file(nucl_gb.accession2taxid/prot.accession2taxid/...)', required=True)
@click.option('-i', '--inputfile', help='search_file(number_nt_virus_blastn.m8/number_nt_bacteria_blastn.m8/...)', required=True)
@click.option('-o1', '--outfile1', help='blastn_match_taxid(number_nt_virus_blastn_match_taxid/...)',required=True)#default='result.fa' ,type=float
@click.option('-o2', '--outfile2', help='blastn_match_taxid_totaxid(number_nt_virus_blastn_match_taxid_totaxid(onlytaxid)/...)',required=True)
@click.option('-o3', '--outfile3', help='blastn_match_taxid_totaxid_lineage(number_nt_virus_blastn_match_taxid_totaxid_lineage.txt/...)',required=True)
@click.option('-o4', '--outfile4', help='number_nt_virus_blastn_match_taxid_lineage/..',required=True)
@click.option('-o5', '--outfile5', help='number_nt_virus_blastn_match_taxid_lineage_add_kgs/...',required=True)
@click.option('-o6', '--outfile6', help='number_nt_virus_blastn_match_taxid_lineage_add_kgs_delrepeat/...',required=True)
@click.option('-o7', '--outfile7', help='number_nt_virus_blastn_match_taxid_lineage_add_kgs_delrepeat1/...',required=True)
@click.option('-o8', '--outfile8', help='number_nt_virus_blastn_match_taxid_lineage_add_kgs_delrepeat1_out.txt/...',required=True)

def get_blastn_match_taxid_lineage(reference,inputfile,outfile1,outfile2,outfile3,outfile4,outfile5,outfile6,outfile7,outfile8):
    accession_file1 = reference
    df_accession = pd.read_csv(accession_file1, sep='\t')
    df_accession1=df_accession.iloc[:,1:3]

    accession_to_taxid_dict=dict(zip(df_accession1['accession.version'],df_accession1['taxid']))

    search_file1=inputfile
    df_search=pd.read_csv(search_file1,sep='\t',header=None)
    df_search.columns=['Query id','accession.version','% identity','alignment length','mismatches','gap openings',
                       'q. start','q. end','s. start','s. end','e-value','bit score']

    df_search['taxid']=df_search['accession.version'].apply(lambda x : accession_to_taxid_dict[x])
    df_search.to_csv(outfile1,index=False,sep='\t',header=None)

    df_search['taxid'].to_csv(outfile2,index=False,sep='\t',header=None)

    val1=os.system("taxonkit lineage %s | taxonkit reformat -f '{k};{g};{s}' |cut -f 1,3 | tee %s |paste %s %s > %s" % (outfile2,outfile3,outfile3,outfile1,outfile4)) #https://blog.csdn.net/njafei/article/details/72764990
    if val1==0:
        print("taxid匹配上lineage成功")
    else:
        print("taxid匹配上lineage失败")

    file1 = outfile3
    df1 = pd.read_csv(file1, sep='\t', header=None)
    # print(df1)

    file2 = outfile4
    df2 = pd.read_csv(file2, sep='\t', header=None)
    # print(df2)

    df3 = df2.iloc[:, 0:df2.shape[1] - 1]
    df3.columns = ['taxid', 'taxnomy_k_g_s', 'Query id', 'accession.version', '% identity', 'alignment length',
                   'mismatches', 'gap openings',
                   'q. start', 'q. end', 's. start', 's. end', 'e-value', 'bit score']
    # print(df3)
    # print(df3.loc[2])

    df4 = df3['taxnomy_k_g_s'].str.split(';', expand=True)
    df4.columns = ['kingdom', 'genus', 'species']
    # print(df4)

    df5 = pd.concat([df3, df4], axis=1, names=['kingdom', 'genus', 'species'])
    # print(df5)
    # print(df5.columns)

    df5.to_csv(outfile5,index=False, sep='\t')

    file1 = outfile5
    df1 = pd.read_csv(file1, sep='\t')
    # print(df1)

    # print(df1.loc[:,'% identity'])
    print(df1['% identity'].value_counts(ascending=False))

    for i in range(len(df1['% identity'])):
        if df1['% identity'][i] > 80:
            i += 1
        else:
            print("unmeet")
            break

    list1 = df1['Query id'].unique().tolist()
    # print(type(df1['Query id'].unique()))
    # print(list1)
    # print(len(list1),df1.shape[0])

    # for i,item in df1.iterrows():
    # print(i,item.Query_id,item.species)
    # print(df1['Query id'].items)

    list3 = []
    for i in range(len(list1)):
        list2 = df1[df1['Query id'] == list1[i]]['species'].tolist()
        # print(list2)
        # m=m+1
        # print(m)
        #   if len(set(list2))==1               怎么判定一个List里面的元素是不是全部一样
        for j in range(1, len(list2)):
            if str(list2[j]).split(' ', 2)[0] != str(list2[0]).split(' ', 2)[0]:
                # print("不一样！")
                # break
                # print(list1[i])
                list3.append(list1[i])
                break
            else:
                # print("一样")
                j = j + 1

    # print(list3)
    for i in range(len(list3)):
        df1 = df1[~df1['Query id'].isin([str(list3[i])])]  # 删除df表中包含指定字符串的行数据
    # print(df1)

    df2 = df1.reset_index(drop=True)  # 重建索引
    # print(df2)
    df2.to_csv(outfile6,index=None, sep='\t')

    file2 =outfile6
    df3 = pd.read_csv(file2, sep='\t')
    df4 = df3.drop_duplicates(subset='Query id').reset_index()
    df4.to_csv(outfile7,index=None, sep='\t')

    #转换输出
    file1 = outfile7
    df1 = pd.read_csv(file1, sep='\t')
    # print(df1)

    df2 = df1.taxid.value_counts()
    # print(type(df2))

    # id_to_taxnomy_dict=link.set_index('taxid')['taxnomy_k_g_s'].to_dict()
    df3 = df1['taxid'].value_counts(ascending=False)
    # print(df3)

    dict_df3 = {'taxid': df3.index, 'count': df3.values}
    # print(type(dict_df3))
    df4 = pd.DataFrame(dict_df3)
    # print(df4)
    # print(df3.iloc[:,1])

    # dict1=df1.to_dict()
    # print(dict1)
    # dict2=df1.set_index('taxid').to_dict()
    # print(dict2)

    dict_id_to_taxnomy = dict(zip(df1['taxid'], df1['taxnomy_k_g_s']))
    # print(dict_id_to_taxnomy)
    df4['taxnomy_k_g_s'] = df4['taxid'].apply(lambda x: dict_id_to_taxnomy[x])
    # print(df4)

    df5 = df4['taxnomy_k_g_s'].str.split(';', expand=True)
    df5.columns = ['kingdom', 'genus', 'species']

    df6 = pd.concat([df4, df5], axis=1, names=['kingdom', 'genus', 'species'])
    # print(df6)

    del df6['taxnomy_k_g_s']
    # print(df6)

    df6['type'] = 0
    # df6['genus_Cname']='0'
    # print(df6.loc[1, 'kingdom'])
    # print(df6.loc[:, 'type'])

    for i in range(len(df6)):
        if df6.loc[i, 'kingdom'] in 'Viruses':
            df6.loc[i, 'type'] = '病毒'
        elif df6.loc[i, 'kingdom'] == 'Bacteria':
            df6.loc[i, 'type'] = '细菌'
        elif df6.loc[i, 'kingdom'] in 'fungi':
            df6.loc[i, 'type'] = '真菌'
    # print(df6)

    df6['Name'] = df6['species']  # col_names.insert(0,'Name')  #name=species name ,level 1=kingdom
    df6['level 1'] = df6['kingdom']
    col_names = df6.columns.tolist()
    # print(col_names)

    col_names_index = ['Name', 'level 1', 'taxid', 'count', 'kingdom', 'genus', 'species', 'type']
    df7 = df6.reindex(columns=col_names_index)
    # print(df7)

    df7.to_csv(outfile8,index=None, sep='\t')


if __name__ == '__main__':
    get_blastn_match_taxid_lineage()