根据bam文件画bin depth统计测序深度-2020-01-

根据bed文件及bam文件，得到每个bin的测序深度：

一、导入包：

import os, subprocess

import numpy as np

import pandas as pd

import matplotlib.pylab as plt

import matplotlib

二、准备文件，更改工作路径：

os.chdir({where you put your data})

bed = {your bed file}

bam = {your bam file}

SampleName = bam.strip(".bam")

三、bedtools统计测序深度depth,输出格式为pandas的dataframe：

cmd = f"bedtools coverage -a {bed} -b {bam} | awk '{{print $1,$2,$3,$4}}'"

with os.popen(cmd) as BedtoolsOut:

BedDF = pd.read_table(BedtoolsOut, delimiter=" ", names=["chrom", "start", "end", "depth"])

BedDF.sort_values(by=['chrom','start'],axis=0,ascending=False,)

四、准备画图参数

1. chrom列表：

ChromNameLis = list(pd.Series(BedDF.iloc[:,0]).unique())

2. 图注文字：

describe = BedDF["depth"].describe()

addtextstr = "\n".join(str(describe).replace(".000000","").split("\n")[:-1])

3. 颜色（多条染色体的画，用不同的颜色区分开）：

color = ['lightskyblue', "yellowgreen"]

4. 字体

matplotlib.rcParams['font.family'] = 'SimHei'

5. 图片设置：

plt.figure(figsize=(20,6))

plt.title(SampleName)

plt.xlabel("chrom")

plt.ylabel("depth")

plt.xticks([])

ax = plt.subplot(1,1,1)

ax.spines['top'].set_visible(False)

ax.spines['right'].set_visible(False)

ax.text(X.max(),BedDF["depth"].max(),addtextstr,fontsize=8,horizontalalignment='left',verticalalignment='top')

注：右侧和上方的边框去除；右上角添加文字说明——pandas对depth列的describe；

五、根据染色体来画图：

xticksArr = [0]

n=0

for index, chrom in enumerate(ChromNameLis):

    Depthlis = BedDF[BedDF.iloc[:,0]==chrom]["depth"]

    meanvalue = int(Depthlis.mean())

    X = np.arange(n,n+len(Depthlis))+1

    Y = np.array(Depthlis)

    Y = Y.astype(int)

    plt.bar(X,Y,width = 1,facecolor = color[index%2],edgecolor = color[index%2])

    plt.text(xticksArr[-1]+0.5*len(Depthlis), max(Depthlis), meanvalue,     fontsize=8,horizontalalignment='center',verticalalignment='bottom' )

    plt.text(xticksArr[-1]+0.5*len(Depthlis),-300, chrom,     fontsize=8,horizontalalignment='center',verticalalignment='top' )

    n += len(Depthlis)

    n += 10

    xticksArr.append(xticksArr[-1]+len(Depthlis))

    xticksArr.append(xticksArr[-1]+10)

注：染色体间隔设置为10；在上方加上depth的均值；下方标注上染色体的编号；

六、图片保存为文件

plt.savefig(f"{SampleName}.png")

RD-CF1578.sorted.dedup.png