简单全面的 Keth-seq 数据分析教程

生信店小二 · 简书 · · 2021-04-17 15:35

前言

今天来跟大家分享一个比较少见的数据类型分析——Keth-seq。不知道大家第一眼看到这个类型是什么感受，反正我的感受就是还有这种数据类型，看来还是在下知识浅薄孤陋寡闻了。不过，不知道没有关系，咱可以学习啊！一顿搜索过后，还是能找到一些资料的，经过一番折腾数据分析流程终于走通了。下面来跟大家做一个分享。
RNA分子具有折叠成复杂结构的内在能力，这些结构对调节许多生物过程非常重要，包括转录、翻译、加工和降解。而对于RNA二级结构的分析可使用icSHAPE(Selective 2′ Hydroxyl Acylation analyzed by Primer Extension)来流程分析，该流程的原理是通过使用2-甲基烟酸咪唑在烟酸环(NAI-N3)上添加叠氮化物，来修饰未配对的RNA核苷酸序列。然后，可以通过生物素来进一步来富集被修饰后的片段以便后续测序分析。Kethoxal (1,1-dihydroxy-3-ethoxy-2-butanone) 在温和的条件下能够与单链RNA (ssRNA)中的鸟嘌呤反应，并诱导逆转录(RT)的停止。与icSHAPE方法类似，azido-kethoxal (N3-kethoxal, 1) 通过特异性地标记ssRNA链中鸟嘌呤链接处的N1和N2位置，然后结合深度测序的方法来探测RNA二级结构，即Keth-seq。由于Keth-seq方法构建的文库与icSHAPE类似，因此，同样可以使用icSHAPE流程的脚本来处理Keth-seq的测序数据。
上面一段文字简单描述了Keth-seq的原理以及该方法的分析目的。简单来说就是Keth-seq能够标记RNA链中没有配对的G碱基，然后通过检测这些G碱基的位置达到探测RNA二级结构的目的。废话就不多说了，下面来说一下具体的数据分析过程。

分析流程

下图是分析流程的示意图：

数据处理

1、Collapsing the reads
使用流程的 readCollapse.pl 脚本来完成这一步骤，默认参数即可。完全相同序列的读码被标记为PCR重复，并在后续分析前会被过滤掉。这里强调一下， 输入的fastq文件需要提前解压好，不能是gz压缩格式 。代码如下所示：

readCollapse.pl  -U sample.fastq -o sample_rmdup.fq -f sample_collapse.fa

2、去除接头序列
使用流程的 trimming.pl 脚本来删除序列中的 adapters 以及低质量碱基。代码如下所示：

trimming.pl -U sample_rmdup.fq -o sample_trimmed.fq -l 13 -t 0 -c phred33 -a adapter.fa -m 0

3、比对
首先，将所有的reads比对到核糖体RNA上，保留未必对上的reads以达到去除rRNA的目的，然后将保留的reads比对参考基因组上。比对使用的软件是 Bowtie 。代码如下所示：

bowtie --sam-nohead --quiet -p 5 -S --un sample_rmrrna.fq sample_trimmed.fq sample_rmrna.sam

bowtie --quiet -p 5 -S genome sample_rmrrna.fq sample.sam

4、计算RT信号
使用流程的 calcRT.pl 脚本来计算RT-stop信号。代码如下所示：

calcRT.pl -i sample.sam -o sample.rt -r sample.rpkm -c 5

5、合并RT信号(可选步骤)
如果有重复的情况下，可使用流程的 combinertreplates .pl 脚本来合并重复间的RT信号。代码如下所示：

combineRTreplicates.pl -i sample_rep1.rt:sample_rep2.rt -o sample_combine.rt

6、标准化RT信号
使用流程的 normalizeRTfile.pl 脚本来分别对处理组和对照组样品进行RT信号的标准化。代码如下所示：

normalizeRTfile.pl -i sample_combine.rt -o sample_normed.rt -m mean:vigintile2 -d 32 -l 32

7、计算 reactivity score
通过使用流程的 calcEnrich.pl 脚本比较处理样品(前景)和对照样品(背景)，然后计算每个转录本中每个核苷酸的 reactivity score 来作为RNA结构的评估得分。代码如下所示：

calcEnrich.pl -f sample_case_normed.rt -b sample_control_normed.rt -o sample_enrich.out -w factor5:scaling1 -x 0.25 -y 10

8、过滤 reactivity score
最后，为了获得高质量的 reactivity score，使用流程的 filterrich .pl 脚本对低质量的结果进行过滤。代码如下所示：

filterEnrich.pl -i sample_enrich.out -o sample_enrich_filter.out -t 200 -T 2 -s 5 -e 30

为了大家能有一个更为直观的感受，我这里展示一下最终得到的文件内容，如下所示：

