基因二代测序-分别对blast和bwa比对结果进行统计比较-20230506

最新推荐文章于 2024-07-20 17:12:48 发布
最新推荐文章于 2024-07-20 17:12:48 发布
阅读量441 收藏 7
点赞数 7
文章标签： python 小程序
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/m0_50542529/article/details/135010145
版权

在这里插入图片描述
例子：
20230506_conut_blast.revise.new

##
import pandas as pd
import argparse

parser=argparse.ArgumentParser(description='input error !')
parser.add_argument('-ip1','--inputpath1',required=True,help='input like NL230413-1_S113.blastn.xls')
parser.add_argument('-r','--reference',required=True,help='taxmap_slv_ssu_ref_nr_138.1.xls')
parser.add_argument('-op1','--outputpath1',required=True,help='NL230413-1_S113.blastn.result.xls')
parser.add_argument('-op2','--outputpath2',required=True,help='NL230413-1_S113.blastn.result.genus.xls')
parser.add_argument('-op3','--outputpath3',required=True,help='NL230413-1_S113.blastn.result.species.xls')
# inputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\NL230413-1_S113.blastn.xls'
# reference=r'C:\Users\Administrator\Desktop\自己测试捕获\taxmap_slv_ssu_ref_nr_138.1.xls'
# outputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\NL230413-1_S113.blastn.result.xls'

args = parser.parse_args()
inputpath1=args.inputpath1
reference=args.reference
outputpath1=args.outputpath1
outputpath2=args.outputpath2
outputpath3=args.outputpath3
# df1=pd.read_csv(inputpath1,sep='\t',header=None)

df_ref=pd.read_csv(reference,sep='\t')
dict_accession_to_genus = dict(zip(df_ref['primaryAccession'], df_ref['genus']))
dict_accession_to_species = dict(zip(df_ref['primaryAccession'], df_ref['species']))

dict_need_genus={}
dict_need_species={}
dict_need_genus1={}
dict_need_species1={}
with open(inputpath1,'r') as f2,open(outputpath1,'w') as f3:
    for line in f2:
        line=line.strip('\n')
        list_in=line.split('\t')
        list_in.append(dict_accession_to_genus[list_in[1]])
        list_in.append(dict_accession_to_species[list_in[1]])
        print(line)
        # if list_in[0] in dict_need_genus:
        #     dict_need_genus[list_in[0]].append(list_in[-2])
        # else:
        #     dict_need_genus[list_in[0]]=[list_in[-2]]
        #
        # if list_in[0] in dict_need_species:
        #     dict_need_species[list_in[0]].append(list_in[-1])
        # else:
        #     dict_need_species[list_in[0]]=[list_in[-1]]

        if list_in[0] in dict_need_genus1:
            pass
        else:
            dict_need_genus1[list_in[0]]=list_in[-2]

        if list_in[0] in dict_need_species1:
            pass
        else:
            dict_need_species1[list_in[0]]=list_in[-2]


# for i in dict_need_genus:
#     # dict_need1[i]=max(dict_need[i],key=dict_need[i].count)
#     dict_need_genus1[i] = dict_need_genus[i][0]
#
# for i in dict_need_species:
#     # dict_need1[i]=max(dict_need[i],key=dict_need[i].count)
#     dict_need_species1[i] = dict_need_species[i][0]




##属
genus_value_counts=pd.value_counts(list(dict_need_genus1.values()))
dict_genus_value_counts=dict(zip(genus_value_counts.index,genus_value_counts.values))

df_genus_value_counts=pd.DataFrame.from_dict(dict_genus_value_counts,orient='index',columns=['count'])
df_genus_value_counts=df_genus_value_counts.rename_axis('genus').reset_index()
df_genus_value_counts.to_csv(outputpath2,sep='\t',index=None)



##种
species_value_counts = pd.value_counts(list(dict_need_species1.values()))
dict_species_value_counts = dict(zip(species_value_counts.index, species_value_counts.values))
print(dict_species_value_counts)
df_species_value_counts = pd.DataFrame.from_dict(dict_species_value_counts, orient='index', columns=['count'])
df_species_value_counts = df_species_value_counts.rename_axis('species').reset_index()
print(df_species_value_counts)
df_species_value_counts.to_csv(outputpath3,sep='\t',index=None)

20230508_bwa.10_genome.count

import pandas as pd
import click

