Comment on page
5.3. CLIP-seq (RNA-Protein Interaction)
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的数据。

Charles Danan et al., 2016
- 在前面的章节中我们已经对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
- 1.install software (already available in Docker,docker images的下载链接如附表所示): CTK, bwa, fastx_toolkit, samtools
- 2.Get data (already available in Docker),我们使用CTK tutorial中的两个PAR-CLIP样本:此处我们提供的是经过mapping的sam文件
SRR048967.sam.gz
和SRR048967.sam.gz
,位于/home/test/files/input/
目录下。为了提高运行速度,我们只选取了原始测序数据的前1,000,000条数据,如果感兴趣的话,也可以通过点击对应链接下载对应的fastq文件从头开始,建立mapping index所需的genome fasta文件可以从清华云下载:download link
方法2: 直接下载软件
- 1.Install software:
- 2.
Format | Description | Notes |
---|---|---|
*sam.gz | mapping之后得到的每一段序列的mapping情况 | - |
Format | Description | Notes |
---|---|---|
*.tag.uniq.peak.sig.bed | 经过peak calling过程之后被认为是显著富集的RBP结合区域 | - |
*.tag.uniq.t2c.CIMS.s13.txt | 经过CIMS分析之后被认为是显著的RBP互作位点 | - |
(这段加在6. RNA regulation里面)
(6) 6.10 CLIP-seq: bioinfo_clip_seq.tar.gzdocker load -i bioinfo_clip_seq.tar.gzdocker run --name=bioinfo_tsinghua -dt -h bioinfo_docker --restart unless-stopped -v ~/Desktop/lulab/gitbook/new_docker/share:/home/test/share zwhbio2017/clip-seq:latestdocker exec -it bioinfo_tsinghua bash
首先进入到容器:
docker exec -it bioinfo_tsinghua bash
以下步骤都在
/home/test/CLIP-seq
目录 下进行:cd /home/test/
mkdir CLIP-seq
cd CLIP-seq
CTK先用
parseAlignment.pl
脚本从sam文件中提取出reads的位置信息和突变信息,方便后续分析。在CTK流程中,parseAlignment.pl这里输入的是gzip压缩后的sam文件,主要是因为sam文件比较占空间。在实际中我们更常见的做法是使用bam格式存储。
for i in SRR048967 SRR048968
do
echo $i
perl /home/test/app/ctk/parseAlignment.pl -v --map-qual 1 --min-len 18 --mutation-file ${i}.mutations.txt /home/test/file/input/${i}.sam.gz ${i}.tag.bed
done
-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中只保留一个用于后续分析:for i in SRR048967 SRR048968
do
echo $i
perl /home/test/app/ctk/tag2collapse.pl -v -big -weight --weight-in-name --keep-max-score --keep-tag-name ${i}.tag.bed ${i}.tag.uniq.bed
done
其中参数的意义如下:
-v
:同上-big
:输入文件较大(采用特定算法)-weight
:在合并时考虑每一个序列的权重--weight-in-name
:序列权重在bed文件的name这一列中--keep-max-score
:在合并时保留权重最高的序列(默认保留长度最长的序列)--keep-tag-name
:保留bed文件中每一个序列原有的name
合并tag.bed中的重复序列之后,对于保留下来的reads,我们用
joinWrapper.py
这个脚本从mutations.txt
选出每一个reads对应的突变信息:for i in SRR048967 SRR048968
do
echo $i
python2 /home/test/app/ctk/joinWrapper.py ${i}.mutations.txt ${i}.tag.uniq.bed 4 4 N ${i}.tag.uniq.mutations.txt
done
joinWrapper.py
实际上是用python对gnu的join
命令进行了一个封装。在mutations.txt和uniq.bed文件中,第4列都对应着read id。 其中4 4代表使用两个文件第四列的read id进行连接,而N表示只会输出第4列相符的行。接下来我们通过peak calling推断RBP的结合位点:
for i in SRR048967 SRR048968
do
echo $i
perl /home/test/app/ctk/tag2peak.pl -big -ss -v --valley-seeking -p 0.05 --valley-depth 0.9 --multi-test --dbkey hg38 ${i}.tag.uniq.bed ${i}.tag.uniq.peak.sig.bed --out-boundary ${i}.tag.uniq.peak.sig.boundary.bed --out-half-PH ${i}.tag.uniq.peak.sig.halfPH.bed
done
这个命令的输入是
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:
for i in SRR048967 SRR048968
do
echo $i
perl /home/test/app/ctk/getMutationType.pl -t t2c --summary ${i}.tag.uniq.mutations.stat.txt ${i}.tag.uniq.mutations.txt ${i}.tag.uniq.t2c.bed
done
然后进行CIMS分析:
for i in SRR048967 SRR048968
do
echo $i
perl /home/test/app/ctk/CIMS.pl -big -n 10 -p -outp ${i}.tag.uniq.t2c.pos.stat.txt -v ${i}.tag.uniq.bed ${i}.tag.uniq.t2c.bed ${i}.tag.uniq.t2c.CIMS.txt
done
这个脚本的输入是去过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的位点进行筛选:for i in SRR048967 SRR048968
do
echo $i
awk '{if($9<=0.05) {print $0}}' ${i}.tag.uniq.t2c.CIMS.txt | sort -k 9,9n -k 8,8nr -k 7,7n > ${i}.tag.uniq.t2c.CIMS.s13.txt
done
前面的介绍中,我们是直接从mapping的结果出发的,以下介绍了mapping以及mapping之前的数据预处理步骤供大家参考。
# 解压缩
gunzip SRR048967.fastq.gz
gunzip SRR048968.fastq.gz
# 1. 去除接头序列
fastx_clipper \
-a AGATCGGAAGAG \ # 指定接头序列。接头序列的不同因实验不同而异,需要结合实际情况进行设置
-l 13 \ # 如果去除接头序列之后长度不足13nt,则舍弃这些read
-Q 33 \ # fastq文件使用Phred33编码描述测序质量
-v \ # 输出更多信息
-i SRR048967.fastq \ # 输入文件
-o SRR048967.clipped.fastq # 输出文件
# 2. 过滤测序质量过低的序列
fastq_quality_filter \
-q 20 \ # 控制测序质量的cutoff
-p 80 \ # 要求测序质量超过cutoff的nucleotide必须超过80% (即只保留超过80% nucleotide有高于20质量的read)
-Q 33 \ # fastq文件使用Phred33编码描述测序质量
-v \ # 输出更多信息
-i SRR048967.clipped.fastq \ # 输入文件
-o SRR048967.clipped.filtered.fastq # 输出文件
# 3. 合并PCR重复
fastx_collapser \
-Q 33 \ # fastq文件使用Phred33编码描述测序质量
-v \ # 输出更多信息
-i SRR048967.clipped.filtered.fastq \ # 输入文件
-o SRR048967.clipped.filtered.collapsed.fa # 输出文件 (注意,此处因为合并重复序列之后无法合并测序质量数值,因此文件格式转变为fasta,不再包含测序质量信息)
CLIP-seq测的是RNA序列,理论上用STAR/hisat2这类spliced aligner会更合适,但是由于CTK推荐使用bwa,我们这里就按他们的习惯用bwa进行mapping。mapping的参数没有特殊的设置。
# 建立基因组索引(使用我们提供的基因组fasta文件)
bwa index -a bwtsw hg38.chr.fa
# 比对序列
bwa aln -t 4 -n 0.06 -q 20 hg38.chr.fa SRR048967.clipped.filtered.collapsed.fa > SRR048967.sai
bwa samse hg38.chr.fa SRR048967.sai SRR048967.clipped.filtered.collapsed.fa | gzip -c > SRR048967.sam.gz
假设我们有一个基因list,并且已知一些RBP的结合位点,我们就可以确定这些是不是会富集某个RBP的结合位点。
Text | target region | background region |
---|---|---|
bound by RBP | a | c |
unbound by RBP | b | d |
Fold enrichment = (a / (a + b)) / (c / (c + d))
Input files
- 1.chr1 113808 113960 human_RBP_eCLIP_ENCODE_1270957 0 - HNRNPL eCLIP K562chr1 172556 172588 human_RBP_eCLIP_ENCODE_487541 0 - DDX24 eCLIP K562chr1 172556 172590 human_RBP_eCLIP_ENCODE_466544 0 - DDX24 eCLIP K562
- 2.
chr1 1167104 1167134 ENSG00000207730.3__1167104__1167134__1 . +
chr1 1167119 1167149 ENSG00000207730.3__1167119__1167149__2 . +
chr1 1167134 1167164 ENSG00000207730.3__1167134__1167164__3 . +
- 1.
chr1 1185074 1185104 ENSG00000162571.13__1185074__1185104__95 . +
chr1 1328816 1328846 ENSG00000224051.6__1328816__1328846__141 . +
chr1 1401351 1401381 ENSG00000224870.6__1401351__1401381__104 . +
**Running steps: **
- 1.用bedtools对peaks的genomic intervals和基因list的genomic intervals取交集:
mkdir /home/test/rbp
cd /home/test/rbp
# overlap between RBP binding sites and target regions
bedtools intersect -a ../file/input/human_ENCODE_eCLIP.bed -b ../file/input/up.bed -wa -wb -sorted -s > up_eCLIP.txt
# overlap between RBP binding sites and background regions
bedtools intersect -a ../file/input/human_ENCODE_eCLIP.bed -b ../file/input/up_background.bed -wa -wb -sorted -s > up_background_eCLIP.txt
- 1.Calculate enrichment:
# (a + b)
cut -f 4 ../file/input/up.bed | sort -u | wc -l
1685
# a
cut -f 13 up_eCLIP.txt | sort -u | wc -l
698
# (c + d)
cut -f 4 ../file/input/up_background.bed | sort -u | wc -l
1685
# c
cut -f 13 up_background_eCLIP.txt | sort -u | wc -l
532
所以Fold enrichment = (698 / 1685) / (532 / 1685) = 1.312
如果使用singularity,则需要使用如下的指令进入singularity:
# 拷贝perl library至自己的工作目录下
mkdir perllib
cp -r /data/zhaoweihao/app/perl/lib/perl5/x86_64-linux/* perllib/
cp -r /data/zhaoweihao/app/perl/lib/czplib/* perllib/
source /WORK/Samples/singularity.sh
singularity run /data/images/bioinfo_clip_seq.simg
export PERL5LIB=`pwd`/perllib
然后再运行程序,运行时需要注意将文件导入到自己所有的目录中,并且将所有运行脚本的
/home/test/app/ctk/
改为/usr/local/bin/
。- 根据上述代码复现CLIP-seq的分析过程,筛选出所有FDR < 0.05的CIMS,并提交对应的文件。
- 查阅资料回答,在CHIP-seq中是不是和CLIP-seq一样也存在PCR duplication导致假阳性peak的情况?如果有,homer和MACS2这类工具是怎么处理这一问题的?
- 1.https://zhanglab.c2b2.columbia.edu/index.php/CTK_Documentation
- 2.https://zhanglab.c2b2.columbia.edu/index.php/PARCLIP_data_analysis_using_CTK
- 3.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.
- 4.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.
- 5.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.
Last modified 1yr ago