2.RNA Regulation - RiboShape

Coupling structure dynamics with translation regulation

0) 背景介绍

RNA的转录后调控对于基因表达来说是至关重要的，单链RNA由于其不稳定性，总是倾向于形成有稳定的碱基互补配对的结构,其中最典型典型的就是茎环结构。不同的位置上的茎环结构的形成会直接影响到表达量，RNA分子的舒展程度可以很明显地影响核糖体沿RNA移动的过程。因此，RNA二级结构在不同条件下的转变最直接影响到的就是翻译过程。在以往的研究中发现，这种影响往往是介导转录调控(可变剪接AS、多聚腺苷酸加尾APA、嵌合体RNA等)进行的。

相较于使用RNA酶选择性切割单链或双链核苷酸的体外酶促RNA结构分析由于RNA酶体积过大且受限于体外反应不能准确的反映活细胞内的RNA真实的折叠状态，基于化学修饰的RNA结构探测方法由于相对体积更小能够更精细更敏锐地体内检测特定的RNA结构。通过引物延伸分析的选择性2’羟基酰化(SHAPE)与硫酸二甲酯化(DMS)是RNA修饰的两种主要的化学加合物。DMS方法的一个显著局限性就是缺失了将近一半的转录组RNA结构信息。因为DMS仅探测As(腺苷酸)和Cs(胞嘧啶)的结构信息。而缺少Us(尿嘧啶)和Gs(鸟苷酸)的配对状态信息。SHAPE方法由于其高分辨率与普适性，在近期研究中被广泛地用来探测mRNA二级结构调控元件以及其调控功能。

该大作业依托于SHAPE-MaP测序数据及Ribo-seq测序数据，测定在外界刺激前后，拟南芥可能产生的各种转录后调控，结构变化以及最终导致的翻译效率的变化。通过各种已有的工具和统计学分析方法，我们希望将两种条件下的差异剪接，差异表达，结构变化等信息与差异翻译联系起来，希望从中找寻到具有生物学意义的结果，例如建立与外界刺激密切相关的RNA二级结构动态模型，探究结构展开动力学与翻译效率之间的关系。

1) 总体流程图

2) 报告要求

报告要求：提交一份完整的工作报告，中英文不限(鼓励英文，可以参考一些发表的文献，如Pervasive Regulatory Functions of mRNA Structure Revealed by High-Resolution SHAPE Probing, 同时提交源代码。请读者使用我们提供的数据，完成下述分析过程，其中汇报关键结果需明确回答。

Part I. Prepare Data Matrix（选做）

完成RNA-seq，Ribo-seq，SHAPE-seq的数据预处理部分。

1 RNA-seq analysis：完成样本的Reads Processing、Remove RNA and Mapping工作，得到Mapped reads(bam)并绘制质量控制相关图，计算RNA-seq reads count matrix。
2 Ribo-seq analysis：完成样本的Reads Processing、Remove RNA and Mapping工作，得到Mapped reads(bam)并绘制质量控制相关图，计算Ribo-seq reads count matrix。
3 SHAPE-seq analysis：利用Shapemapper 进行数据预处理并计算每个转录本的reactivity。

汇报关键结果：
给定样本中AT1G09530.3的平均reactivity值为多少。

Part II. Data analysis

完成RNA-seq，Ribo-seq和SHAPE-seq的数据分析，得到外界刺激前后差异表达、差异翻译、RNA结构发生改变基因。

1 RNA-seq analysis：完成differential expression分析和differential splicing分析(选做)
- 1.1) 利用DESeq2和edgeR进行差异表达分析
- 1.2) 对差异表达基因进行通路注释并作图
- 1.3) 利用rMATS计算splicing events
- 1.4) 对差异剪接基因进行通路注释并作图
  汇报关键结果：
  外界刺激前后发生差异表达(DESseq)的基因总共有多少个，其中多少个基因上调，多少个基因下调；请给出外界刺激前后表达量上调和下调变化程度最大的基因id,其log2FoldChange值为多少。
  绘制差异表达基因的GO/KEGG通路分析。
2 Ribo-seq analysis：完成差异翻译效率分析。
- 2.1) 利用Xtail基于Ribo-seq的count matrix 和RNA-seq的count matrix计算差异翻译效率。为了方便对比RNAseq只统计CDS区域的reads。
  汇报关键结果：
  外界刺激前后发生差异翻译的基因总共有多少，其中多少基因上调，多少基因下调；请给出外界刺激前后差异翻译变化程度最大(上调和下调)的基因id,其差异翻译log2FC_TE_final值为多少。
  绘制差异翻译基因的GO/KEGG通路分析图