@click.command()
@click.option('-ip1','--inputpath1',help='inputpath1(contig.mapped.mapped.sam|..)',required=True)
@click.option('-ip2','--inputpath2',help='inputpath2(NL230413-1_S113.inter.xls|..)',required=True)
# @click.option('-op1','--outputpath1',help='outputpath1(contig.mapped.mapped.inter.xls|..)',required=True)
# @click.option('-op1','--outputpath1',help='outputpath1(genus.count.xls|..)',required=True)
@click.option('-op2','--outputpath2',help='outputpath2(species.count.xls|..)',required=True)

def count_genus_species(inputpath1,inputpath2,outputpath2):

    # inputpath1 = r'C:\Users\Administrator\Desktop\自己测试捕获\contig.mapped.mapped.20000.sam'
    #
    # inputpath2 = r'C:\Users\Administrator\Desktop\自己测试捕获\NL230413-1_S113.inter.xls'
    # # outputpath1 = r'C:\Users\Administrator\Desktop\自己测试捕获\contig.mapped.mapped.20000.inter.xls'
    #
    # outputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\genus.count.xls'
    # outputpath2=r'C:\Users\Administrator\Desktop\自己测试捕获\species.count.xls'

    df2 = pd.read_csv(inputpath2, sep='\t')
    # dict_contig_to_genus = dict(zip(df2['contig_name'], df2['genus']))
    # dict_contig_to_species = dict(zip(df2['contig_name'], df2['species']))
    # print(dict_contig_to_genus)
    dict_contig_to_species = dict(zip(df2['accession'], df2['species']))

    df1=pd.read_csv(inputpath1,sep='\t',quoting=3,usecols=[0,1,2,3,4,5,6,7,8,9,10,11,12],header=None)
    df1.columns=['0','1','2','3','4','5','6','7','8','9','10','11','12']
    # print(df1)
    df1=df1.drop_duplicates(subset=['0'],keep='first').reset_index(drop=True)

    df1['2']=df1['2'].apply(lambda x:x.split(' ')[0])
    df1['6']=df1['6'].apply(lambda x:x.split(' ')[0])

    # list2=[item for item in dict_contig_to_genus]
    # df1=df1[df1['2'].isin(list2)]
    # list2.append('=')
    # df1=df1[df1['6'].isin(list2)]
    #
    # print(df1['6'].tolist())
    #
    ## part1
    df_part1=df1[df1['2']==df1['6']]#.reset_index(drop=True)
    # df_part1['genus']=df_part1['2'].apply(lambda x:dict_contig_to_genus[x.split(' ')[0]])
    df_part1['species']=df_part1['2'].apply(lambda x:dict_contig_to_species[x.split(' ')[0]])
    # df_part1=df_part1[['genus','species']]
    print(df_part1)
    #
    ## part2
    df_part2=df1[df1['6']=='=']#.reset_index(drop=True)
    # df_part2['genus']=df_part2['2'].apply(lambda x:dict_contig_to_genus[x.split(' ')[0]])
    df_part2['species']=df_part2['2'].apply(lambda x:dict_contig_to_species[x.split(' ')[0]])
    # df_part2=df_part2[['genus','species']]
    print(df_part2)
    #
    #
    ## part3
    df_part12=pd.concat([df_part1,df_part2])
    # print(df_part12)
    #
    df_part3=df1[~df1['0'].isin(df_part12['0'])].reset_index(drop=True)
    #
    # # list2=[item+' ' for item in dict_contig_to_genus]
    # df_part3=df_part3[df_part3['6'].isin(list2)]
    #
    #
    # df_part3['genus1']=df_part3['2'].apply(lambda x:dict_contig_to_genus[x.split(' ')[0]])
    # df_part3['genus2']=df_part3['6'].apply(lambda x:dict_contig_to_genus[x.split(' ')[0]])
    df_part3['species1']=df_part3['2'].apply(lambda x:dict_contig_to_species[x.split(' ')[0]])
    df_part3['species2']=df_part3['6'].apply(lambda x:dict_contig_to_species[x.split(' ')[0]])
    print(df_part3)
    # df_part3_1=df_part3[df_part3['genus1']==df_part3['genus2']]
    # df_part3_1['genus']=df_part3_1['genus1']
    #
    df_part3_2=df_part3[df_part3['species1']==df_part3['species2']]
    df_part3_2['species']=df_part3_2['species1']
    # df_part3_2=df_part3_2[['genus','species']]
    #
    #
    df_part3_3=df_part3[df_part3['species1']!=df_part3['species2']]
    df_part3_3['species']=''

    df_part1=df_part1[['species']]
    df_part2=df_part2[['species']]
    df_part3_2=df_part3_2[['species']]
    df_part3_3=df_part3_3[['species']]
    #
    df_last=pd.concat([df_part1,df_part2,df_part3_2,df_part3_3])

    print(df_last)

    df2=df_last.copy(deep=True)
    ## 统计属
    # value1=pd.value_counts(df2['genus'])
    # print(value1.index)
    # dict_value1=dict(zip(value1.index,value1.values))
    # print(dict_value1)
    #
    # df_value1=pd.DataFrame.from_dict(dict_value1,orient='index',columns=['count'])
    # df_value1=df_value1.rename_axis('genus').reset_index()
    # print(df_value1)

    ## 统计种
    value2=pd.value_counts(df2['species'])
    print(value2.index)
    dict_value2=dict(zip(value2.index,value2.values))
    print(dict_value2)

    df_value2=pd.DataFrame.from_dict(dict_value2,orient='index',columns=['count'])
    df_value2=df_value2.rename_axis('species').reset_index()
    print(df_value2)

    # df_value1.to_csv(outputpath1,index=None,sep='\t')
    df_value2.to_csv(outputpath2,index=None,sep='\t')

