用python处理测序数据

一个华大实习的小作业,要求如下:

  1. 统计该Scaffold 序列的Scaffold N50长度。

  2. 按格式输出每一条染色体信息

Chromosome Total Length Effective Length N Region GC Length GC Rate (%)

 

于是刘爷又开始现学现卖……

#-*-coding:utf-8-*-
import sys
#from pettytable import PrettyTable
#打开fasta文件
a = sys.argv[0]
fasta = sys.argv[1]
fa = open(fasta,'r')
#用字典保存文件中的内容,包括名字和序列,名字就是fasta中大于号后面的东西
Fasta_dict = dict()
seq = ''
name = ''
for line in fa.readlines():
    if line.startswith('>'):
        line = line.strip('\n')
        name = line.replace('>','') 
        seq = ''
    else:
        seq += line
    Fasta_dict[name] = seq
#按照要求进行一个排序,按照序列的长度来做一个降序
#找到Scaffold N50
half = 0
total_len = 0
for value in Fasta_dict.values():
    length = len(value)
    total_len = total_len + length
    #print(total_len)


for key,value in sorted(Fasta_dict.items(),key = lambda item:len(item[1]),reverse=True):
    #print('{key}:{value}'.format(key = key,value = len(value)))


    if half <total_len/2:
        half = half +len(value)
    else:
        print("")
        print("Scaffold N50 length is: " + "{value}".format(value = len(value)))
        break
list = [["Chromosome","Total_length","Effective_length","N_region","GC_length","GC_rate"]]
for key,value in sorted(Fasta_dict.items(),key = lambda item:len(item[1]),reverse=True):


    countN = 0
    countGC = 0
    for i in value:
        if i == "N":
            countN = countN +1
    effective_value = len(value) - countN
    for i in value:
        if i == "G":
            countGC += 1
        elif i == "C":
            countGC += 1
    GC_percent = round(countGC / len(value),2)


    list.append([key,len(value),effective_value,countN,countGC,GC_percent])
#打印部分:
for i in list:
    for j in i:
        print("{:^13}".format(str(j)),end='   ')
    print("")
fa.close()

输出效果如下图

这个文件中有一个坑,拼接得到的scaffold中,虽然绝大多数都是
>序列名称及信息
ATCG序列
这样的格式,但是有些片段DNA序列并不是在一行的。在最开始写的版本中,先读>开头的一行,作为字典的key,然后把它的下一行作为value,最终导致我的N50比正确的结果要短一些。于是这个部分更改为
for line in fa.readlines():
    if line.startswith('>'):
        line = line.strip('\n')
        name = line.replace('>','') 
        seq = ''
    else:
        seq += line
    Fasta_dict[name] = seq

这种方式把“>”后面的内容到下一个“>”前的所有行都写到value中,得到正确的结果。

用python处理测序数据》有1个想法

  1. skelviper 文章作者

    另外这个列表输出对齐的事情可真是恶心死我了……最后我也没给他搞对齐了,本来看网上有个包叫prettytable,但是也不知道是作者不维护了还是咋回事,GitHub上搜不到,按网上的教程也不work……如果有什么好的方法请务必告诉我。

发表评论

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据