基因测序结果计算N50

μμ豆奶

已于 2022-09-07 22:31:06 修改

阅读量1.4k

点赞数 2

分类专栏：个人随笔文章标签： python

于 2022-09-05 23:08:01 首次发布

本文链接：https://blog.csdn.net/qq_53666171/article/details/126715454

版权

这篇博客介绍了如何使用Python读取基因测序文件，计算总长度、条数、平均长度、最长序列以及N50等统计指标。通过从大到小排序序列长度，确定不同N值对应的序列比例，最终输出结果。作者强调了在Python环境下运行脚本的方法，并提到将参数提取到外部以提高灵活性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

读取文件

import sys
fn = sys.argv[1]
fi = open(fn,'r')
chr_len = {}
for line in infos:      #还可以写成for line in fi 对句柄进行迭代
  line = line.strip()   #取消数据两端的空格。或者两端的数据。
  infos2 = line.split('\t')  #分割
  chr_len[infos2[0]] =int(infos2[1])

计算

简单信息：总长、条数、平均长度

n= len(chr_len)
total_len = sum(chr_len.values())
mean_len =total_len/n

最长序列

sorted(chr_len.items(),key=lambda x:x[1],reverse=True)
longest_id,longest_len = max(chr_len.items(),key=lambda x:x[1])

N50

思路，我们需要的是列表长度，先要从大到小排序，

l =chr_len.values()  #所有染色体长度
l = sorted(l,reverse=True) #从大到小排序
total_len = sum(l)
n50 = 0
idx = 0  #用于存储变量n的值
cumsum = 0
for i,v in enumerate(l,start=