BETA整合ChIPseq和RNAseq

前些天，生信技能树 表观转录调控之ChIP-seq和RNA-Seq联合分析 介绍了一篇文献取ChIP-seq和RNA-seq数据的交集进行联合分析，小编在底下留言提到了刘Shirley实验室出品的几款整合分析工具，其中有一个BETA软件。本文就此工具做一个使用介绍。

BETA简介

BETA（Binding and Expression Target Analysis）软件整合了ChIP-seq和转录表达水平来探究基因表达调控的机理。此前，一些研究用于转录因子靶基因预测，但是鲜有高效的工具供人们使用。

Home: http://cistrome.dfci.harvard.edu/BETA/ 这个网站可以下载示例的数据

Github: https://github.com/hanfeisun/BETA

文章发表在Nature protocols: https://www.nature.com/articles/nprot.2013.150

作者开发了名为BETA的软件来研究结合位点与转录表达的关系，基于结合位点信息和差异表达信息可以进行如下三个主要分析：

• 用于预测转录因子具有激活还是抑制的功能；
• 推断识别转录因子的直接靶基因；
• 用于鉴定转录因子的motif及其结合者。

一. BETA软件环境配置

BETA的安装建议严格按照以下脚本（不要使用mamba）：

conda create -y -n beta_chip python=2.7.15
conda install -y -c hcc beta 
conda install -y libiconv

个人觉得最难的部分就是软件安装了。

二. 软件参数及运行

对应BETA的三种功能，软件主要有三种子命令：

• BETA Basic: BETA Basic和BETA plus需要2个输入文件:
  • 输入文件一，TF结合位点，即peak的位置；
  • 输入文件二，敲降或敲除或过表达或激活TF后所有基因的变化倍数和P值;
  • BETA Basic能预测转录因子的激活/抑制功能，并识别直接靶基因。
• BETA plus: 除了Basic的两个功能外，还能鉴定转录因子motif及其collaborator，相当于Basic的扩展。
• BETA minus: 只有ChIPseq的数据的情况下，BETA minus可以根据bed文件计算出TF对靶基因的调控潜力。

测试数据在http://cistrome.org/BETA/src/BETA_test_data.zip

1. BETA basic

转录因子激活或抑制功能预测，以及直接靶基因检测。

BETA basic \
    -p 3656_peaks.bed \
    -e AR_diff_expr.xls  \
    -k LIM  \
    -g hg19 \
    --da 500 \
    -n AR \
    -o basic_output_dir 

BETA basic

Ubuntu系统可能会遇到“cut: write error: Broken pipe”，官方回复可以忽略这个报错https://groups.google.com/g/cistromebeta/c/IZ_GttVR6bg/m/2iYPDkRl0U4J

涉及参数及其他可选参数解释如下：

• -p， peaks结合位点文件，BED格式，支持3列或5列格式，其格式为：

 chrom   start   end     name（可选）  score（可选）

使用基本的前3列信息即可。

本例示例的bed文件前三行如下：

$ cat 3656_peaks.bed | head -n 3
chr1    1208689 1209509 AR_LNCaP_2      51.58
chr1    1334246 1335348 AR_LNCaP_7      54.55
chr1    2179351 2180790 AR_LNCaP_9      257.72

• -e，差异表达文件，可以使用LIMMA和Cuffdiff的结果，本例示例是limma的结果，差异表达文件前三行如下：

$ cat AR_diff_expr.xls | head -n 3
#ID     logFC   AveExpr t       P.Value adj.P.Val       B
NR_045762_at    3.16711734      9.140369116     35.91057535     6.99E-11        4.18E-07        14.13456018
NM_001002231_at 3.214550493     9.169929883     35.32505807     8.07E-11        4.18E-07        14.05227211

Gene symbol和DESeq2的结果当然也是可以作为input data，详见下面的参数设置：

