引言：本期文章是记录一下如何利用一对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.path2 = "/media/zy/Elements/zy/组蛋白/最后的图"
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/zy/组蛋白"
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):
a1=str(i)+'_noP53_sig.txt'
a2=str(i)+'_p53_sig.txt'
Picture(a1,a2)
思路：python主要是进行了数据整理和画图，首先要把我们得到的信号文件传入列表，对列表里的元素要进行两两配对（例如：a_p53和a_noP53为一对）,因为两个文件画一张图。
难点：在对数据整理的时候用到了chunksize这个参数，是用来对文件切块，上手快效果好。（具体：数据处理是先每列求均值，得到的结果再每隔十个取均值）
－－－－－－－－－－－－－－－－－－－－－－－－－－
结果展示：
image
image
