2024年3月1日更新

推薦直接用最近的新軟件PanDepth，C語(yǔ)言寫(xiě)的，速度很快。

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

以下為舊流程或僅供原理了解：

群體分析一般會(huì)對(duì)所有個(gè)體進(jìn)行深度和覆蓋度的統(tǒng)計(jì)，以便確認(rèn)個(gè)體的測(cè)序質(zhì)量，及時(shí)去除質(zhì)量差的個(gè)體。

計(jì)算深度和覆蓋度的軟件有很多，PanDepth（推薦）、samtools、mosdepth、bedtools、bamdst等等，可以計(jì)算個(gè)體、染色體、某窗口上的深度和覆蓋度。以下文章都有很好的總結(jié)：
測(cè)序數(shù)據(jù)的深度、覆蓋度等計(jì)算 - 簡(jiǎn)書(shū) (jianshu.com)
Bedtools genomecov 計(jì)算覆蓋度 - 簡(jiǎn)書(shū) (jianshu.com)
基因組質(zhì)量評(píng)估：（五）mapping法：7. 用軟件mosdepth統(tǒng)計(jì)BAM文件的深度 - 知乎 (zhihu.com)

但針對(duì)于群體分析，批量計(jì)算所有個(gè)體的深度和覆蓋度，這些軟件則顯得命令繁瑣，本文通過(guò)python腳本一步到位多進(jìn)程同時(shí)計(jì)算深度和覆蓋度，省心省力。

一、生成depth.gz文件

這一步需要排序好、去除pcr重復(fù)的bam文件，通過(guò)samtools生成深度文件。

samtools depth -aa input.bam
#-aa計(jì)算所有位點(diǎn)的深度，包括深度為0的位點(diǎn)
#默認(rèn)輸出input.bam.depth.gz

image.png

二、統(tǒng)計(jì)深度和覆蓋度

1. 原理

深度 = 每個(gè)位點(diǎn)深度相加/有深度的位點(diǎn)總數(shù)
覆蓋度 = 有深度的位點(diǎn)/所有位點(diǎn)

bedtools

2.python多進(jìn)程加速計(jì)算

需要安裝tqdm庫(kù)：pip install tqdm
使用方法：

python3 cal_genome_avg.py filepath out.txt cpu_number
#第一個(gè)參數(shù)為路徑，則會(huì)匹配路徑下該所有depth.gz文件，并計(jì)算深度
python3 cal_genome_avg.py file.depth.gz out.txt cpu_number
#第一個(gè)參數(shù)為depth.gz文件，則會(huì)自動(dòng)計(jì)算該文件深度，追加到out.txt中

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

#1、創(chuàng)建函數(shù)執(zhí)行計(jì)算每個(gè)文件平均深度
#2、主函數(shù)多進(jìn)程

def cal_ave_depth(file):
    """
    輸入：文件名
    輸出：個(gè)體名稱、被覆蓋到的位點(diǎn)的深度、沒(méi)有被覆蓋的位點(diǎn)/所有位點(diǎn)
    """
    counter = 0   # 所有位點(diǎn)的總數(shù) 預(yù)期是基因組的全長(zhǎng)
    dp = 0        # 所有位點(diǎn)的深度
    gap = 0       # 深度為0的位點(diǎn)數(shù)
    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)# 被覆蓋到的位點(diǎn)的深度。按照原來(lái)的計(jì)算方法，需要減去0的地方
    genome_cov = float(gap)/float(counter)  # 沒(méi)有被覆蓋的位點(diǎn)/所有位點(diǎn)
    indiv_name=file.split('/')[-1].split('.')[0]    # 獲取個(gè)體名稱
    j=round(avg,3)                   # 被覆蓋到的位點(diǎn)的深度，保留三位有效數(shù)字 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-沒(méi)有被覆蓋的位點(diǎn)/所有位點(diǎn)
    dp =0
    counter = 0
    print(file+" caluate down!")
    return indiv_name,j,k     # 個(gè)體名稱   被覆蓋到的位點(diǎn)的深度   沒(méi)有被覆蓋的位點(diǎn)/所有位點(diǎn)
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)  # 個(gè)體的名字，深度，覆蓋度
        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:# 每計(jì)算完成一個(gè)，寫(xiě)入outfile
            f.write(str(i)+"\t"+str(j)+"\t"+str(k)+"\n")
    pool.close()
    pool.join()   
    filelist=[]

3.結(jié)果

每列依次為個(gè)體名字、深度、覆蓋度,這個(gè)圖片是第v1版本腳本，需要1-第三列才為真實(shí)的覆蓋度。第三列計(jì)算結(jié)果就是覆蓋度80.4%。

results