# 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](https://zhanglab.c2b2.columbia.edu/index.php/CTK_Documentation)是[Chaolin Zhang lab](https://zhanglab.c2b2.columbia.edu/index.php/Main_Page)开发的一个处理CLIP-seq数据的工具,支持多种CLIP-seq实验的数据分析。在本章中，我们以PAR-CLIP为例介绍如何使用CTK处理CLIP的数据。

![Charles Danan et al., 2016](https://4115668567-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-LPVsf5VZbQ7h14X29qW%2Fuploads%2Fgit-blob-ad59ea44c5193b364243686918163a16a5ebf5a3%2FPAR-CLIP.png?alt=media)

## 2) Workflow

* 本章主要参考CTK文档中的[PARCLIP data analysis using CTK](https://zhanglab.c2b2.columbia.edu/index.php/PARCLIP_data_analysis_using_CTK)一节。有兴趣的同学可以直接阅读文档。
* CLIP-seq在数据分析上和我们[前面](https://book.ncrnalab.org/teaching/part-iii.-ngs-data-analyses/3.chip-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，只保留其中的一个。

![](https://4115668567-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2F-LPVsf5VZbQ7h14X29qW%2Fuploads%2Fgit-blob-2466730dce8cb9dddfbefeda4f87401108d9a3ef%2FPAR-CLIP-pipeline.png?alt=media)

## 3) Data structure

### 3a) getting software and data

**方法1: 使用docker**

1. install software (already available in Docker，docker images的下载链接如[附表](https://lulab2.gitbook.io/teaching/appendix/appendix-iv.-teaching#teaching-docker)所示): CTK, bwa, fastx\_toolkit, samtools
2. Get data (already available in Docker)，我们使用CTK tutorial中的两个PAR-CLIP样本：

   [SRR048967](https://www.ebi.ac.uk/ena/browser/view/SRR048967), [SRR048968](https://www.ebi.ac.uk/ena/browser/view/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](https://lulab2.gitbook.io/teaching/appendix/appendix-iv.-teaching#teaching-docker)

**方法2: 直接下载软件**

1. Install software:

   CTK: [download link](https://zhanglab.c2b2.columbia.edu/index.php/CTK)

   bwa: [download link](http://bio-bwa.sourceforge.net/)

   fastx\_toolkit: [download link](http://hannonlab.cshl.edu/fastx_toolkit/download.html)

   samtools: [download link](http://samtools.sourceforge.net/)
2. Get data: [download link](https://lulab2.gitbook.io/teaching/appendix/appendix-iv.-teaching#teaching-docker)

### 3b) input

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

### 3c) output

| Format                        | Description                        | Notes |
| ----------------------------- | ---------------------------------- | ----- |
| `*.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
>
> ```bash
> 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
> ```

首先进入到容器：

```bash
docker exec -it bioinfo_tsinghua bash
```

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

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

### 4a) Parsing alignment & Collapsing PCR duplicate

CTK先用`parseAlignment.pl`脚本从sam文件中提取出reads的位置信息和突变信息，方便后续分析。

{% hint style="info" %}
在CTK流程中，parseAlignment.pl这里输入的是gzip压缩后的sam文件，主要是因为sam文件比较占空间。在实际中我们更常见的做法是使用bam格式存储。
{% endhint %}

```bash
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`文件。

{% hint style="info" %}
在CLIP-seq和CHIP-seq中，"tag"是"read"的同义词。
{% endhint %}

接下来我们使用`tag2collapse.pl`，在起始位置相同的reads中只保留一个用于后续分析:

```bash
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对应的突变信息:

```bash
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的结合位点:

```bash
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:

```bash
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分析:

```bash
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的位点进行筛选:

```bash
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

```bash
# 解压缩

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的参数没有特殊的设置。

```bash
# 建立基因组索引（使用我们提供的基因组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                 |

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

**Input files**

1. RBP binding sites / CLIP-seq peaks: [human\_ENCODE\_eCLIP.bed](https://lulab2.gitbook.io/teaching/appendix/appendix-iv.-teaching#teaching-docker)

   ```
   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
   ```
2. 目标基因list对应的genomic intervals: [up.bed](https://github.com/lulab/teaching_book/blob/master/part-iii.-ngs-data-analyses/5.network/%E6%B8%85%E5%8D%8E%E4%BA%91%E9%93%BE%E6%8E%A5/README.md)

```
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. 背景基因list对应的genomic intervals: [up\_background.bed](https://github.com/lulab/teaching_book/blob/master/part-iii.-ngs-data-analyses/5.network/%E6%B8%85%E5%8D%8E%E4%BA%91%E9%93%BE%E6%8E%A5/README.md)

```
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取交集:

```bash
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:

```bash
# (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:

```bash
# 拷贝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

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.