• 其他参数：
  • -k 差异表达文件的类型，即LIM, CUF, BSF, O，分别对应LIMMA、cuffdiff标准输出格式、BSF格式或其他格式的文件需要使用--info指定。
  • -g 基因组，支持任何小鼠，hg38, hg19, hg18, mm10, mm9，也可使用-r指定其他基因组文件。
  • --gname2 给定次参数说明文件中的基因名称或转录本名称是官方的基因symbols，默认FALSE。
  • --info 指定表达数据中的geneID, up/down status和statistcal values列。
  • -r 参考基因组文件。
  • -o 输出路径。
  • --bl 是否使用CTCF边界过滤基因周围的peaks，默认FALSE。
  • --bf CTCF保守peakBED文件，需要同时指定--bl参数且基因组不能指定为hg19和mm9。
  • --pn 计划分析的peaks数目，默认10000。
  • --method 支持score和distance，指定TF/CR给你预测的方法，score：调控潜能的打分，disrance：邻近的结合peak的距离。
  • -n 输出文件前缀名称。
  • -d 从TSS一定范围内获取peaks，默认100000（100kb）。
  • --df 输入一个0~1范围内的阈值来筛选最显著的差异表达基因，默认1，全部输入基因。
  • --da 通过比例或数目选取最显著的差异表达基因，默认0.5。如果想使用diff_fdr请设置此参数为1。
  • -c 设置0~1阈值，通过p值（单尾KS检验）选择最近的靶基因，默认0.001。

2. BETA plus

进行激活/抑制功能预测， 直接靶基因预测和靶向区域motif分析，比基础版本增加了motif分析。

BETA plus \
    -p 3656_peaks.bed \
    -e AR_diff_expr.xls  \
    -k LIM \
    -g hg19 \
   --gs ~/reference/genome/hg19/hg19.fa \
    -n AR \
    -o plus_output_dir 

该模块运行时间比较久。

该模块与basic相似，以下参数功能有差异：

• --gs, 基因序列文件。
• --mn, 指定0~1范围内数值作为p-value或大于1作为数目来获取motif的阈值， 默认10。

其他参数参考上述BETA Basic的描述。

三.结果展示

1.结果总览

basic有的，plus都有：

$ tree -h -L 2
.
├── [308K]  3656_peaks.bed
├── [2.9M]  AR_diff_expr.xls
├── [4.0K]  basic_output_dir
│   ├── [118K]  basic_function_prediction.pdf
│   ├── [ 13K]  basic_function_prediction.R
│   ├── [ 47K]  basic_uptarget_associate_peaks.bed
│ └── [ 19K]  basic_uptarget.txt
└── [4.0K]  plus_output_dir
    ├── [111K]  AR_function_prediction.pdf
    ├── [196K]  AR_function_prediction.R
    ├── [330K]  AR_uptarget_associate_peaks.bed
    ├── [177K]  AR_uptarget.txt
    └── [4.0K]  motifresult

2.功能预测结果

结果显示转录因子AR具有显著的激活功能：

3.直接靶基因

top的几个target基因也明显是和AR相关的：

$ cat plus_output_dir/AR_uptarget.txt | head
    Chroms txStart txEnd refseqID rank product Strands GeneSymbol
    chr19 51376688 51383823 NM_001256080 2.192e-07 + KLK2
    chr19 51376688 51383823 NM_005551 2.192e-07 + KLK2
    chr19 51376688 51383823 NR_045762 2.192e-07 + KLK2
    chr19 51376688 51383823 NR_045763 2.192e-07 + KLK2
    chr19 51376688 51383823 NM_001002231 2.192e-07 + KLK2
    chr1 207191865 207206101 NM_001083924 8.847e-07 - C1orf116
    chr1 207191865 207206101 NM_023938 8.847e-07 - C1orf116
    chr21 42836477 42880085 NM_005656 1.024e-06 - TMPRSS2
    chr21 42836477 42879992 NM_001135099 1.032e-06 - TMPRSS2

4.plus的motif结果

相较于basic，plus功能除了生成上/下调靶标文件和靶基因相关的peak文件，还额外生成了靶基因相关motif的结果文件，存放于motifresult文件夹，包含一个*motif.html的motif分析文件。

鉴定到AR相关motif：

Plus motif.html- END -