if __name__=='__main__':
    count_genus_species()

20230509_blast.10_genome.count.revise

##
import pandas as pd
import argparse

parser=argparse.ArgumentParser(description='input error !')
parser.add_argument('-ip1','--inputpath1',required=True,help='input like NL230413-1_S113.blastn.xls')
parser.add_argument('-r','--reference',required=True,help='taxmap_slv_ssu_ref_nr_138.1.xls')
parser.add_argument('-op1','--outputpath1',required=True,help='NL230413-1_S113.blastn.result.xls')
# parser.add_argument('-op2','--outputpath2',required=True,help='NL230413-1_S113.blastn.result.genus.xls')
parser.add_argument('-op3','--outputpath3',required=True,help='NL230413-1_S113.blastn.result.species.xls')
# inputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\NL230413-1_S113.blastn.xls'
# reference=r'C:\Users\Administrator\Desktop\自己测试捕获\taxmap_slv_ssu_ref_nr_138.1.xls'
# outputpath1=r'C:\Users\Administrator\Desktop\自己测试捕获\NL230413-1_S113.blastn.result.xls'

args = parser.parse_args()
inputpath1=args.inputpath1
reference=args.reference
outputpath1=args.outputpath1
# outputpath2=args.outputpath2
outputpath3=args.outputpath3
# df1=pd.read_csv(inputpath1,sep='\t',header=None)

df_ref=pd.read_csv(reference,sep='\t')
# dict_accession_to_genus = dict(zip(df_ref['primaryAccession'], df_ref['genus']))
dict_accession_to_species = dict(zip(df_ref['accession'], df_ref['species']))

# with open(inputpath1,'r') as f2,open(outputpath1,'w') as f3:
#     for line in f2:
#         line=line.strip('\n')
#         list_in=line.split('\t')
#         f3.write(line+'\t'+dict_accession_to_species[list_in[1]]+'\n')

# list_need=[]
dict_need={}
dict_need1={}
with open(inputpath1,'r') as f2,open(outputpath1,'w') as f3:
    for line in f2:
        line=line.strip('\n')
        list_in=line.split('\t')
        list_in.append(dict_accession_to_species[list_in[1]])
        if list_in[0] in dict_need:
            dict_need[list_in[0]].append(list_in[-1])
        else:
            dict_need[list_in[0]]=[list_in[-1]]

# list1=[1,23,4,1,1,4,45,4,4,6]
# print(max(list1,key=list1.count))
for i in dict_need:
    # dict_need1[i]=max(dict_need[i],key=dict_need[i].count)
    dict_need1[i] = dict_need[i][0]

species_value_counts = pd.value_counts(list(dict_need1.values()))
dict_species_value_counts = dict(zip(species_value_counts.index, species_value_counts.values))
print(dict_species_value_counts)

df_species_value_counts = pd.DataFrame.from_dict(dict_species_value_counts, orient='index', columns=['count'])
df_species_value_counts = df_species_value_counts.rename_axis('species').reset_index()
print(df_species_value_counts)
df_species_value_counts.to_csv(outputpath3,sep='\t',index=None)



