在进行一些特定的功能基因组学数据分析时，我们可能需要对 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 即可。请根据研究需求或者数据来源文章指定恰当的值。