二代测序，测试验证bwtagger和blast软件比对结果

最新推荐文章于 2024-06-10 22:18:41 发布
最新推荐文章于 2024-06-10 22:18:41 发布
阅读量383 收藏 6
点赞数 7
文章标签： python windows 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/m0_50542529/article/details/134992015
版权

1.比较bowtie2和bmtagger的种属差异

import pandas as pd
import random
inputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\NL230601-1B_S45_L002.2.classify.last.xls'

inputpath2=r'C:\Users\Administrator\Desktop\自己测试捕获\20230619\第三张图-new20230629\宏基因组\22-12373F.2.classify.last.xls'

df1=pd.read_csv(inputpath1,sep='\t')
df2=pd.read_csv(inputpath2,sep='\t')
print(df1)

##genus
df1_genus=df1[(df1['count'] >=2) & (df1['G/S']=='species')].copy().reset_index(drop=True)
list1_genus=df1_genus['taxid'].tolist()
dict1=dict(zip(df1_genus['taxid'],df1_genus['Ename']))
# print(list1_genus)

df2_genus=df2[(df2['count'] >=2) & (df2['G/S']=='species')].copy().reset_index(drop=True)
list2_genus=df2_genus['taxid'].tolist()
dict2=dict(zip(df2_genus['taxid'],df2_genus['Ename']))
# print(list2_genus)

list_1_to_2=list(set(list2_genus).difference(set(list1_genus)))
print(list_1_to_2)

# list_2_to_1=list(set(list1_genus).difference(set(list2_genus)))
# print(list_2_to_1)

list_intersection=list(set(list2_genus).intersection(set(list1_genus)))
print(list_intersection)

# list_intersection_10=random.sample(list_intersection,10)
# print(list_intersection_10)
list_intersection_10=[74426,28116,2528203,48074,2025876,208479,821,29442,2492837,357276]

# list_1_to_2_10=random.sample(list_1_to_2,10)
# print(list_1_to_2_10)
list_1_to_2_10=[1784,1423,2023057,83656,2005458,1357916,247,1969841,562,2584122]

inputpath3=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\bmtagger.kraken2'

inputpath4=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\bowtie2.kraken2'


df3=pd.read_csv(inputpath3,sep='\t',header=None)

df_add=pd.DataFrame(columns=[0])
for i in list_intersection_10:
    df3_part1=df3[df3[2]==i].copy().reset_index(drop=True)
    df3_part1['Ename']=df3_part1[2].apply(lambda x:dict1[x])
    # print(df3_part1)
    df_add=df_add.append(df3_part1)
