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