3 SHAPE-seq analysis：完成结构改变区域计算。
- 3.1）计算hit level，根据hit level分布确定阈值
- 3.2) 计算归一化因子，对Reactivity进行预处理和归一化。
- 3.3）计算滑动窗口Gini index，利用滑动窗口判断每个转录本外界刺激前后的结构改变区域，对连续的滑动窗口进行合并，汇总结构。
  汇报关键结果：
  绘制转录本随hit level变化趋势图，思考应该选择何值作为hit level阈值，为什么.
  用于归一化的转录本为哪几个。
  给定delta Gini index绝对值大于0.1为结构改变区域，外界刺激后有多少转录本发生了结构改变，多少转录本外界刺激后结构变多，多少转录本外界刺激后结构变少。

Part III. Data Integration

（1）转录本丰度变化和翻译效率(TE)变化程度之间关系，绘制Log2FoldChange(TE)和Log2FoldChange(RNA-seq)的关系图。
汇报关键结果：
- 绘制Log2FoldChange(TE)和Log2FoldChange(RNA-seq)的关系图，计算相关系数。
- 思考说明转录本丰度变化和TE负相关的原因。(选做)
（2）转录本结构改变程度和翻译效率（TE）变化程度之间关系。
- 2.1) 通过假设检验计算结构改变程度是否影响翻译效率TE变化程度；
- 2.2) 结构改变程度和翻译效率TE变化程度关系图。
  汇报关键结果：
  绘制结构改变程度和翻译效率TE变化程度关系图。
（3）翻译效率发生变化的基因的motif分析。
汇报关键结果：
- 翻译效率（TE）发生变化的基因的3‘UTR和5'UTR分别富集到的motif是什么？有什么不同？
- 寻找文献支持，尝试说明这些motif如何影响TE，可能的生物过程。
- 利用 MEME分析的是 sequence motif，我们能否利用一些预测structure motif 的软件和方法探索一些影响和调控翻译的structure motif？(选做，可以参考该教程：4.2 Structure Motif)

3）数据介绍

SHAPE-MaP测序数据及Ribo-seq测序数据。实验设计如下

实验使用Columbia(Col-0)作为野生型拟南芥植株。分别设置接受外界刺激的实验组和不接受外界刺激的对照组。随后对实验组和对照组植株进行SHAPE-MaP实验以及不进行化学修饰的对照实验以及Ribo-seq实验。SHAPE-MaP 使用NAI作为化学修饰物，而对照实验加入等量的DMSO溶液处理。每组实验均进行三次生物重复。实验数据如下：

表1 样本名称中，C代表Col-0植株;D代表SHAPE-MaP实验中使用DMSI对照处理，N代表NAI处理；0代表接受外界刺激，1代表未接受外界刺激。
表2 样本名称中，C代表Col-0植株;R代表Ribo-seq;0代表接受外界刺激，1代表未接受外界刺激。
由于SHAPE-MaP实验中作为对照组的，未加入小分子修饰的组别，与正常进行的RNA-seq得到的测序数据没有明显差异，因此用来做转录后调控分析的RNA-seq数据即为SHAPE-MaP实验中的对照组数据CD。

4）作业数据及软件存放地址

4.0) 数据集下载方式

本次大作业所用到的SHAPE-MaP数据及Ribo-seq数据的获取方式如下：

方法 1，有T集群账号的同学可以直接调用T集群的/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4目录下文件。
~~方法 2，远程下载这些数据： “Download Links of Assignments" @清华云盘。~~

4.1）Part I.Prepare Data Matrix相关文件

Part I. Prepare Data Matrix 包括RNA-seq，Ribo-seq，SHAPE-seq 的数据预处理三部分。这里为了节省计算时间，每个部分仅给出示例样本进行计算。

RNA-seq：选用WT型未受外界刺激不加NAI的SHAPE-seq的三个重复样本(CD1_1，CD1_2,CD1_3)的fastq文件进行数据预处理
Ribo-seq：选用WT型未受外界刺激的Riboseq的三个重复样本(CR1_1,CR1_2,CR1_3)的fastq文件进行数据预处理
SHAPE-seq：选用WT型未受外界刺激加NAI的CN1_1以及未加NAI的CD1_1的fastq文件用于shapemapper的计算。target fasta文件选用(AT4G21960.1,AT1G20620.4,AT1G29930.1)这三个转录本。

