利用Python脚本计算基因组测序数据Nx0

    import sys
    
    fn = sys.argv[1]
    fi = open(fn,'r')
    
    chr_len = {}
    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
      
    longest_id,longest_len = max(chr_len.items(),key=lambda x:x[1])
    
    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 = 'C:/Users/ASUS/Documents/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))

    C:\Users\ASUS\AppData\Local\Programs\Python\Python310>C:/Users/ASUS/Documents/Python20201122/test.py  C:/Users/ASUS/Documents/Python20201122/data/TAIR10/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.fai
    stats for Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.fai
    sum = 119667750,n =7,ave =17095392.85714286,largest =1
    N50=23459830,n=3
    N60=23459830,n=3
    N70=19698289,n=4
    N80=19698289,n=4
    N90=18585056,n=5
    N100=154478,n=7