在进行一些特定的功能基因组学数据分析时,我们可能需要对 bedgraph 文件中每个 bin 的值进行一定的平滑操作,以降低随机噪声的影响并提供更好的可视化效果。例如:

  • Repli-seq / BrdU-seq 中量化得到的 Replication Timing
  • OK-seq / Pu-seq 中量化得到的 Replication Fork Directionality

以下是一个用于对 bedgraph 进行 lowess 平滑操作的 python script:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import pandas as pd
import statsmodels.api as sm
import argparse

def get_args():

parser = argparse.ArgumentParser(description='Perform lowess smooth on bedgraph.')

parser.add_argument('--input', '-i', help='Input bedgraph file.', required=True)
parser.add_argument('--output', '-o', help='Output smoothed bedgraph file.', required=True)
parser.add_argument('--span', help='Span size of loess smoothing.', type=int, required=True)
parser.add_argument("--chr", required=True,
help="The chrom to input. e.g. 1-22,X,Y")

args = parser.parse_args()

return args


def parse_range(value):
result = []
for part in value.split(','):
if '-' in part:
start, end = part.split('-')
result.extend(range(int(start), int(end) + 1))
else:
result.append(int(part))
return result


def main():

args = get_args()

data = pd.read_csv(args.input, sep='\t', header=None, names=['chrom', 'start', 'end', 'value'])
data['midpoint'] = (data['start'] + data['end']) / 2

span_size = args.span
smoothed_data = []
chr_list = [f'chr{i}' for i in parse_range(args.chr)]

for chrom in data['chrom'].unique():
chrom_data = data[data['chrom'] == chrom]
if chrom not in chr_list:
continue

chrom_length = chrom_data['end'].iloc[-1] - chrom_data['start'].iloc[0]
span_frac = span_size / chrom_length

lowess = sm.nonparametric.lowess(chrom_data['value'], chrom_data['midpoint'], frac=span_frac)
smoothed_df = pd.DataFrame(lowess, columns=['midpoint', 'smoothed_value'])
chrom_data = chrom_data.reset_index(drop=True)
chrom_data['smoothed_value'] = smoothed_df['smoothed_value']
smoothed_data.append(chrom_data)

smoothed_data = pd.concat(smoothed_data)
smoothed_data[['chrom', 'start', 'end', 'smoothed_value']].to_csv(args.output, sep='\t', header=False, index=False)

if __name__ == '__main__':
main()

下载地址:

https://github.com/JuseTiZ/PyScript-for-CT/blob/main/bedgraph_lowess.py

依赖 module 安装:

1
2
pip install pandas
pip install statsmodels

参数详情:

  • --input / -i,指定需要进行 lowess 平滑操作的 bedgraph 文件,该文件应当仅具有四列,且第一列染色体编号应以 chr 开头。
  • --output / -o,指定输出的平滑后 bedgraph 文件,输出的新 bedgraph 中第四列为 LOWESS smooth 后的值。
  • --span,指定平滑操作时使用的长度,脚本将根据每条染色体的总长确定用于平滑的数据比例。
  • --chr,指定进行平滑操作的染色体编号,可使用 - 指定数字范围,也可使用逗号分隔,例如 1-22,X,Y

应用示例:

假设目前有一个通过 Repli-seq 计算得到的小鼠 RT(Replication Timing) bedgraph 文件 RT.bedgraph,通过以下命令进行 30w bp 的平滑:

1
python bedgraph_lowess.py -i RT.bedgraph -o RT.lowess.bedgraph --span 300000 --chr 1-19

平滑前后 IGV track 示例:

请注意,平滑操作在减少噪音的同时,也损失了部分信息量,因此请根据自己当前使用的数据进行权衡,合理设置 --span 参数,一般情况下:

  • 数据的分辨率越高,该参数指定的值应当越低。反之亦然。
  • 对信息精细程度的要求越高,该参数指定的值应当越低。反之亦然。

例如对于某些 OK-seq 数据而言只需要 60kb 左右的 span size 即可。请根据研究需求或者数据来源文章指定恰当的值。