data

path

RNA-seq raw data

/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/SHAPE-MaP/WT/control/nouvb/CD1_*.clean.[1/2].fastq.gz

Ribo-seq raw data

/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/Ribo-seq/CR*.fq.gz

SHAPE-seq raw data

/data/TA_QUIZ_RNA_regulation/data/shapemapper_test/control/CD1_1.clean.part1.R[1/2].fastq.gz; /data/TA_QUIZ_RNA_regulation/data/shapemapper_test/modified/CN1_1.clean.part[1/2].R1.fastq.gz

SHAPE-seq target file

/data/TA_QUIZ_RNA_regulation/data/shapemapper_test/Arabidopsis_thaliana.TAIR10.34.transcripts_new_2.fa

dna.toplevel.fa

/data/TA_QUIZ_RNA_regulation/data/ATH/DNA/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa

gtf

/data/TA_QUIZ_RNA_regulation/data/ATH/GTF/Arabidopsis_thaliana.TAIR10.34.gtf

STAR genome index

/data/TA_QUIZ_RNA_regulation/data/ATH/index/STAR/genome/

rRNA index

/data/TA_QUIZ_RNA_regulation/data/ATH/index/bowtie1/rRNA/

RNA-seq raw data：CD1_1,CD1_2,CD1_3的原始fastq文件
Ribo-seq raw data：CR1_1,CR1_2,CR1_3的原始fastq文件
SHAPE-seq raw data：截取部分的CD1_1与CN1_1的原始fastq文件
SHAPE-seq target file：shapemapper target file
dna.toplevel.fa：拟南芥参考基因组序列文件
gtf：拟南芥参考基因组注释文件
STAR genome index：构建好的STAR genome index
rRNA index：rRNA的bowtie index

4.2）Part II. Data analysis 相关文件

Part II. Data analysis 包括RNA-seq，Ribo-seq，SHAPE-seq 的数据分三部分。我们给出每个类型数据预处理后的6个样本的完整结构，后续分析利用此处给出数据。

RNA-seq的Bam文件：包含6个样本
RNA-seq表达矩阵：分为exon区域计数和CDS区域计数，包含6个样本
Ribo-seq的Bam文件：包含6个样本
Ribo-seq表达矩阵：CDS区域计数，包含6个样本
SHAPE-seq：转录本信息整合文件，包含WT接受外界刺激和WT未接受外界刺激的两个文件。

data

path

RNA-seq Mapped reads

/data/TA_QUIZ_RNA_regulation/result/PartI.RNA-seq_analysis/differential_expression/3.mapping_expression/control/

RNA-seq expression count matrix(exon)

/data/TA_QUIZ_RNA_regulation/result/PartI.RNA-seq_analysis/differential_expression/exp_count_matrix/count_all.txt

RNA-seq expression count matrix(CDS)

/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/6.Differential_translation_efficiency/RNA-seq-CDS/count_CDS.txt

Ribo-seq Mapped reads

/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/3.mapping/

Ribo-seq read count

/data/TA_QUIZ_RNA_regulation/result/PartII.Ribo-seq_analysis/6.Differential_translation_efficiency/Ribo-seq/

SHAPE-seq merge files(WT未受外界刺激)

/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/final.modified_umodified/col_nouv/final.modified_unmodified

SHAPE-seq merge files(WT接受外界刺激)

/data/TA_QUIZ_RNA_regulation/data/riboshape_liulab_batch4/final.modified_umodified/col_uv/final.modified_unmodified

RNA-seq Mapped reads：RNA-seq所有样本预处理后的bam文件
RNA-seq expression count matrix(exon)：RNA-seq所有样本的表达矩阵数据(以exon计算),每一行为一个样本，每一列为一个基因。
RNA-seq expression count matrix(CDS)：RNA-seq所有样本的表达矩阵数据(以CDS计算),每一行为一个样本，每一列为一个基因。
Ribo-seq Mapped reads : Ribo-seq 所有样本预处理后的bam文件
Ribo-seq read count :Ribo-seq 所有样本的read count文件所在路径
SHAPE-seq merge files :我们把同一个实验条件三个replica合并后计算reactivity，得到final.modified_unmodified文件。文件每一行是一个转录本，每列内容为transcript_id,Nucleotide,Modified_mutations,Modified_effective_depth,Untreated_mutations,Untreated_effective_depth。