ENSMUST00000117219.2    1545    2031.080251985  NULL    NULL    NULL    NULL    NULL    NULL    0.337   1.000   0.273   0.145   0.243   0.278   0.103   0.092   0.03
ENSMUST00000221295.2    556     187.220512465203        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000025271.17   1362    619.179886358961        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000133910.3    849     169.099803070766        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000119383.2    555     189.536742599081        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000153590.2    834     347.609558826383        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000031243.15   1410    107.555258591905        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000235468.2    593     84.462331623428 NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000088336.3    294     667.237103151603        NULL    NULL    NULL    NULL    NULL    NULL    0.084   0.240   0.230   0.444   0.340   0.070   0.290   0.12
ENSMUST00000172132.10   1607    131.447561137884        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000119822.2    798     291.101195382799        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000077915.10   364     259.473892898157        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000117557.8    2019    162.738982819819        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000074353.6    495     179.75851041931 NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000173844.8    382     276.731461610335        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000082402.1    1545    7076.29904084836        NULL    NULL    NULL    NULL    NULL    NULL    1.000   1.000   0.613   0.254   0.900   0.535   0.043   0.34
ENSMUST00000118499.2    750     213.193368235271        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000129380.2    845     1245.55315764876        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000117757.9    1273    134.028421025852        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000173253.2    979     180.173196018338        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000119790.2    793     172.560595436876        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000016463.4    1240    245.385465919021        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000099371.5    276     172.660742860779        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000136346.2    260     139.892599922161        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000106599.8    457     190.224308764774        NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL    NULL
ENSMUST00000082419.1    519     1088.91055899109        NULL    NULL    NULL    NULL    NULL    NULL    0.000   0.000   0.106   0.450   0.000   0.409   0.123   0.19

列说明：
第一列：转录本的Ensembl ID；
第二列：转录本的长度；
第三列：FPKM值；
第四-最后；每个碱基的 reactivity score。

其实前面的1-8步可以使用icSHAPE流程一步来完成，前面之所以展示分步式是想大家可以了解流程中都做了哪些数据处理，一键分析前需要在“icshape.conf”文件配置好流程需要的软件和参考基因等信息，一键式分析代码如下所示：

icSHAPE_pipeline.pl -i notreat1.fastq:notreat2.fastq -t kethrep1.fastq:kethrep2.fastq -o icshape_output -c icshape.conf

9、结果可视化
可利用 IGV 软件来展示Keth-seq的结果以显示RNA的结构信号。在此之前，需要得到可以用来展示的数据。处理数据的过程如下所示：

enrich2Bedgraph.pl -i sample.out -o sample.bdg -g gtf -a fasta

sort -k1,1 -k2,3n sample.bdg >sample_sorted.bdg

uniqueTrack.pl sample_sort.bdg sample_uniq.bdg

cut -f1-4 sample_uniq.bdg | grep -v NULL > sample_sim.bdg

bedGraphToBigWig sample_sim.bdg genome.chr.size sample_sim.bw

最后可用得到的bigwig文件在IGV中进行可视化展示，这里展示一个效果图，如下所示：

最后

使用icSHAPE流程的脚本来处理 Keth-seq 的数据还是相当快捷方便的。今天就分享到这里，后面附上了一些 Keth-seq 相关的参考文献，需要的朋友可以看一看。

参考文献

[1] X Weng, Gong J , Chen Y , et al. Keth-seq for transcriptome-wide RNA structure mapping[J]. Nature Chemical Biology, 2020, 16(5):1-4.
[2] Li P, Shi R, Zhang Q C. icSHAPE-pipe: A comprehensive toolkit for icSHAPE data analysis and evaluation[J]. Methods, 2019.
[3] Spitale, R. C. et al. Structural imprints in vivo decode RNA regulatory mechanisms. Nature 519, 486–490 (2015).
[4] R.C. Spitale, R.A. Flynn, Q.C. Zhang, P. Crisalli, B. Lee, J.W. Jung, H.Y. Kuchelmeister, P.J. Batista, E.A. Torre, E.T. Kool, H.Y. Chang, Structural imprints in vivo decode RNA regulatory mechanisms, Nature 519 (7544) (2015) 486–490.
[5] Xu, Z. & Culver, G.M. In Methods in Enzymology; Biophysical, Chemical, and Functional Probes of Rna Structure, Interactions and Folding, Pt A (ed. Herschalag, D.) Vol 468, 47–165 (Academic Press, 2009).
[6] Weng X , Gong J , Chen Y , et al. Keth-seq for transcriptome-wide RNA structure mapping[J]. Nature Chemical Biology, 2020, 16(5):1-4.

原文地址：访问原文地址
快照地址：访问文章快照

分享到微博