生物信息系列：根据成对bed文件以及多个BigWig文件画出信号

引言：本期文章是记录一下如何利用一对bed文件以及多个BigWig文件提取信号值并进行画图的实现过程。

Materials：两个bed文件，多个ENCODE数据库下载的组蛋白CHip-seq数据的bigwig文件。

bigwig:

image

两个bed以及部分内容:

image

Method：bwtool

shell文件：

#! /bin/bash
cd /media/zy/Elements/组蛋白/
a='/home/zy/p53_enh.txt'
b='/home/zy/no_p53_enh.txt'
d='_p53.txt'
e='_noP53.txt'
k='_sig.txt'
c=0
f=1
for i in $(ls);
do
name1=${i%.*}${d}
name2=${i%.*}${e}
bwtool extract bed $a $i $name1
bwtool extract bed $b $i $name2
filelist[$c]=$name1
filelist[$f]=$name2
let c=c+2
let f=f+2
done
for y in ${filelist[*]};
do
cat $y |awk '{print $5}' |tr ',' '\t' > ${y%.*}${k}
done

思路：bwtool这个工具的extract参数可以将bed对应范围内的信号值提取出来。

对于生成的文件我们都放进了一个列表里，然后利用for对每一项进行了awk提取特定行操作。

难点：shell和python的写法可谓大相径庭，像字符串的切割，循环等等

运行完成得到的结果：

image

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

python文件：

# -*- coding:utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pandas import DataFrame
import os
from tqdm import tqdm
class Picture():
    def __init__(self,a,b):
        self.path1="/media/zy/Elements/张茵/组蛋白"
        self.path2 = "/media/zy/Elements/张茵/组蛋白/最后的图"
        self.a = os.path.join(path,a)
        self.b = os.path.join(path,b)
        self.nameHeader=a.split('_')[0]
        self.savePic=os.path.join(self.path2,self.nameHeader)
        self.picName1=a.split('_')[0]+'_'+a.split('_')[1]
        self.picName2=b.split('_')[0]+'_'+b.split('_')[1]
        self.signalPic(self.a,self.b)
    def signalPic(self,x,y):
        df1=pd.read_table(x,header=None)
        df2 =pd.read_table(y, header=None)
        v1=df1.mean(axis=0)
        v2=df2.mean(axis=0)
        v1.to_csv('/home/zy/1.txt',index=False,header=False,sep='\t')
        v2.to_csv('/home/zy/2.txt', index=False, header=False, sep='\t')
        chunk1=pd.read_table('/home/zy/1.txt',header=None,chunksize=10)
        chunk2 = pd.read_table('/home/zy/2.txt', header=None, chunksize=10)
        l1=[]
        l2=[]
        for i in chunk1:
            l1.append(i.mean())
        for y in chunk2:
            l2.append(y.mean())
        self.showPic(l1,l2)
        # 下面是画图操作
    def showPic(self,p1,p2):
        d1 = list(range(-100, 100))
        fig, ax = plt.subplots()
        ax.plot(d1, p2, color="blue", label=self.picName2)
        ax.plot(d1, p1, color="IndianRed", label=self.picName1)
        q = [-100, 0, 100]
        ax.set_xticks(q)
        ax.set_xticklabels(labels=['-1kb', 'center', '1kb'])
        ax.spines['right'].set_color('none')
        ax.spines['top'].set_color('none')
        ax.grid(axis='y', linestyle="-.")
        ax.legend(loc="best")
        ax.set_title(self.nameHeader+'_VALUE')
        if not os.path.exists(self.path2):
            os.makedirs(self.path2)
        plt.savefig(self.savePic+'.png')
if __name__=='__main__':
    path="/media/zy/Elements/张茵/组蛋白"
    li=[]
    li_head=[]
    # 对路径下的文件筛选出sig.txt结尾的文件
    for i in os.listdir(path):
        if i[-7:]=="sig.txt":
            li.append(i)
    # 把筛选出来的文件的‘头’构建成一个新列表
    for y in li:
        y1=y.split('_')
        li_head.append(y1[0])
    # 利用set去重
    l1=list(set(li_head))
    # 利用‘头’列表把头一样的文件筛选出来
    for i in tqdm(l1):
        for y in li:
            if i==y.split('_')[0]:
                if 'no' in y:
                    a1=y
                else:
                    a2=y
                    Picture(a1,a2)

思路：python主要是进行了数据整理和画图，首先要把我们得到的信号文件传入列表，对列表里的元素要进行两两配对（例如：a_p53和a_noP53为一对）,因为两个文件画一张图。

难点：在对数据整理的时候用到了chunksize这个参数，是用来对文件切块，上手快效果好。（具体：数据处理是先每列求均值，得到的结果再每隔十个取均值）

－－－－－－－－－－－－－－－－－－－－－－－－－－

结果展示：

image image