#         f3.write(line+'\t'+dict_accession_to_species[list_in[1]]+'\n')
#
#
#
#
#
# # df1['genus']=df1[1].apply(lambda x:dict_accession_to_genus[x])
# # df1['species']=df1[1].apply(lambda x:dict_accession_to_species[x])
# df1=pd.read_csv(outputpath1,sep='\t',header=None)
# print(df1)
#
# list_contig=list(set(df1[0].tolist()))
# list_contig.sort(key=list(set(df1[0].tolist())).index)
# print(len(list_contig))
#
# # dict_contig_to_genus={}
# dict_contig_to_species={}
# for item in list_contig:
#     df_part=df1[df1[0]==item].copy()
#     ## 种
#     species_value_counts=pd.value_counts(df_part[12].tolist())
#     dict_species_value_counts=dict(zip(species_value_counts.index,species_value_counts.values))
#     print(dict_species_value_counts)
#     df_species_value_counts=pd.DataFrame.from_dict(dict_species_value_counts,orient='index',columns=['count'])
#     df_species_value_counts=df_species_value_counts.rename_axis('species').reset_index()
#     print(df_species_value_counts)
#     dict_contig_to_species[item]=df_species_value_counts.loc[0,'species']
#
# # print(dict_contig_to_genus)
# print(dict_contig_to_species)
#
# with open(outputpath1,'w') as f1:
#     f1.write('contig_name\tspecies\n')
#     for i in dict_contig_to_species:
#         f1.write(i+'\t'+dict_contig_to_species[i]+'\n')
#
# df_last=pd.read_csv(outputpath1,sep='\t')
#
#
# ##种
# species_value_counts = pd.value_counts(df_last['species'].tolist())
# dict_species_value_counts = dict(zip(species_value_counts.index, species_value_counts.values))
# print(dict_species_value_counts)
# df_species_value_counts = pd.DataFrame.from_dict(dict_species_value_counts, orient='index', columns=['count'])
# df_species_value_counts = df_species_value_counts.rename_axis('species').reset_index()
# print(df_species_value_counts)
# df_species_value_counts.to_csv(outputpath3,sep='\t',index=None)

20230508_匹配基因组的名称和染色体id

import os
inputpath1=r'/data/lijing/assemble_capture_test/species_whole_genome/ZymoBIOMICS.STD.refseq.v2/Genomes'
outputpath1=r'/data/lijing/assemble_capture_test/species_whole_genome/10_species.accession.tax.xls'
list1=os.listdir(inputpath1)

with open(outputpath1,'w') as f2:
    f2.write('accession\tspecies\n')
    for i in list1:
        with open(inputpath1+'/'+i,'r') as f1:
            for line in f1:
                line=line.strip('\n')
                if line.startswith('>'):
                    f2.write('{}'.format(line.split(' ')[0][1:]+'\t'+i.split('_')[0]+' '+i.split('_')[1]+'\n'))

20230509_拆分tax分为genus和species

import pandas as pd

df1=pd.read_csv(r'D:\pycharm\project\20230506_使用\taxonomy.tsv',sep='\t')
print(df1)
df1.columns=['primaryAccession','Taxon']
df1['genus']=df1['Taxon'].apply(lambda x:x.split('g__')[-1].split(';')[0])
df1['species']=df1['Taxon'].apply(lambda x:x.split('s__')[-1].replace('_',' '))
print(df1)

df1.to_csv(r'D:\pycharm\project\20230506_使用\taxonomy.change.tsv',index=None,sep='\t')

