用python进行基因组组装结果统计

可以说写的是非常蹩脚了，后面再慢慢改进。。。。

#!/usr/bin/python
# -*- coding: utf-8 -*-
# 需要先conda activate python3, 一定要是python3版本！！！
"""
    作者：徐诗芬
    内容：统计基因组序列特征:
    Total，Longest，Shortest, Mean, Length>=1kb， Length>=2kb, Length>=5kb, N50, N60, N70, N80, N90
    日期：2020.11.20
    版本信息：适用于单行或多行的序列数据
"""
import sys
def usage():
    print('Usage: python genome_stat.py [input_genome.fa] [genome_stat]')


def read_fasta(input_file):
    #  with open(input_file, 'r') as f:  # with open 相当于读取一个文件后并关闭了，因此后面的代码模块要在with里面
    fasta = {}  # 定义一个空的字典：存放序列数据，key为序列名，value为序列
    for line in input_file:
        line = line.strip()  # 去除末尾换行符
        if line.startswith('>'):
            line = line.strip()
            header = line[1:]
            # name.append(header)  #
        else:
            sequence = line
            fasta[header] = fasta.get(header, '') + sequence
            # header是一个变量表示key, .get()获得当前key的值,
            # 此处一开始键hearder是没有值的，赋予一个'',如果fasta.get(header)会报错
    return fasta

def seq_length(fasta):
    seq_len = []
    for seq_i in list(fasta.values()):
        a = len(seq_i)
        seq_len.append(a)
    return seq_len


def seq_len_nkb(seq_len, size):
    """
    根据阈值求contig_number和contig_length
    """
    seq_len_kb = []
    for x in seq_len:
        if x >= size:
            seq_len_kb.append(x)
    num_seq = len(seq_len_kb)
    sum_len = sum(seq_len_kb)
    return num_seq, sum_len

def seq_len_N(seq_len, size):
    """
        根据阈值求N50, N60, N70, N80的contig_number和contig_length
    """

    seq_len_sort = sorted(seq_len, reverse=True)
    sl = []
    for y in seq_len_sort:
        sl.append(y)  # 循环体，循环添加符合条件的序列到列表中
        if sum(sl) / sum(seq_len) * 100 > size:
            break  # 设置跳出循环的条件
        # m = y
    num_seq = len(sl)
    return num_seq, y



def main():

    input_file = open(sys.argv[1], 'rt')
    output_file = open(sys.argv[2], 'wt')
    #output_file = open("genome_stat.txt", 'wt')
    fasta = read_fasta(input_file)  # 读取fasta文件
    # print(fasta)
    name = list(fasta.keys())  # 提取字典里的键（序列名），生成list
    # print(name)
    # seq = list(fasta.values())  # 提取字典里的值（序列），生成list
    # print(seq)
    # GC_content =
    seq_len = seq_length(fasta)
    # print(seq_len)
    # seq_len_sort = sorted(seq_len, reverse=True)
    # print(seq_len_sort)

    # 根据列表创建字典，用于储存对应的各个序列长度数据
    # sequence_distribution = dict(zip(name, seq_len))
    # print('序列长度分布为：{}'.format(sequence_distribution))
    # print('**************************************************')
    # 序列长度的统计1: Total，Longest，Shortest, Mean
    seq_total = sum(seq_len)
    total_num = len(name)
    # print('序列总长度：{:,} 序列总数目：{}'.format(seq_total, total_num))

    seq_longest = max(seq_len)
    # print('最长的序列长度：{:,}'.format(seq_longest))
    seq_shortest = min(seq_len)
    # print('最短的序列长度：{:,}'.format(seq_shortest))
    seq_mean = sum(seq_len) / len(name)
    # print('平均序列长度：{:,}'.format(int(seq_mean)))
    # print('**************************************************')
    # 序列长度的统计2：Length>=1kb， Length>=2kb, Length>=5kb
    seq_len_1kb = seq_len_nkb(seq_len, 1000)
    seq_len_2kb = seq_len_nkb(seq_len, 2000)
    seq_len_5kb = seq_len_nkb(seq_len, 5000)
    # print('Length>=1kb的序列总长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_1kb))
    # print('Length>=2kb的序列总长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_2kb))
    # print('Length>=5kb的序列总长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_5kb))
    # print('**************************************************')
    # 序列长度的统计3：N50, N60, N70, N80, N90
    seq_len_N50 = seq_len_N(seq_len, 50)
    seq_len_N60 = seq_len_N(seq_len, 60)
    seq_len_N70 = seq_len_N(seq_len, 70)
    seq_len_N80 = seq_len_N(seq_len, 80)
    seq_len_N90 = seq_len_N(seq_len, 90)
    # print('N50的序列长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_N50))
    # print('N60的序列长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_N60))
    # print('N70的序列长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_N70))
    # print('N80的序列长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_N80))
    # print('N90的序列长度：{0[1]:,} 序列数：{0[0]}'.format(seq_len_N90))

    output_file.write('StatType\tContigLength\tContigNumber\n')
    output_file.write('Total\t{0:,}\t{1}\n'.format(seq_total, total_num))
    output_file.write('Longest\t{0:,}\t{1}\n'.format(seq_longest, 1))
    output_file.write('Shortest\t{0:,}\t{1}\n'.format(seq_shortest, 1))
    output_file.write('Mean\t{0:,.0f}\t{1}\n'.format(seq_mean, 1))
    output_file.write('Length>=1kb\t{0[1]:,}\t{0[0]}\n'.format(seq_len_1kb))
    output_file.write('Length>=2kb\t{0[1]:,}\t{0[0]}\n'.format(seq_len_2kb))
    output_file.write('Length>=5kb\t{0[1]:,}\t{0[0]}\n'.format(seq_len_5kb))
    output_file.write('N50\t{0[1]:,}\t{0[0]}\n'.format(seq_len_N50))
    output_file.write('N60\t{0[1]:,}\t{0[0]}\n'.format(seq_len_N60))
    output_file.write('N70\t{0[1]:,}\t{0[0]}\n'.format(seq_len_N70))
    output_file.write('N80\t{0[1]:,}\t{0[0]}\n'.format(seq_len_N80))
    output_file.write('N90\t{0[1]:,}\t{0[0]}\n'.format(seq_len_N90))

    input_file.close()
    output_file.close()


try:
    main()
except IndexError:
    usage()