print(df_add)
df_add.to_csv(r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\intersection_10.kraken2.xls',sep='\t',index=None,header=None)

df4=pd.read_csv(inputpath4,sep='\t',header=None)
df_add=pd.DataFrame(columns=[0])
for i in list_1_to_2_10:
    df4_part1=df4[df4[2]==i].copy().reset_index(drop=True)
    df4_part1['Ename']=df4_part1[2].apply(lambda x:dict2[x])
    # print(df3_part1)
    df_add=df_add.append(df4_part1)
print(df_add)
df_add.to_csv(r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.kraken2.xls',sep='\t',index=None,header=None)

2.extract_seq

## [83561, 74426, 28116, 214856, 2528203, 48074, 1307, 2025876, 208479, 821]
## [1969841, 1784, 277, 1423, 585, 63363, 646637, 53462, 732, 1580]

import pyfastx
import pandas as pd
inputpath1=r'/data/lijing/Cust_NAD_WGS/NL230601-1B_S45_L002_R1_001.fastq.gz'
inputpath2=r'/data/lijing/Cust_NAD_WGS/NL230601-1B_S45_L002_R2_001.fastq.gz'

inputpath3=r'/data/lijing/20230629_test_Cust_NAD_WGS/intersection_10.kraken2.xls'
inputpath4=r'/data/lijing/20230629_test_Cust_NAD_WGS/1_to_2_10.kraken2.xls'
# inputpath3=r'/data/lijing/20230629_test_Cust_NAD_WGS/NL230601-1B_S45_L002.out'
#
# outputpath1=r'/data/lijing/20230629_test_Cust_NAD_WGS/1_NL230601-1B_S45_L002.human.fq.1.gz'
# outputpath2=r'/data/lijing/20230629_test_Cust_NAD_WGS/1_NL230601-1B_S45_L002.human.fq.2.gz'

## seq
fq1 = pyfastx.Fastx(inputpath1, comment=True)
dict1 = {}
for name, seq, qual, comment in fq1:
    # print(name)
    # print(seq)
    # print(qual)
    # print(comment)
    dict1[name] = [seq, qual, comment]

fq2 = pyfastx.Fastx(inputpath2, comment=True)
dict2 = {}
for name, seq, qual, comment in fq2:
    dict2[name] = [seq, qual, comment]

df3=pd.read_csv(inputpath3,sep='\t',header=None)
with open(r'/data/lijing/20230629_test_Cust_NAD_WGS/intersection_10.kraken2.fa','w') as f1:
    for i in range(len(df3)):
        f1.write('>'+df3.loc[i,1]+'_'+str(int(df3.loc[i,2]))+':'+str(df3.loc[i,5].replace(' ','_'))+':r1'+'\n'+dict1[df3.loc[i,1]][0]+'\n'+\
                 '>'+df3.loc[i,1]+'_'+str(int(df3.loc[i,2]))+':'+str(df3.loc[i,5].replace(' ','_'))+':r2'+'\n'+dict2[df3.loc[i,1]][0]+'\n')

df4=pd.read_csv(inputpath4,sep='\t',header=None)
with open(r'/data/lijing/20230629_test_Cust_NAD_WGS/1_to_2_10.kraken2.fa','w') as f1:
    for i in range(len(df4)):
        f1.write('>'+df3.loc[i,1]+'_'+str(int(df4.loc[i,2]))+':'+str(df4.loc[i,5].replace(' ','_'))+':r1'+'\n'+dict1[df4.loc[i,1]][0]+'\n'+\
                 '>'+df3.loc[i,1]+'_'+str(int(df4.loc[i,2]))+':'+str(df4.loc[i,5].replace(' ','_'))+':r2'+'\n'+dict2[df4.loc[i,1]][0]+'\n')

3.针对在线blast比对-拆分json比对文件

import json
# dict_part={}
with open(r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.json','r') as f1,open(r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.xls','w') as f2:
    f2.write('query_title\tquery_len\taccession\ttitle\ttaxid\tsciname\tidentity\tquery_from\tquery_to\thit_from\thit_to\talign_len\tgaps\tqseq\thseq\tmidline\n')
    f1.readline()
    f1.readline()
    line_part=''
    for line in f1:
        line=line.strip('\n')
        if line.startswith(','):
            # line_part +=']}'
            # print(line_part)

            dict_part=json.loads(line_part)
            # print(dict_part['report']['results']['search'])
            # print(len(dict_part['report']['results']['search']['hits']))
            # print(dict_part['report']['results']['search']['hits'][0])
            #
            # print(dict_part['report']['results']['search']['hits'][0]['description'][0])
            # print(dict_part['report']['results']['search']['hits'][0]['description'][0]['accession'])
            # print(dict_part['report']['results']['search']['hits'][0]['description'][0]['title'])
            #
            # print(dict_part['report']['results']['search']['hits'][0]["hsps"][0]['title'])
            # print(dict_part['report']['results']['search']['hits'][1])
            for i in range(len(dict_part['report']['results']['search']['hits'])):
                for j in range(len(dict_part['report']['results']['search']['hits'][i]["hsps"])):

                    f2.write(dict_part['report']['results']['search']['query_title']+'\t'+
                             str(dict_part['report']['results']['search']['query_len'])+'\t'+
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["accession"]+'\t'+
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["title"]+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]['description'][0]["taxid"]) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["sciname"] + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["identity"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_from"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_to"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_from"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_to"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["align_len"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["gaps"]) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["qseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["midline"]+'\n')

            # break
            # f2.write(dict_part['report'])

            line_part=''

        elif line.startswith(']'):
            line_part=line_part[:-1]+'}'
            # print(line_part)
            dict_part=json.loads(line_part)
            print(type(dict_part))
            for i in range(len(dict_part['report']['results']['search']['hits'])):
                for j in range(len(dict_part['report']['results']['search']['hits'][i]["hsps"])):

                    f2.write(dict_part['report']['results']['search']['query_title']+'\t'+
                             str(dict_part['report']['results']['search']['query_len'])+'\t'+
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["accession"]+'\t'+
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["title"]+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]['description'][0]["taxid"]) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["sciname"] + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["identity"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_from"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_to"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_from"])+'\t'+
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_to"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["align_len"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["gaps"]) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["qseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["midline"]+'\n')

        else:
            line_part += line.strip()
            # line_part +=line

# print(dict_part)

4.针对在线blast比对-判断bowtie2中有但bmtagger没有的10物种中假阳性（可能为人源序列或其他）的比例

import pandas as pd
inputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.xls'
df1=pd.read_csv(inputpath1,sep='\t')


df1['query_title_id']=df1['query_title'].apply(lambda x:x.split('_')[0])
df1['query_title_taxid']=df1['query_title'].apply(lambda x:int(x.split('_')[1].split(':')[0]))
print(df1)

list1=list(set(df1['query_title'].tolist()))
print(list1)

for i in list1:
    df1_part=df1[df1['query_title']==i].copy().iloc[0:2,:].reset_index(drop=True)
    # print(df1_part)
    df1_part=df1_part[df1_part['query_title_taxid']==df1_part['taxid']]
    # df1_part=df1_part[df1_part['qseq']==df1_part['hseq']].copy()
    # print(df1_part)
    if len(df1_part)==2:
        print(df1_part)

5.blast_alignment-extract_json_result

import json

# dict_part={}
with open(r'intersection_10.kraken2.blastn.json', 'r') as f1, open(
        r'intersection_10.kraken2.blastn.xls', 'w') as f2:
    f2.write(
        'query_title\tquery_len\taccession\ttitle\tidentity\tquery_from\tquery_to\thit_from\thit_to\talign_len\tgaps\tqseq\thseq\tmidline\n')
    f1.readline()
    f1.readline()
    line_part = ''
    for line in f1:
        line = line.strip('\n')
        if line.startswith(','):
            # line_part +=']}'
            # print(line_part)

            dict_part = json.loads(line_part)
            # print(dict_part['report']['results']['search'])
            # print(len(dict_part['report']['results']['search']['hits']))
            # print(dict_part['report']['results']['search']['hits'][0])
            #
            # print(dict_part['report']['results']['search']['hits'][0]['description'][0])
            # print(dict_part['report']['results']['search']['hits'][0]['description'][0]['accession'])
            # print(dict_part['report']['results']['search']['hits'][0]['description'][0]['title'])
            #
            # print(dict_part['report']['results']['search']['hits'][0]["hsps"][0]['title'])
            # print(dict_part['report']['results']['search']['hits'][1])
            for i in range(len(dict_part['report']['results']['search']['hits'])):
                for j in range(len(dict_part['report']['results']['search']['hits'][i]["hsps"])):
                    f2.write(dict_part['report']['results']['search']['query_title'] + '\t' +
                             str(dict_part['report']['results']['search']['query_len']) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["accession"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["title"] + '\t' +
                             # str(dict_part['report']['results']['search']['hits'][i]['description'][0][
                             #         "taxid"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["identity"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_from"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_to"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_from"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_to"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["align_len"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["gaps"]) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["qseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["midline"] + '\n')

            # break
            # f2.write(dict_part['report'])

            line_part = ''

        elif line.startswith(']'):
            line_part = line_part[:-1] + '}'
            # print(line_part)
            dict_part = json.loads(line_part)
            print(type(dict_part))
            for i in range(len(dict_part['report']['results']['search']['hits'])):
                for j in range(len(dict_part['report']['results']['search']['hits'][i]["hsps"])):
                    f2.write(dict_part['report']['results']['search']['query_title'] + '\t' +
                             str(dict_part['report']['results']['search']['query_len']) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["accession"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]['description'][0]["title"] + '\t' +
                             # str(dict_part['report']['results']['search']['hits'][i]['description'][0][
                             #         "taxid"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["identity"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_from"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["query_to"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_from"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hit_to"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["align_len"]) + '\t' +
                             str(dict_part['report']['results']['search']['hits'][i]["hsps"][j]["gaps"]) + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["qseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["hseq"] + '\t' +
                             dict_part['report']['results']['search']['hits'][i]["hsps"][j]["midline"] + '\n')

        else:
            line_part += line.strip()

6.determine_accurate

import pandas as pd
inputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.kraken2.blastn.xls'
inputpath2=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.kraken2.fa'
outputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.kraken2.blastn.last.top1.xls'
## outputpath2=r'C:\Users\Administrator\Desktop\自己测试捕获\20230704_比较bowtie2和bwtagger\1_to_2_10.kraken2.blastn.inter2.xls'

# inputpath1=r'/data/lijing/20230629_test_Cust_NAD_WGS/intersection_10.kraken2.blastn.xls'
# inputpath2=r'/data/lijing/20230629_test_Cust_NAD_WGS/intersection_10.kraken2.fa'
# outputpath1=r'/data/lijing/20230629_test_Cust_NAD_WGS/intersection_10.kraken2.blastn.last.xls'
# ## outputpath2=r'/data/lijing/20230629_test_Cust_NAD_WGS/intersection_10.kraken2.blastn.inter2.xls'

list_all=[]
dict_all={}
with open(inputpath2,'r') as f3:
    for line in f3:
        line=line.strip('\n')
        # list_line=line.split(':')
        if line.startswith('>'):
            list_all.append(line.split('_',1)[1].split(':')[1].replace('_',' '))
            dict_all[line.split('_',1)[1].split(':')[1].replace('_',' ')]=line.split('_',1)[1].split(':')[0]
print(dict_all)

value_counts_all=pd.value_counts(list_all)
dict_value_count_all=dict(zip(value_counts_all.index,value_counts_all.values))
print(dict_value_count_all)
df_value_count_all=pd.DataFrame.from_dict(dict_value_count_all,orient='index',columns=['raw_reads'])
df_value_count_all=df_value_count_all.rename_axis('Ename').reset_index()

with open(inputpath1,'r') as f1,open(outputpath1,'w') as f2:
    f1.readline()
    f2.write('query_title\tdescription\tequal\n')
    for line in f1:
        line=line.strip('\n')
        list_line=line.split('\t')
        if list_line[0].split(':')[-2].replace('_',' ') in list_line[3]:
            f2.write(list_line[0]+'\t'+list_line[3]+'\t'+'yes'+'\n')
        elif list_line[0].split(':')[-2].replace('_',' ')=='[Clostridium] bolteae' and 'Enterocloster bolteae' in list_line[3]:
            f2.write(list_line[0] + '\t' + list_line[3] + '\t' + 'yes' + '\n')
        else:
            f2.write(list_line[0] + '\t' + list_line[3] + '\t' + 'no' + '\n')

## 第3次写

df2=pd.read_csv(outputpath1,sep='\t')
df2['Ename']=df2['query_title'].apply(lambda x:x.split('_', 1)[1].split(':')[1].replace('_',' '))
df2=df2[df2['equal']=='yes'].copy().reset_index(drop=True)
print(df2)
# df2=df2.groupby(["query_title"]).head(1).reset_index(drop=True)
df2=df2.groupby(["query_title"]).head(2).reset_index(drop=True)
print(df2)
print("**************"*5)

value_counts_have=pd.value_counts(df2['query_title'])
dict_value_count_have=dict(zip(value_counts_have.index,value_counts_have.values))
print(dict_value_count_have)
df_value_count_have=pd.DataFrame.from_dict(dict_value_count_have,orient='index',columns=['counts'])
df_value_count_have=df_value_count_have.rename_axis('query_title').reset_index()

# df_value_count_have=df_value_count_have[df_value_count_have['counts']==1]
df_value_count_have=df_value_count_have[df_value_count_have['counts']==2]
print(df_value_count_have)


list_have1=list(set(df_value_count_have['query_title'].tolist()))
df_value_count_have['Ename']=df_value_count_have['query_title'].apply(lambda x:x.split('_', 1)[1].split(':')[1].replace('_',' '))
# df2=df2[df2['query_title'].isin(list_have1)]
value_counts_have=pd.value_counts(df_value_count_have['Ename'])
dict_value_count_have=dict(zip(value_counts_have.index,value_counts_have.values))
print(dict_value_count_have)
df_value_count_have=pd.DataFrame.from_dict(dict_value_count_have,orient='index',columns=['real_reads'])
df_value_count_have=df_value_count_have.rename_axis('Ename').reset_index()
print(df_value_count_have)

df_merge=pd.merge(df_value_count_all,df_value_count_have,how='outer',on=['Ename'])

df_merge.fillna(0,inplace=True)
print(df_merge)
df_merge['taxid']=df_merge['Ename'].apply(lambda x:dict_all[x])
df_merge=df_merge[['Ename','taxid','raw_reads','real_reads']]
df_merge.to_csv(outputpath1,sep='\t',index=None,header=True)



# ## 第二次写
# list_have=[]
# with open(outputpath1,'r') as f3,open(outputpath2,'w') as f4:
#     f3.readline()
#     f4.write('query_title\tdescription\tequal\n')
#     for line in f3:
#         line = line.strip('\n')
#         list_line = line.split('\t')
#
#         if list_have.count(list_line[0]) >= 2:
#             pass
#         else:
#             f4.write(line+'\n')
#         list_have.append(list_line[0])
#
# df2=pd.read_csv(outputpath2,sep='\t')
# df2['Ename']=df2['query_title'].apply(lambda x:x.split('_', 1)[1].split(':')[1].replace('_',' '))
# df2=df2[df2['equal']=='yes'].copy().reset_index(drop=True)
# print(df2)
#
# value_counts_have=pd.value_counts(df2['query_title'])
# dict_value_count_have=dict(zip(value_counts_have.index,value_counts_have.values))
# print(dict_value_count_have)
# df_value_count_have=pd.DataFrame.from_dict(dict_value_count_have,orient='index',columns=['counts'])
# df_value_count_have=df_value_count_have.rename_axis('query_title').reset_index()
#
# df_value_count_have=df_value_count_have[df_value_count_have['counts']==2]
# print(df_value_count_have)
#
#
# list_have1=list(set(df_value_count_have['query_title'].tolist()))
# df_value_count_have['Ename']=df_value_count_have['query_title'].apply(lambda x:x.split('_', 1)[1].split(':')[1].replace('_',' '))
# # df2=df2[df2['query_title'].isin(list_have1)]
# value_counts_have=pd.value_counts(df_value_count_have['Ename'])
# dict_value_count_have=dict(zip(value_counts_have.index,value_counts_have.values))
# print(dict_value_count_have)
# df_value_count_have=pd.DataFrame.from_dict(dict_value_count_have,orient='index',columns=['real_reads'])
# df_value_count_have=df_value_count_have.rename_axis('Ename').reset_index()
# print(df_value_count_have)
#
# df_merge=pd.merge(df_value_count_all,df_value_count_have,how='outer',on=['Ename'])
#
# df_merge.fillna(0,inplace=True)
# print(df_merge)
# df_merge['taxid']=df_merge['Ename'].apply(lambda x:dict_all[x])
# df_merge=df_merge[['Ename','taxid','raw_reads','real_reads']]
# df_merge.to_csv(outputpath1,sep='\t',index=None,header=True)

## 第一次写
# df1=pd.read_csv(outputpath1,sep='\t')
#
# list1=list(set(df1['query_title'].tolist()))
# print(list1)
#
# list_have=[]
# for i in list1:
#     df1_part=df1[df1['query_title']==i].copy().iloc[0:2,:]
#     df1_part=df1_part[df1_part['equal']=='yes'].reset_index(drop=True)
#
#     if len(df1_part)==2:
#         # print(df1_part)
#         # list_have.append(df1_part.loc[0,'query_title'].split('_',1)[1].rsplit(':',1)[0])
#         list_have.append(df1_part.loc[0, 'query_title'].split('_', 1)[1].split(':')[1].replace('_',' '))
# print(list_have)
#
# value_counts_have=pd.value_counts(list_have)
# dict_value_count_have=dict(zip(value_counts_have.index,value_counts_have.values))
# print(dict_value_count_have)
# df_value_count_have=pd.DataFrame.from_dict(dict_value_count_have,orient='index',columns=['real_reads'])
# df_value_count_have=df_value_count_have.rename_axis('Ename').reset_index()
# print(df_value_count_have)
#
# df_merge=pd.merge(df_value_count_all,df_value_count_have,how='outer',on=['Ename'])
#
# df_merge.fillna(0,inplace=True)
# print(df_merge)
# df_merge.to_csv(outputpath1,sep='\t',index=None,header=True)

7.extract_seq_base_human_id

import pyfastx
import gzip
import argparse

parser=argparse.ArgumentParser(description='input error !')
parser.add_argument('-ip1','--inputpath1',required=True,help='input like NL230601-1B_S45_L002.r1.paired.fq.gz')
parser.add_argument('-ip2','--inputpath2',required=True,help='input like NL230601-1B_S45_L002.r2.paired.fq.gz')
parser.add_argument('-ip3','--inputpath3',required=True,help='input like NL230601-1B_S45_L002.out')
parser.add_argument('-op1','--outputpath1',required=True,help='NL230601-1B_S45_L002.human.fq.1.gz')
parser.add_argument('-op2','--outputpath2',required=True,help='NL230601-1B_S45_L002.human.fq.2.gz')
args = parser.parse_args()

inputpath1=args.inputpath1
inputpath2=args.inputpath2
inputpath3=args.inputpath3
outputpath1=args.outputpath1
outputpath2=args.outputpath2


# inputpath1=r'/data/lijing/20230629_test_Cust_NAD_WGS/bmtagger_minikraken2_test/NL230601-1B_S45_L002.r1.paired.fq.gz'
# inputpath2=r'/data/lijing/20230629_test_Cust_NAD_WGS/bmtagger_minikraken2_test/NL230601-1B_S45_L002.r2.paired.fq.gz'
#
# inputpath3=r'/data/lijing/20230629_test_Cust_NAD_WGS/NL230601-1B_S45_L002.out'
#
# outputpath1=r'/data/lijing/20230629_test_Cust_NAD_WGS/1_NL230601-1B_S45_L002.human.fq.1.gz'
# outputpath2=r'/data/lijing/20230629_test_Cust_NAD_WGS/1_NL230601-1B_S45_L002.human.fq.2.gz'

## seq
fq1 = pyfastx.Fastx(inputpath1, comment=True)
dict1 = {}
for name, seq, qual, comment in fq1:
    # print(name)
    # print(seq)
    # print(qual)
    # print(comment)
    dict1[name] = [seq, qual, comment]

fq2 = pyfastx.Fastx(inputpath2, comment=True)
dict2 = {}
for name, seq, qual, comment in fq2:
    dict2[name] = [seq, qual, comment]

## human id
list_human_id=[]
with open(inputpath3,'r') as f1:
    for line in f1:
        line=line.strip('\n')
        list_human_id.append(line)
f1.close()
# print(list_human_id)
# dict_remove_human1={ item:dict1[item] for item in dict1 if item not in list_human_id}
# dict_remove_human2={ item:dict2[item] for item in dict2 if item not in list_human_id}
for i in list_human_id:
    del dict1[i]
dict_remove_human1=dict1
# print(dict_remove_human1)
for i in list_human_id:
    del dict2[i]
dict_remove_human2=dict2


with gzip.open(outputpath1,'wb') as f1:
    for i in dict_remove_human1:
        line='@'+i+' '+dict_remove_human1[i][2]+'\n'+dict_remove_human1[i][0]+'\n+\n'+dict_remove_human1[i][1]+'\n'
        line=line.encode(encoding='utf-8')
        f1.write(line)

with gzip.open(outputpath2,'wb') as f2:
    for i in dict_remove_human2:
        line='@'+i+' '+dict_remove_human2[i][2]+'\n'+dict_remove_human2[i][0]+'\n+\n'+dict_remove_human2[i][1]+'\n'
        line=line.encode(encoding='utf-8')
        f2.write(line)
lj_learn
关注
  • 7
    点赞
  • 踩
  • 6
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
kneaddata安装及使用
dys的博客
09-28 6233
文章目录简介安装下载数据库创建定制数据库(可选)运行参考资料 简介 安装 conda create -n kneaddata conda activate kneaddata conda install -c biobakery kneaddata 下载数据库 mkdir kneaddata_database cd kneaddata_database/ kneaddata_database --download human_genome bowtie2 ./ kneaddata_database --d
基因二代测序-分别对blast和bwa比对结果进行统计比较-20230506
12-15
基因二代测序-分别对blast和bwa比对结果进行统计比较-20230506
病毒组学数据分析 -01 kneaddata数据质控
Kun_98的博客
01-12 1934
kneaddata软件介绍，使用bowtie2去除人体DNA污染，并通过trimmomatic对序列进行筛选，最后fastqc报告序列质控结果
Microbiome：宏基因组分箱流程MetaWRAP分析实战和结果解读
刘永鑫的博客——宏基因组公众号
10-17 2万+
文章目录MetaWRAP—a flexible pipeline for genome-resolved metagenomic data analysis分析实战0.下载肠道宏基因组数据1.read_qc质控和去宿主2. 组装3. kraken物种注释(可选)4. 三种方法Bin5.Bin提纯6. Blobology可视化bin7. Bin定量8.重组装9.bin物种注释10.Bin功能注释我们...
二代测序手册
09-11
对生信分析相当有用，值得推荐，如果有R语言基础则更好，希望对从事相关研究的人有不小的帮助，本来打算免费下载，可考虑到自己也没积分，就当是创收吧。呵呵
基因组二代测序数据的自动化分析流程
05-10
基因组二代测序数据的自动化分析流程基因组二代测序数据的自动化分析流程
父母基因型验证对靶向二代测序结果解读的重要性
01-04
父母基因型验证对靶向二代测序结果解读的重要性，何龙霞，杨涛，目的：探索父母基因型验证对靶向二代测序结果解读的重要性。方法：收集21例无家族史非综合征型耳聋先证者（已排除三大常见基因GJB2
二代测序NGS实验方案和应用.doc
10-02
二代测序技术是指第二代高通量测序技术，基于新的测序平台和仪器，能够进行高通量测序，具有广泛的应用，包括全基因组从头测序、目标序列重测序、转录组分析、微生物组研究和基因调控研究等。 二代测序仪器有很多种...
Nature Method：Bioconda解决生物软件安装的烦恼
刘永鑫的博客——宏基因组公众号
08-30 1万+
软件管理器conda 选择一、Anaconda安装 603Mb 选择二、Miniconda安装 38Mb Codna配置 Conda安装软件 指定虚拟环境、python版本和软件版本安装 查看己安装软件的版本和依赖关系 Reference 猜你喜欢 写在后面 其它之前很多人都用过Conda，特别是生信科学家近两年都在用Bioconda，那叫一个妙不可言、根本停不下来。这里...
宏病毒组分析常见的分析软件
E_gene的博客
01-13 3211
病毒是地球上数量最多的生物实体，其中细菌病毒(即噬菌体)约有1031个类群，从海洋到陆地再到人体几乎都是它们的栖息地。研究者将病毒视为调节人类生态系统的重要成员，人体内主要包括真核病毒和噬菌体，包括双链DNA (double-stranded DNA, dsDNA)，单链DNA (single-stranded DNA, ssDNA)和RNA病毒。随着对病毒研究的广泛开展，“病毒组”与“病毒组学”的概念也应运而生，这些术语分别涵盖了栖息在生态系统中的所有病毒及其基因组和对它们的研究(Lefkowitz EJ
Bootstrap Manifest Generator(BMG) 生成的product.xml 试分析
murphey的专栏
12-25 1049
BMG是一个应用于ClickOnce发布中自制系统必备包的图形化工具。ClickOnce的系统必备包放在C:\Program Files\Microsoft SDKs\Windows\v6.0A\Bootstrapper\（VS2008）。通过BMG我们可以自行定义包括第三方组件在内的安装方案。以下是通过比对VS2008自带Packages，对product.xml的试分析。 定
MetaWRAP 安装与配置
HuFeiHu
10-31 6988
MetaWRAP 安装 1、安装Python环境Minconda2 https://repo.anaconda.com/miniconda/Miniconda2-latest-Linux-x86_64.sh bash Miniconda2-latest-Linux-x86_64.sh 2、安装MetaWrap # 配置，按照下面的顺序配置，顺序很重要 conda config --add channels defaults conda config --add channels conda-forge c
bowtie和bowtie2使用条件区别及用法
soyabean555999的博客
03-15 2万+
一、转录组还是基因组？ map常用的工具有bowtie/bowtie2, BWA,SOAP1/SOAP2等。这个问题又会被分成两个问题，是基因组测序（DNA-seq）还是转录组测序(mRNA-seq)。其中的区别是对于真核生物而言，mRNA序列与DNA序列并不完全相同，在经历了后剪切之后，成熟的mRNA可能是原基因的一部分，甚至顺序及个别碱基会产生变化。如果是mRNA测序，那map工作就会在
生信软件 | bowtie2（测序序列与参考序列比对）
热门推荐
Baimoc
04-06 5万+
1. 介绍 2. 安装 2.1下载 2.2 解压 2.3 设置环境变量 3. 使用 3.1 命令 3.1.1 必需参数 3.1.2 可选参数（常用） 3.2 构建索引 3.2.1 官方索引 3.2.2 自建索引 3.3 例子 3.3.1 例子：M. musculus, UCSC mm10 1. 介绍 Bowtie2 是将测序reads与长参考序列比对工具。适...
bowtie和bowtie2用法详解
soyabean555999的博客
03-15 3万+
bowtie 短序列比对工具详解 常见的短序列比对工具有很多，如fasta、blast、bowtie、shrimp、soap等。每个工具都有其自身的优点，但同时也具备了一些缺点。权衡利弊，我选择bowtie作为主要的短序列比对工具。它速度很快，比对结果也容易理解。 现在举个例子来探讨bowtie的使用方法：现在有GENOME.fa、高通量测序数据Reads.fa，我们希望将Reads.fa比对到...
每日一题——Python实现PAT乙级1042 字符统计（举一反三+思想解读+逐步优化）
weixin_44915521的博客
06-09 1105
这些编程思想和哲学理念在编写高效、健壮的代码时至关重要。通过理解和应用这些原则，我们不仅能够写出性能更好的代码，还能在不同情况下做出更加明智的设计决策。通过理解和应用这些哲学和编程思想，可以在不同场景中灵活运用，写出更加高效、健壮和灵活的代码。评估需求：根据实际需求选择最合适的数据结构和算法。避免冗余：尽量避免不必要的重复计算和操作。权衡利弊：在空间和时间之间找到最佳平衡点。灵活应变：根据输入数据的特性，灵活使用动态数据结构和提前决策策略。
基于JavaScript 如何实现爬山算法以及优化方案
最新发布
乐闻世界
06-10 1311
爬山算法（Hill Climbing Algorithm）是一种常见的启发式搜索算法，常用于解决优化问题。其核心思想是从一个初始状态出发，通过逐步选择使目标函数值增大的邻近状态来寻找最优解。接下来，我们将通过 JavaScript 实现一个简单的爬山算法，帮助大家理解其原理和应用。从一个初始状态开始。评估当前状态的目标函数值。在当前状态的邻居中选择一个目标函数值更大的状态。如果找到了更优的邻居，则移动到该邻居并重复步骤2和步骤3。如果没有更优的邻居，则算法结束，当前状态即为局部最优解。
二代测序数据处理序列比对
03-30
二代测序数据处理中的序列比对是将测序产生的短序列（reads）与一个参考基因组进行比对，以确定每个read的来源和位置。序列比对是整个数据分析流程的重要部分，其结果会影响后续的变异检测、转录本定量等分析结果。 ...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值