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

最新推荐文章于 2025-03-28 17:20:59 发布

转载最新推荐文章于 2025-03-28 17:20:59 发布 · 2.6k 阅读

7 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/liuhui0622/p/6284577.html

文章标签：

#python #awk #r语言

本文介绍了一种使用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