统计 fasta 文件序列长度及 GC 含量

本文介绍了一种使用Python脚本快速计算Fasta格式基因序列中GC含量的方法,并提供了两种实现方式:一种是通过手动处理文件,另一种是利用Biopython库进行序列解析。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

注:该脚本适用于序列不断开的情况

可用一下命令将折行的序列合并为一行

awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' < assembly.fasta | egrep -v '^$' | tr "\t" "\n"  > assembly.fas 

运行脚本

python script.py assembly.fas

import sys

def count_bases(seq):
    counts = {}
    for base in seq:
        if base not in counts:
            counts[base] = 1
        else:
            counts[base] = counts[base] + 1
    return counts

def gc_content(seq):
    counts = count_bases(seq)
    return len(seq), (counts["G"] + counts["C"]) / float(len(seq))

out_file = open(sys.argv[2], 'w')

for line in open(sys.argv[1]):
    line = line.rstrip()
    if line.startswith(">"):
        print >> out_file, line.strip(">"),
    else:
        print >> out_file, "%s %s" % gc_content(line)

out_file.close()

升级版,输入文件是 fasta 格式即可。用 Bio 中的 Seq.IO 解析 fasta 文件,

用 python 的内置函数 count() 的计算速度更快。

import sys
from Bio import SeqIO

out_file = open(sys.argv[2], 'w')

for rec in SeqIO.parse(open(sys.argv[1]), 'fasta'):
    print >> out_file, rec.id, len(rec.seq), (rec.seq.count("C") + rec.seq.count("G")) / float(len(rec.seq))

out_file.close()

转载于:https://www.cnblogs.com/liuhui0622/p/6284577.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值