# df1=pd.read_csv(r'D:\pycharm\project\20230506_使用\taxonomy.tsv',sep='\t')
# print(df1)
# df1.columns=['accession','Taxon']
# df1['genus']=df1['Taxon'].apply(lambda x:x.split('g__')[-1].split(';')[0])
# df1['species']=df1['Taxon'].apply(lambda x:x.split('s__')[-1].replace('_',' '))
# print(df1)
lj_learn
关注
  • 7
    点赞
  • 踩
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
生信知识与软件小结：BWA、wtdbg
BaldStrong's Log Cabin
10-05 4328
三代组装软件:canu、smartdenovo、wtdbg 三代矫正软件:Racon、Nanopolish Falcon： 一个实验性的二倍体组装工具，测试multi Gb genomes。 Canu ：Celera Assembler的一个分支，专门用于高噪音单分子测序。 blast ：的全称是 basic local alignment search tool 是一种极其常见的序列比对工具。其...
bwa比对软件的使用以及其结果文件（sam）格式说明
weixin_30828379的博客
05-18 2133
一、bwa比对软件的使用 1、对参考基因组构建索引 bwa index -a bwtsw hg19.fa # -a 参数：is[默认] or bwtsw，即bwa构建索引的两种算法，两种算法都是基于BWT的（BWT search while the CIGAR string by Smith-Waterman alignment.）。-a bwtsw对于短的参考序列是不工作的，必须要大于...
BWA比对及Samtools提取目标序列
weixin_56827666的博客
05-13 4917
今天想看一下自己的序列里面会不会有某细菌基因组存在，主要使用BWA和Samtools：bwa主要用于将低差异度的短序列与参考基因组进行比对。主要包含三种比对算法：backtrack、SW和MEM，第一种只支持短序列比对（<100bp），后两种支持长序列比对(70bp~1M)，并支持分割比对（split alignment）。MEM算法是最新的也是官方推荐的。BWA-MEM 是一种新的比对算法，用于将测序 reads 或者组装后 contigs 比对至大型参考基因组，例如人参考基因组。它会自动选择局部比
序列比对-BLAST
wangprince2017
05-19 2万+
一.BWA BWA主要是将reads比对到大型基因组上，主要功能是：序列比对。首先通过BWT(Burrows-Wheeler Transformation，BWT压缩算法)为大型参考基因组建立索引，然后将reads比对到基因组。特点是快速、准确、省内存。由三种类似算法组成： BWA-backtrack，BWA-SW和BWA-MEM。首推BWA-MEM。 三种算法的使用范围 BWA-back...
【bioinfo】bwa mem 比对分值参数测试
lillianqdzpw
06-08 7044
测试bwa BWA mem -M -k10 -T10 ./ref/ref.fa tst1.fq > tst1.sam mem命令比对： 比对分值阈值默认是[30]，即read比对得分值>=30，输出比对结果。 参数 默认 比对情况 分值 说明 -A [1] Match 1 1bp比对得1分 -B [4] Mismatch -4 1bp错配扣4分 -O [6,6] gap(ins,del) -6,-6 1bp的ins扣6分，del扣6分 -E [1] gap e
宏基因组实战7. bwa序列比对, samtools查看, bedtools丰度统计
刘永鑫的博客——宏基因组公众号
11-13 1万+
前情提要如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章 宏基因组分析理论教程 微生物组入门圣经+宏基因组分析实操课程 1背景知识-Shell入门与本地blast实战 2数据质控fastqc, Trimmomatic, MultiQC, khmer 3组装拼接MEGAHIT和评估quast 4基因注释Prokka 5基于Kmer比较数据集sourmash 5基于Kmer比
生信入门：序列比对之blast在线和本地使用
m0_46757774的博客
04-02 2万+
主要内容 1 背景 2 在线blast 3 本地blast 3.1 老版本blast 3.2 新版本blast 背景 序列比对（Sequence Alignment）的基本问题是比较两个或两个以上序列的相似性。 blast作为一种序列相似性比对工具，是生物信息分析最常用的一款软件，必须掌握。不管是做两序列相似性的简单比对，还是引物特异性、序列的来源等个性化分析，都会用到blast比对。许...
行业分类-设备装置-发现新基因的方法和使用的计算机系统平台以及新基因.zip
08-25
3. 生物信息学工具：如BLAST、Bowtie、BWA等，用于序列比对；BEDTools、Samtools等用于处理和分析测序结果；还有GATK、VarScan等用于变异检测。 4. 数据挖掘与机器学习：利用AI和深度学习算法，如卷积神经网络（CNN...
chap二代测序数据分析实用学习教案.pptx
11-22
【二代测序数据分析】是指利用现代生物信息技术对高通量测序产生的大量短序列数据进行处理、解析和解释的过程。这种技术在生物学、医学、遗传学等领域广泛应用，以揭示基因功能、变异信息以及物种间的进化关系等。 ...
ads1-slides:DNA测序算法的幻灯片Coursera类
05-09
3. **比对算法**：比对是将测序读取与参考基因组进行匹配的过程，例如Smith-Waterman算法和Burrows-Wheeler Transform（BWT）为基础的 Bowtie 和 BWA 算法。 4. **变体检测**：比对后，我们需要寻找序列间的差异，...
基因组工程和分析软件_Python_JavaScript_下载.zip
04-04
3. 对齐：将测序数据比对到参考基因组或组装的基因组上，如BWA-MEM、Bowtie2等。 4. 变异检测：检测与参考基因组的差异，包括SNP（单核苷酸多态性）、INDEL（插入/缺失）等，如GATK、FreeBayes。 5. 功能注释：识别...
bwa序列比对
qq_27390023的博客
10-06 4853
BWA主要用于把基因组测序数据比对到参考基因组上，由三种类似算法组成：BWA-backtrack，BWA-SW和BWA-MEM。首推BWA-MEM。 1. 建立索引 bwa index hg38.fa 注：在同一目录下生成以hg38.fa为前缀的5个文件，输入不能是.gz压缩文件。索引建立好了，bwa比对时不再需要原始的hg38.fa 文件。 2.用默认参数比对 # SE （single end） bwa mem hg38.fa test.fq.gz > test_se.sam bw.
生信：2：sam格式文件解读
热门推荐
genome_denovo的博客
12-04 5万+
第二章：生物信息分析第一节：解读sam格式文件1，SAM文件格式介绍SAM（The Sequence Alignment / Map format）格式，即序列比对文件的格式，详细介绍文档：http://samtools.github.io/hts-specs/SAMv1.pdfSAM文件由两部分组成，头部区和主体区，都以tab分列。 头部区：以’@’开始，体现了比对的一些总体信息。比如比对的SA
bwa的使用方法
weixin_34416649的博客
03-05 1059
bwa的使用需要两中输入文件：    Reference genome data（fasta格式 .fa, .fasta, .fna）    Short reads data (fastaq格式 .fastaq, .fq)step 1: 建立 Index根据reference genome data(e.g. reference.fa) 建立 Index File    bwa index -a ...
BWA MEM比对人类参考基因组详解
睿智如水
11-10 3万+
在获得下机数据后，做的第一步是质控。质控工具有很多，这里就不做一一介绍了。这里讲如何使用BWA MEM将质控合格的数据比对到参考基因组上。         BWA是一款基于BWT的快速比对工具，其由三个算法组成。这三个算法分别是：BWA backtrack, BWA SW and BWA MEM。其中，BWA MEM是最新的，其更快更准确，更适合用于人重数据分析。对于上述三种算法，首先需要使用索
bwa mem比对结果错误，sam文件不能被samtools识别的原因之一
fethin的博客
09-21 7626
解决[M::bwa_idx_load_from_disk] read 0 ALT contigs的问题 总结如下，bwa mem比对结果错误，sam文件不能被samtools识别的原因之一是bwa安装的问题！ 写这个初级的帖子，为后来人遇到同样问题的人，在百度搜索的时候能够找到能解决问题的帖子！
使用Python的Turtle库绘制动态风车
最新发布
爱写代码的小朋友
07-20 118
python绘制风车
bwa比对结果怎么统计比对率
01-19
BWA（Burrows-Wheeler Alignment）是一种常用的基因组比对工具，用于将测序数据与参考基因组进行比对。在使用BWA对样本进行比对后，可以根据比对结果来计算比对率。下面是一种常用的统计方法： 首先，从BWA比对输出的文件中提取比对信息。BWA比对结果一般以SAM（Sequence Alignment/Map）格式保存，可以使用SAMtools等相关工具对比对结果进行解析和提取。 从SAM文件中，可以获得每个比对结果的相关信息，如比对的状态（比对成功或失败）、比对上的read数量等。 比对率通常被定义为成功比对的read数量占总read数量的比例。因此，可以通过计算成功比对的read数量除以总read数量得到比对率。 具体的计算步骤如下： 1. 从SAM文件中按行读取每个比对结果。 2. 对于每行结果，判断比对的状态是否为成功（比对状态通常在每行的FLAG列中标注）。 3. 如果比对成功，则成功比对的read数量加1。 4. 继续读取下一行并重复上述步骤，直到读取完所有比对结果。 5. 计算成功比对的read数量除以总read数量，并乘以100%得到比对率（以百分比表示）。 需要注意的是，统计比对率时，应该使用有效的read数量。在有些情况下，比如双端测序，一对read中的其中一条无法比对成功，但另一条可以成功比对，这种情况下，只统计比对成功的read数量。 总之，通过解析BWA比对结果的SAM文件，并统计成功比对的read数量除以总read数量，就可以得到比对率。

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

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

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值