5.3. CLIP-seq (RNA-Protein Interaction)

1) Background

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的数据。

2) Workflow

本章主要参考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，只保留其中的一个。

3) Data structure

3a) getting software and data

方法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

3b) input

Format Description Notes

Format	Description	Notes
`*sam.gz`	mapping之后得到的每一段序列的mapping情况	-

*sam.gz

mapping之后得到的每一段序列的mapping情况

3c) output

Format Description Notes

Format	Description	Notes
`*.tag.uniq.peak.sig.bed`	经过peak calling过程之后被认为是显著富集的RBP结合区域	-
`*.tag.uniq.t2c.CIMS.s13.txt`	经过CIMS分析之后被认为是显著的RBP互作位点	-

*.tag.uniq.peak.sig.bed

经过peak calling过程之后被认为是显著富集的RBP结合区域

*.tag.uniq.t2c.CIMS.s13.txt

经过CIMS分析之后被认为是显著的RBP互作位点

4) Running steps

(这段加在6. RNA regulation里面)

(6) 6.10 CLIP-seq: bioinfo_clip_seq.tar.gz

docker load -i bioinfo_clip_seq.tar.gz

docker run --name=bioinfo_tsinghua -dt -h bioinfo_docker --restart unless-stopped -v ~/Desktop/lulab/gitbook/new_docker/share:/home/test/share zwhbio2017/clip-seq:latest

docker exec -it bioinfo_tsinghua bash

首先进入到容器：

docker exec -it bioinfo_tsinghua bash

以下步骤都在/home/test/CLIP-seq目录下进行：

cd /home/test/
mkdir CLIP-seq
cd CLIP-seq

4a) Parsing alignment & Collapsing PCR duplicate

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列相符的行。

4b) Peak calling

接下来我们通过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：输出峰一半高度处的边界。

4c) CIMS analysis

我们知道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

5) Tip/Utilities

前面的介绍中，我们是直接从mapping的结果出发的，以下介绍了mapping以及mapping之前的数据预处理步骤供大家参考。

5a) Pre-processing

# 解压缩

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，不再包含测序质量信息）

5b) Mapping

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

5c) Enrichment of binding sites

假设我们有一个基因list，并且已知一些RBP的结合位点，我们就可以确定这些是不是会富集某个RBP的结合位点。

	target region	background region
bound by RBP	a	c
unbound by RBP	b	d

target region

background region

bound by RBP

unbound by RBP

Fold enrichment = (a / (a + b)) / (c / (c + d))

Input files

RBP binding sites / CLIP-seq peaks: human_ENCODE_eCLIP.bed

chr1	113808	113960	human_RBP_eCLIP_ENCODE_1270957	0	-	HNRNPL	eCLIP	K562
chr1	172556	172588	human_RBP_eCLIP_ENCODE_487541	0	-	DDX24	eCLIP	K562
chr1	172556	172590	human_RBP_eCLIP_ENCODE_466544	0	-	DDX24	eCLIP	K562

目标基因list对应的genomic intervals: up.bed

chr1	1167104	1167134	ENSG00000207730.3__1167104__1167134__1	.	+
chr1	1167119	1167149	ENSG00000207730.3__1167119__1167149__2	.	+
chr1	1167134	1167164	ENSG00000207730.3__1167134__1167164__3	.	+

背景基因list对应的genomic intervals: up_background.bed

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: **

用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

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

5d) Tips when using singularity

如果使用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/。

6) Homework

根据上述代码复现CLIP-seq的分析过程，筛选出所有FDR < 0.05的CIMS，并提交对应的文件。
查阅资料回答，在CHIP-seq中是不是和CLIP-seq一样也存在PCR duplication导致假阳性peak的情况？如果有，homer和MACS2这类工具是怎么处理这一问题的？

7) References

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.

Previous5.2.miRNA Targets Next6.RNA Regulation - I

Last updated 2 years ago