第一次尝试N50计算，走进生信的第一步

作为个人学习记录，也分享给有需要的uu们！

不理解N50的uu们点这里N50是什么https://blog.csdn.net/u010608296/article/details/102771217?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522163483133716780274188627%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=163483133716780274188627&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduend~default-2-102771217.pc_search_result_hbase_insert&utm_term=N50&spm=1018.2226.3001.4187https://blog.csdn.net/u010608296/article/details/102771217?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522163483133716780274188627%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=163483133716780274188627&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduend~default-2-102771217.pc_search_result_hbase_insert&utm_term=N50&spm=1018.2226.3001.4187

本来是学习Perl的课堂作业，还是尝试了一下换成python，但是python确实没有Perl处理的速度快。小蟒蛇输给大骆驼了 QAQ

写的过程也超级曲折，不断修改完善代码，熬了个大夜！(小白哭哭)

##Author：Griffy
##Date:2021-10-21

##Version1:①使用正则，没必要用；②不能够巧妙地使用字典，如何同时遍历字典的键和值
##Version2：①改正了对基因组的理解；②边读文件边处理速度更快(处理大文本的正确方法)
##Version3:设置了函数
##Version4:修改了列表创建方式seq_length_list
##Version5：修改失败，速度更慢
##Version6：将header进行排序，方便看出染色体DNA；②使用time.perf_counter()计算程序运行时间。如果用time()则会受其他程序的干扰
import time
def N50(list_):
    add=0
    list_.sort(reverse=True)
    for length in list_:
        add+=length
        if add >= total_length*0.5:
            print("N50 =",length)
            break

start=time.perf_counter()
seq_dict={}
total_length=0
with open("griffy.fa","r") as f:   ##这里处理的是水稻基因组文件
    for line in f:
        line=line.strip("\n")
        if ">" in line:
            header=line
            seq_dict[header]=0
        else:
            seq_dict[header]+=len(line)
            total_length+=len(line)
head_list=[h for h in seq_dict]
head_list.sort()
for d in head_list:
    print(d,"\t",seq_dict[d])

seq_length_list=[seq_dict[d] for d in seq_dict]
N50(seq_length_list)
end=time.perf_counter()
print(end-start,"s")