生信算法6 - 比对reads碱基数量统计及百分比统计

比对reads碱基数量统计及百分比统计

输出bam文件路径，需要安装samtools

from addict import Dict
import collections
import os

def count_base(bam_path: str):
    read_bam = os.popen(f"samtools view {bam_path}|less -S").readlines()

    sequence, dict_base = '', Dict()
    total_base = dict_base['A'] = dict_base['T'] = dict_base['C'] = dict_base['G'] = dict_base['N']= 0
    for idx, line in enumerate(read_bam):
        sequence = line.strip().split('\t')[9]
        # reads碱基计数
        dict_seq = collections.Counter(sequence)
        
        dict_base['A'] += dict_seq.get('A', 0)
        dict_base['T'] += dict_seq.get('T', 0)
        dict_base['C'] += dict_seq.get('C', 0)
        dict_base['G'] += dict_seq.get('G', 0)
        dict_base['N'] += dict_seq.get('N', 0)
        total_base += len(sequence)

    # 碱基百分比
    dict_base['a_pct'] = round(int(dict_base['A']) / total_base * 100, 4)
    dict_base['t_pct'] = round(int(dict_base['T']) / total_base * 100, 4)
    dict_base['c_pct'] = round(int(dict_base['C']) / total_base * 100, 4)
    dict_base['g_pct'] = round(int(dict_base['G']) / total_base * 100, 4)
    dict_base['n_pct'] = round(int(dict_base['N']) / total_base * 100, 4)
    dict_base['total_base'] = total_base

    # 整理输出字典
    dict_statistics = Dict()
    dict_statistics['a_stat'] = '{0}({1}%)'.format(dict_base['A'], dict_base['a_pct'])
    dict_statistics['t_stat'] = '{0}({1}%)'.format(dict_base['T'], dict_base['t_pct'])
    dict_statistics['c_stat'] = '{0}({1}%)'.format(dict_base['C'], dict_base['c_pct'])
    dict_statistics['g_stat'] = '{0}({1}%)'.format(dict_base['G'], dict_base['g_pct'])
    dict_statistics['n_stat'] = '{0}({1}%)'.format(dict_base['N'], dict_base['n_pct'])
    dict_statistics['total_base'] = total_base
        
    return dict_statistics

dict_statistics = count_base(bam_path="/path/sample.bam")
print(dict_statistics)

# 输出结果
# {'a_stat': '9592(30.6816%)', 't_stat': '8675(27.7485%)', 'c_stat': '8435(26.9808%)', 'g_stat': # '4554(14.5667%)', 'n_stat': '7(0.0224%)', 'total_base': 31263}