二代测序通过bam文件统计深度和覆盖度

2024年3月1日更新

推荐直接用最近的新软件PanDepth，C语言写的，速度很快。

#命令示例
pandepth -i input.bam -o test
#软件结果最后一行
##RegionLength: 1416543111  CoveredSite: 1239625885 Coverage(%): 87.51  MeanDepth: 26.74

以下为旧流程或仅供原理了解：

群体分析一般会对所有个体进行深度和覆盖度的统计，以便确认个体的测序质量，及时去除质量差的个体。

计算深度和覆盖度的软件有很多，PanDepth（推荐）、samtools、mosdepth、bedtools、bamdst等等，可以计算个体、染色体、某窗口上的深度和覆盖度。以下文章都有很好的总结：
测序数据的深度、覆盖度等计算 - 简书 (jianshu.com)
Bedtools genomecov 计算覆盖度 - 简书 (jianshu.com)
基因组质量评估：（五）mapping法：7. 用软件mosdepth统计BAM文件的深度 - 知乎 (zhihu.com)

但针对于群体分析，批量计算所有个体的深度和覆盖度，这些软件则显得命令繁琐，本文通过python脚本一步到位多进程同时计算深度和覆盖度，省心省力。

一、生成depth.gz文件

这一步需要排序好、去除pcr重复的bam文件，通过samtools生成深度文件。

samtools depth -aa input.bam
#-aa计算所有位点的深度，包括深度为0的位点
#默认输出input.bam.depth.gz

image.png

二、统计深度和覆盖度

1. 原理

深度 = 每个位点深度相加/有深度的位点总数
覆盖度 = 有深度的位点/所有位点

bedtools

2.python多进程加速计算

需要安装tqdm库：pip install tqdm
使用方法：

python3 cal_genome_avg.py filepath out.txt cpu_number
#第一个参数为路径，则会匹配路径下该所有depth.gz文件，并计算深度
python3 cal_genome_avg.py file.depth.gz out.txt cpu_number
#第一个参数为depth.gz文件，则会自动计算该文件深度，追加到out.txt中

# -*- coding: utf-8 -*-
import multiprocessing as mp
import gzip
import os
from tqdm import tqdm
import sys

#1、创建函数执行计算每个文件平均深度
#2、主函数多进程

def cal_ave_depth(file):
    """
    输入：文件名
    输出：个体名称、被覆盖到的位点的深度、没有被覆盖的位点/所有位点
    """
    counter = 0   # 所有位点的总数 预期是基因组的全长
    dp = 0        # 所有位点的深度
    gap = 0       # 深度为0的位点数
    o = gzip.open(file,'r')
    print("open "+file+" successfully!")
    for line in o:
        site_depth = int(line.rstrip().split()[-1])
        if site_depth == 0:
            gap += 1
        else:
            dp += site_depth
        counter +=1
    o.close()
    avg = float(dp)/float(counter-gap)# 被覆盖到的位点的深度。按照原来的计算方法，需要减去0的地方
    genome_cov = float(gap)/float(counter)  # 没有被覆盖的位点/所有位点
    indiv_name=file.split('/')[-1].split('.')[0]    # 获取个体名称
    j=round(avg,3)                   # 被覆盖到的位点的深度，保留三位有效数字 The round() function returns a floating point number that is a rounded version of the specified number, with the specified number of decimals.
    k=1-round(genome_cov,3)            # 1-没有被覆盖的位点/所有位点
    dp =0
    counter = 0
    print(file+" caluate down!")
    return indiv_name,j,k     # 个体名称   被覆盖到的位点的深度   没有被覆盖的位点/所有位点
if __name__ == "__main__":
    if(len(sys.argv) < 3):
        print("This script is to count average depth.\nUsage:InDir \nAttention: InDir includes all your .depth.gz file.")
        exit(1)
  
    filepath=sys.argv[1]
    outfile=sys.argv[2]
    cpu_num=int(sys.argv[3])
    print(filepath,outfile,cpu_num)
    filelist = []
    if os.path.isfile(filepath):
        print("Your input is only one file.")
        i,j,k = cal_ave_depth(filepath)  # 个体的名字，深度，覆盖度
        with open(outfile,"a+") as f:
            f.write(str(i)+"\t"+str(j)+"\t"+str(k)+"\n")
    elif os.path.exists(filepath):
        for root,curdirs,files in os.walk(filepath):
            for file in files:
                if ("depth.gz" in file):
                    filelist.append(filepath+file)
                else:
                    pass
    else:
        print("input filepath is error! Please check it!")
    #print(len(filelist))
    #print(filelist)
    pool=mp.Pool(cpu_num)
    for indiv_name,j,k in tqdm(pool.imap(cal_ave_depth,filelist)):
        #print(indiv_name)
        with open(outfile,"a+") as f:# 每计算完成一个，写入outfile
            f.write(str(i)+"\t"+str(j)+"\t"+str(k)+"\n")
    pool.close()
    pool.join()   
    filelist=[]

3.结果

每列依次为个体名字、深度、覆盖度,这个图片是第v1版本脚本，需要1-第三列才为真实的覆盖度。第三列计算结果就是覆盖度80.4%。

results