使用idseq组装contig写说明并标注鉴定_20211206

import os
import pandas as pd
import subprocess
# subprocess.call(["tshark", "-r", "ethernet-0-0.pcap", "-Tfields", "-e", "frame.time_epoch", "-e", "ip.src", "-e", "ip.dst", ">",'linux.txt'])
# os.system('tshark -r ethernet-0-0.pcap -Tfields -e frame.time_epoch -e ip.src -e ip.dst > linux.txt')
def get(val):
    if val != 0:
        print("There are errors in steps!"+":"+str(val))
        exit()

val_inital=os.system("echo '进入容器' && docker exec -it 823fa78ed660 /bin/bash")
get(val_inital)

# gsnap建立nt库
val0=os.system("echo 'gsnap建立nt库' && gmap_build -d nt_virus_k16 -k 16 -D /idseqflow-dockerfile-container-share/20211109/GMAP-GSNAP_data/nt_virus_db /idseqflow-dockerfile-container-share/20211109/GMAP-GSNAP_data/genomeFastaFiles/Viruses_nt_gb.fa")
get(val0)

# 使用nt库进行gsnap比对，得到nt库比对结果
val1=os.system("echo '使用nt库进行gsnap比对，得到nt库比对结果' && gsnapl -A m8 --batch=0 --use-shared-memory=0 --gmap-mode=none --npaths=100 --ordered -t 36 --max-mismatches=40 -D /idseqflow-dockerfile-container-share/20211109/GMAP-GSNAP_data/nt_virus_db -d nt_virus_k16 /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_priceseq_dedup_lzw_bowtie2_subsample_gsnap.fasta > /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_priceseq_dedup_lzw_bowtie2_subsample_gsnapfilter_gasnap.m8")
get (val1)

# 进行组装
val2=os.system("echo '进行组装' && spades.py -s /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_priceseq_dedup_lzw_bowtie2_subsample_gsnap.fasta -o /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_spades_out/ -m 250 -t 32 --only-assembler")
get (val2)

# 用nt库gsnap比对的结果去建立contig的blast库-提取nt的gsnap比对结果的accession访问号
val3 = os.system("echo '用nt库gsnap比对的结果去建立contig的blast库' && cat /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_priceseq_dedup_lzw_bowtie2_subsample_gsnapfilter_gasnap.m8 | cut -f 2 |sort | uniq > /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_priceseq_dedup_lzw_bowtie2_subsample_gsnapfilter_gasnap.m8_accession")
get (val3)

# 根据accession访问号提取nt库的相关序列
val4=os.system("echo '提取nt库的相关序列' && exit && docker cp 823fa78ed660:/idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_priceseq_dedup_lzw_bowtie2_subsample_gsnapfilter_gasnap.m8_accession /home/lijing/nt_nr_database_division/ && cd /home/lijing/nt_nr_database_division/ && seqtk subseq nt /home/lijing/K200005nt_nr_database_division/K200005137_L01_126_priceseq_dedup_lzw_bowtie2_subsample_gsnapfilter_gasnap.m8_accession > /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus.refseq.fasta")
get (val4)

# 使用nt库的序列建blast库
val5= os.system("echo '使用nt库的序列建blast库' && docker exec -it 823fa78ed660 /bin/bash && mkdir -p /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/blastdb/K200005137_L01_126__nt_virus.refseq_blast_index/ && makeblastdb -in /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_nt_virus.refseq.fasta -dbtype nucl -out /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/blastdb/K200005137_L01_126__nt_virus.refseq_blast_index/K200005137_L01_126__nt_virus.refseq")
get (val5)

# 比对contig(使用blastn比对nt，blastx比对nr)
val6=os.system("echo '比对contig(使用blastn比对nt，blastx比对nr)' && blastn -query /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_spades_out/contigs.fasta -db /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/blastdb/K200005137_L01_126__nt_virus.refseq_blast_index/K200005137_L01_126__nt_virus.refseq -out /idseqflow-dockerfile-container-share/20211109/testdata/K200005137_L01_126/K200005137_L01_126_nt_virus_blastn.m8 -outfmt 6 -num_alignments 5 -num_threads 32 && exit")
get (val6)

# accession匹配taxid,K200005137_L01_126_nt_virus_blastn_match_taxid.py
accession_file1 = '/home/lijing/nt_nr_database_division/nucl_gb.accession2taxid'
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="/home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn.m8"
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('/home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid',index=False,sep='\t',header=None)

#K200005137_L01_126_nt_virus_blastn_match_taxid_lineage,匹配id，自定义谱系
val7=os.system("echo '匹配id，自定义谱系' && cat /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid | cut -f 13 > /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_totaxid && taxonkit lineage /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_totaxid | taxonkit reformat -f '{k};{g};{s}' | cut -f 1,3 | tee /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_totaxid_lineage.txt && paste /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_totaxid_lineage.txt /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid > /home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_lineage")
get (val7)

# 形成所有可用信息,K200005137_L01_126_nt_virus_blastn_match_taxid_lineage_k_g_s.py
file1="/home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_totaxid_lineage.txt"
df1=pd.read_csv(file1,sep='\t',header=None)
# print(df1)

file2="/home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_lineage"
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('/home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_lineage_add_kgs', index=False,sep='\t')

# 转换输出
file1 = "/home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_lineage_add_kgs"
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("/home/lijing/nt_nr_database_division/K200005137_L01_126_nt_virus_blastn_match_taxid_lineage_add_kgs_output.txt",index=None,sep='\t')