5.3. CLIP-seq (RNA-Protein Interaction)
Last updated
Last updated
RBP(RNA-binding protein)在RNA调控中发挥了重要的作用。CLIP(CrossLinking and ImmunoPrecipitation)-seq是一种高通量试验方法,用于鉴定特定RBP在转录组上的结合位点。CLIP-seq的做法是用紫外光对RBP和RNA进行共价交联,对所关注的RBP用抗体富集后,对RBP-RNA complex中的RNA组分建库测序。这样我们就能用和CHIP-seq数据分析比较类似的方法确定我们所关注的RBP在RNA上的结合位点。
CLIP-seq实验有很多的变种,其中有一种被称为PAR-CLIP (PhotoActivatable-Ribonucleoside-enhanced CrossLinking and ImmunoPrecipitation)。在PAR-CLIP实验中,对于培养的细胞系, 会在培养基中加入photoactivatable ribonucleotides,如4-thiouridine(4-SU)。4-SUs在转录过程中会被掺入细胞内的RNA中。同样经过紫外处理后,在建库的过程中,和蛋白发生交联4-SU会在逆转录中被识别为胞嘧啶,最终在交联的位点引入T->C的转变。 这样除了read coverage的信号,T->C conversion的信号也可以为确定RBP结合位点提供有用的信息。
CTK是Chaolin Zhang lab开发的一个处理CLIP-seq数据的工具,支持多种CLIP-seq实验的数据分析。在本章中,我们以PAR-CLIP为例介绍如何使用CTK处理CLIP的数据。
本章主要参考CTK文档中的PARCLIP data analysis using CTK一节。有兴趣的同学可以直接阅读文档。
CLIP-seq在数据分析上和我们前面介绍的CHIP-seq非常相似,同样是对reads进行预处理和mapping后,再通过peak calling确定所关注的调控因子的结合位点。
在前面的章节中我们已经对mapping进行了讨论,所以这里提供了mapping的结果,大家可以直接从sam文件出发进行后续分析。
clip-seq实验的主要目的是确定RBP的结合位点,这和CHIP-seq数据分析一样,也是通过peak calling来实现的。
二代测序在PCR扩增的过程中常常会产生一些duplicates,这些duplicate序列基本是相同的,所以会mapping到基因组上的同一个位置,这就很容易产生一些假阳性的peak。为了降低duplication导致的假阳性,CTK进行了非常严格的deduplication,具体做法是采用下面两个策略:
在数据预处理阶段,就对reads进行collapsing,对于完全identical的reads,只保留其中的一个。
CTK提供了一个tag2collapse.pl
脚本。即使不是完全identical的reads,如果它们mapping到基因组上的起始位置相同,这个脚本也认为它们是duplication,只保留其中的一个。
方法1: 使用docker
install software (already available in Docker,docker images的下载链接如附表所示): CTK, bwa, fastx_toolkit, samtools
Get data (already available in Docker),我们使用CTK tutorial中的两个PAR-CLIP样本:
SRR048967, SRR048968: HEK293 cells performed PAR-CLIP with 4-SU to PUM2 (pumilio2)
此处我们提供的是经过mapping的sam文件
SRR048967.sam.gz
和SRR048967.sam.gz
,位于/home/test/files/input/
目录下。为了提高运行速度,我们只选取了原始测序数据的前1,000,000条数据,如果感兴趣的话,也可以通过点击对应链接下载对应的fastq文件从头开始,建立mapping index所需的genome fasta文件可以从清华云下载:download link
方法2: 直接下载软件
Install software:
CTK: download link
bwa: download link
fastx_toolkit: download link
samtools: download link
Get data: download link
(这段加在6. RNA regulation里面)
(6) 6.10 CLIP-seq: bioinfo_clip_seq.tar.gz
首先进入到容器:
以下步骤都在/home/test/CLIP-seq
目录下进行:
CTK先用parseAlignment.pl
脚本从sam文件中提取出reads的位置信息和突变信息,方便后续分析。
在CTK流程中,parseAlignment.pl这里输入的是gzip压缩后的sam文件,主要是因为sam文件比较占空间。在实际中我们更常见的做法是使用bam格式存储。
-v
是"verbose"的意思,即输出更详细的日志信息
--map-qual 1
:只考虑mapping quality > 1的reads
--min-len 18
:只考虑长度大于等于18的reads
--mutation-file
:将突变信息存储在哪个文件中
这一步我们的输入是sam.gz文件,输出是记录read mapping到的位置的${i}.tag.bed
文件和记录突变信息的${i}.mutations.txt
文件。
在CLIP-seq和CHIP-seq中,"tag"是"read"的同义词。
接下来我们使用tag2collapse.pl
,在起始位置相同的reads中只保留一个用于后续分析:
其中参数的意义如下:
-v
:同上
-big
:输入文件较大(采用特定算法)
-weight
:在合并时考虑每一个序列的权重
--weight-in-name
:序列权重在bed文件的name这一列中
--keep-max-score
:在合并时保留权重最高的序列(默认保留长度最长的序列)
--keep-tag-name
:保留bed文件中每一个序列原有的name
合并tag.bed中的重复序列之后,对于保留下来的reads,我们用joinWrapper.py
这个脚本从mutations.txt
选出每一个reads对应的突变信息:
joinWrapper.py
实际上是用python对gnu的join
命令进行了一个封装。在mutations.txt和uniq.bed文件中,第4列都对应着read id。 其中4 4代表使用两个文件第四列的read id进行连接,而N表示只会输出第4列相符的行。
接下来我们通过peak calling推断RBP的结合位点:
这个命令的输入是tag.uniq.bed
文件,即去过duplication后的reads mapping到基因组上的位置。tag2peak.pl
只用到了read coverage的信息,并没有用到突变的信息。tag.uniq.peak.sig.bed
是输出的peaks的对应的genomic intervals。
其他参数含义如下:
-big
:输入文件较大(采用特定算法)
-ss
:将正负链分开考虑
-v
:输出更详细的信息
--valley-seeking
:通过valley seeking的方式来确定具体的结合区域
-p
:p value cutoff
--valley-depth
:要求的峰与峰之间的“谷”相对于峰的高度
--multi-test
:使用Bonferroni方法对多次检测的p值进行修正
--dbkey
:基因组和版本。CTK内部自带了dm6,hg19,hg38,mm10四个基因组的基因注释,如果提供了这个参数,tag2peak.pl
会把peak注释到相应的基因上。
--out-boundary
和--out-half-PH
都是可选的参数,可以获得一些额外的输出:
--out-boundary
:输出峰的边界
--out-half-PH
:输出峰一半高度处的边界。
我们知道PAR-CLIP会在RNA和RBP cross-linking的位点引入T->C的转变,这种conversion也可以为确定RNA-RBP的相互作用提供有用的信息。CTK提供了一种他们称之为CIMS(cross-link-induced mutation sites)的分析来利用PAR-CLIP实验中T->C conversion提供的信息,确定哪些位置是比较可信的cross linking位点。
首先,我们从去过duplication的mutation中挑选出T->C的conversion:
然后进行CIMS分析:
这个脚本的输入是去过duplication的read位置信息tag.uniq.bed
,T->C conversion的信息tag.uniq.t2c.bed
,输出是推断出的cross linking位点和tag.uniq.t2c.CIMS.txt
和mutation位点的统计信息tag.uniq.t2c.pos.stat.txt
。
其他参数的意义如下:
outp
: output mutation position summary
-big
:输入文件较大(采用特定算法)
-n
:迭代次数
-p
:记录突变位点相对于read起始位点的距离
-v
:输出更详细的信息
CIMS.pl
输出中,每行对应一个推断出的cross linking的位点,第9列提供了FDR值。我们可以根据FDR对这些cross linking的位点进行筛选:
前面的介绍中,我们是直接从mapping的结果出发的,以下介绍了mapping以及mapping之前的数据预处理步骤供大家参考。
CLIP-seq测的是RNA序列,理论上用STAR/hisat2这类spliced aligner会更合适,但是由于CTK推荐使用bwa,我们这里就按他们的习惯用bwa进行mapping。mapping的参数没有特殊的设置。
假设我们有一个基因list,并且已知一些RBP的结合位点,我们就可以确定这些是不是会富集某个RBP的结合位点。
Fold enrichment = (a / (a + b)) / (c / (c + d))
Input files
RBP binding sites / CLIP-seq peaks: human_ENCODE_eCLIP.bed
目标基因list对应的genomic intervals: up.bed
背景基因list对应的genomic intervals: up_background.bed
**Running steps: **
用bedtools对peaks的genomic intervals和基因list的genomic intervals取交集:
Calculate enrichment:
所以Fold enrichment = (698 / 1685) / (532 / 1685) = 1.312
如果使用singularity,则需要使用如下的指令进入singularity:
然后再运行程序,运行时需要注意将文件导入到自己所有的目录中,并且将所有运行脚本的/home/test/app/ctk/
改为/usr/local/bin/
。
根据上述代码复现CLIP-seq的分析过程,筛选出所有FDR < 0.05的CIMS,并提交对应的文件。
查阅资料回答,在CHIP-seq中是不是和CLIP-seq一样也存在PCR duplication导致假阳性peak的情况?如果有,homer和MACS2这类工具是怎么处理这一问题的?
https://zhanglab.c2b2.columbia.edu/index.php/CTK_Documentation
https://zhanglab.c2b2.columbia.edu/index.php/PARCLIP_data_analysis_using_CTK
Shah, A., Qian, Y., Weyn-Vanhentenryck, S. M., & Zhang, C. (2017). CLIP Tool Kit (CTK): a flexible and robust pipeline to analyze CLIP sequencing data. Bioinformatics (Oxford, England), 33(4), 566–567.
Hafner, M., Landthaler, M., Burger, L., Khorshid, M., Hausser, J., Berninger, P., Rothballer, A., Ascano, M., Jr, Jungkamp, A. C., Munschauer, M., Ulrich, A., Wardle, G. S., Dewell, S., Zavolan, M., & Tuschl, T. (2010). Transcriptome-wide identification of RNA-binding protein and microRNA target sites by PAR-CLIP. Cell, 141(1), 129–141.
Zhao, W., Zhang, S., Zhu, Y., Xi, X., Bao, P., Ma, Z., Kapral, T. H., Chen, S., Zagrovic, B., Yang, Y. T., & Lu, Z. J. (2021). POSTAR3: an updated platform for exploring post-transcriptional regulation coordinated by RNA-binding proteins. Nucleic acids research, gkab702. Advance online publication.
Format | Description | Notes |
---|---|---|
Format | Description | Notes |
---|---|---|
target region | background region | |
---|---|---|