JZXBK的个人博客分享 http://blog.sciencenet.cn/u/JZXBK

博文

python3 计算 基因组测序结果文件 各碱基数目(个人练习)

已有 2049 次阅读 2020-2-13 20:03 |系统分类:科研笔记

    基因组测学回来的结果后,从assembly(组装)里找到序列文件,格式可能是:.fasta、.fastq、.seq、和.contig。fastq要转化为fasta,转化方法网上一大把哈。我的基因组序列是一个.seq文件:

image.png

$首先先呈上代码

seq=open('SK007.seq','r')

count={}

for i in seq:

    if i[0]=='>':

        continue

    else:

        s=i.split()

        s1=s[0]

        for c in s1:

            count[c] = count.get(c,0) + 1

seq.close()                                                   

for k, v in count.items():

    print(k + ' '+ str(v))

#序列文件的内部样子,共有119个scaffold,即scaffold1-scaffold119

image.png

#首先读取文件,将其赋值给seq,并创造一个空的字典。由于我的序列文件就在我的家目录下,所以读取文件时,可以写作“open('SK007.seq','r')”,如果序列文件在别的路径就需要写出详细的文件路径,例如“open('/data/SK007基因组/Assembly/SK007.seq','r')”,这样python才能找到你的文件。

seq=open('SK007.seq','r')

count={}

image.png

#接着使用‘遍历’的方法读取基因组序列。读取时以一行一行的字符串赋值给i,为了去除">scafford"行,选着条件判断,凡是以‘>’开头的字符串都跳过。其它的字符串通过str.split()方法,将每一个字符串一个一个的放在列表中,去掉末尾的的'\n'(直接用方法str.strip()就能很好的达到目的)。

for i in seq:

    if i[0]=='>':

        continue

    else:

        s=i.split()

        s1=s[0]

        for c in s1:

            count[c] = count.get(c,0) + 1

#如果print(s)话会得到这样一组列表,然后将列表中的每一个字符串取出,对每一个字符串进行一个字母一个字母的'遍历',dict.get(key,default=None),key--字典中要查找的键,如果字典中有该键,则返回相应的值,default--如果没有该键,返回的该参数值。

image.png

#这样就将基因组序列众多碱基序列过了一遍,并利用字典将相应的碱基作为键(key),碱基数目作为值(value),形成了一个字典(dict)。再将字典中的键和值读取出来。由于序列中混有小写字母碱基片段所以结果是这样子。

for k, v in count.items():

    print(k + ' '+ str(v))

image.png



https://wap.sciencenet.cn/blog-3419243-1218367.html

上一篇:α多样性指数变化趋势
下一篇:python3 fasta txt seq contig等纯文本文件的读取 写入

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-12-3 03:17

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部