基因测序结果计算N50

读取文件

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=1):  #可以帮我们多迭代出当前循环的第几个元素
  cumsum += v              # 累加               
  if cumsum >= total_len*0.5:
    idx = i
    n50 = v
    break     #结束所有循环
 
print('N50 = {}，n={}' .format(n50,idx))
  
  
def Nx0(l,x):
  '''
  该函数用于计算Nx0，接受两个参数：
  l : 长度列表
  x : N50|N60|N70...
  返回两个值：
  idx : Nx0 的编号
  Nx0 : Nx0 的长度
  '''
  l = chr_len.values()  #所有染色体长度
  l = sorted(l,reverse=True) #从大到小排序
  p = int (x[1:]）/100     #x[1:]去掉Nx0中的N，int()转换为数字
  total_len = sum(l)
  nx0 = idx =cumsum =0
  for i,v in enumerate(l,start=1):  #可以帮我们多迭代出当前循环的第几个元素
    cumsum += v              # 累加               
    if cumsum >= total_len * p:
      idx = i
      nx0 = v
    break
       

  return  idx, nx0

Nx0(chr_len.values(),'N100')

打印输出

fn = ‘Python20201122/data/TAIR10/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.fai’
fi_name = fn.split(‘/’)[-1]
print(‘stats for {}’.format(fi_name))
print(‘sum = {},n ={},ave ={},largest ={}’.format(total_len,n,mean_len,longest_id))

for n in (‘N50’, ‘N60’, ‘N70’, ‘N80’, ‘N90’, ‘N100’):
idx,nx0 = Nx0(chr_len.values(),n)
print(‘{}={},n={}’.format(n,nx0,idx))

运行程序

在Python路径来运行

直接点Source Script

把参数提取到外部

import sys

print(sys.argv) #列表 程序后的命令行存写运行程序的时候，把文件路径加到程序后就ok了。

小结

今天好累，加油，努力，明天也要充满干劲！