RNA sequencing data exploration and quality assessment.Rmd

---
title: 'RNA sequencing data exploration and quality assessment'
author: 'Yifei Wan'
date: 'Apirl/5/2018'
output: html_document
---

## 1. 测序数据读入与整理

```{r}
#rawdata <- read.csv2('LncRNA_all.txt', sep = '\t', row.names = 1) 
library('pasilla')
pascount <- system.file('extdata', 'pasilla_gene_counts.tsv', package='pasilla', mustWork=TRUE)
pasmeta <- system.file('extdata', 'pasilla_sample_annotation.csv', package='pasilla', mustWork=TRUE)
rawdata <- as.matrix(read.csv2(pascount,sep='\t',row.names='gene_id'))
coldata <- read.csv2(pasmeta, row.names=1)
coldata <- coldata[,c('condition','type')]
```
  使用`pasilla`包的内置测序数据作为本文的数据源。`rawdata`内为保存测序读数的数字矩阵。`coldata`中保存实验条件。在函数read.csv2中指定制表符tab为分隔符并指明行名。

接下来让我查看一下数据的基本结构：

```{r}
head(rawdata)
```
  可以清楚的从控制台返回的结果看出这是一个以转录本名称为行（row），以样本为列（column）的表格。

```{r, warning = FALSE}
conditions <- c(rep('control', 3), rep('treatment', 3)) ## 以实验条件作为列名
cols <- paste(conditions, 1:3, sep = '_')               ## 为各组编号
colnames(rawdata) <- cols                               ## 将新列名赋予原始表格
colnames(rawdata)                                       ## 打印新列名到控制台
```
    
## 2. 数据的初步探索
  我们会利用`ggplot2`包对原始数据进行可视化的探索，以求对数据的分布的有一个初步认识。可视化探索的结果会影响我们对后续策略的选择。
  载入`ggplot2`:
```{r}
library(ggplot2)
```

### 2.1 直方图显示单个样本分布
  直方图是最简单的作图方式之一，可以非常直观的展示数据的分布。下面我们以`control_1`为例绘制直方图：
```{r}
ggplot(data = rawdata, aes(x = control_1)) + 
  geom_histogram(fill = 'purple3', binwidth = 5000)
```
  
  直方图呈现出一种奇怪的结果：几乎所有数据都挤在0附近的一根bar上面，而且肉眼可见的bar只有这一根。这就说明绝大多数的的测序结果都是很小的数字。但是如果对水平坐标轴进行观察，我们可以发现坐标轴横跨了从0到150,000的巨大区间。这暗示我们有稀少种类的RNA有远高于其他RNA的表达水平。这也符合我们大部分RNA表达水平较低而小部分RNA表达水平极高的假设（此假设是进行表达差异分析的基本假设）。
  为了将种类少而表大量高的RNA也展示在直方图上，我们需要对数据整体进行转化（transformation），最常见的方式是取对数。但是因为数据中有0存在，所以为了避免出错需要对每一个数据点` + 1`。当然也可以加其他合理小量。
    完成转化之后重新绘图：

```{r}
logdata <- log2(rawdata + 1)                                        ## 对原始数据取对数
ggplot(data = logdata, aes(x = control_1)) +
  ylab(expression(log[2](count + 1))) +                             ## 使用新y-lab
  geom_histogram(colour = 'white', fill = 'purple', binwidth = 0.7) ## 添加白色包线
```
  
  新直方图清晰的揭示了测序结果的分布情况也证明了我们上文的推测：表达水平高的RNA种类很少，绝大多数RNA都处于较低的表达水平换言之数据点极其分散。结合前文`head`命令所显示的原始数据，我们可以知道测序结果会是一种方差极大的离散数据（都是整数），这与microarray的连续数据（有非整数）有本质差异，也决定了后续分析方法将会不同。并且巨大的方差（过离散：overdispersion）导致过去测序分析中常见的泊松分布逐渐被负二项分布取代（泊松分布的期望和方差均等于其参数lambda, 这导致方差大小受限与数据实际情况不符）。
  基于以上的结果，我们需要在后续的分析中对大量低表达数据进行过滤（filter）与移除，因为低表达数据无法与噪音相区隔会对分析过程造成误导。
  
### 2.2 箱形图探索组间分布差异
  直方图展示了单组数据的分布，接下来我们需要探索组间的差异。箱形图可以将各组的的总体分布情况展示在同一途中，对于初步评估方差差异或观察实验与分析中的各种处理效应都有助益。需要注意的是ggplot的箱形图只接受两列的数据，所以需要对数据结构进行重构（reshape）。这里使用`reshape`包的`melt`函数：
```{r}
library(reshape)
boxdata <- melt(data = logdata, variable_name = 'Count')
head(boxdata)
```
  可以看出数据被重构为两列，第一列是数据来源的组别，第二列是对应的拷贝数。接下来作出箱形图：
```{r}
groups <- substr(boxdata$Count, 1, 7)                               ## 标记数据组别来源的标签：treatment没用完整拼写
ggplot(data = boxdata, aes(x = Count, y = value, fill = groups)) +
  geom_boxplot() +
  ylab(expression(log[2](count + 1))) +
  scale_fill_manual(values = c('steelblue3', 'hotpink1'))
```
  
  箱形图中的方块囊括了从25%~75%的数据点，箱体中的黑线是平均值。显然各组的数据在75%之后都拖有长长的尾巴。这再次印证了上文基于直方图做出的结论：大量的低表达RNA需要被移除。同时我们可以看出同种处理条件下的箱体对的并不是很齐，比如`control_1`与`control_2`之间箱体有很明显的差异。这种差异可能是测序过程中出现的。如果原始数据有提供测序方法（如单端测序或双端测序）等额外信息就可以通过聚类等手段对误差出现原因进行估计。但是遗憾的是本数据并未提供此类数据。
### 2.3 密度图
  箱形图虽然给出了各组平均数的位置，但是毕竟只是算术平均数。如果数据存在多于一个众数，那么平均数就无法用来描述数据的集中趋势了。利用密度图我们可以更加详细观察各组数据的分布：
```{r}
boxdata <- data.frame(boxdata, condition = groups)                       ## 添加一列，将数据分为对照与处理两类
ggplot(data = boxdata, aes(x = value, colour = Count, fill = Count)) +   ## 以带序号的组名为准填充颜色和勾线
  ylim(c(0, 0.5)) +                                                      
  geom_density(alpha =0.1, size = 1) +                                   ## 参数alpha设置图像透明度，size设置勾线粗细
  facet_wrap(~ condition) +                                              ## 以刚才添加的分类为准分割图像为两框
  theme(legend.position = 'top') +                                       ## 设置图注的位置
  xlab(expression(log[2](count + 1)))                                    ## 设置x轴标注
```
  
  密度图清晰的显示出各个组都有两个峰值，但是较高的峰值出现在低表达区（接近于0）。所以我们可以期待在过滤低表达数据后数据将形成一个只有单一集中趋势的分布，这对后期的拟合与DE分析是有利的。另外，虽然在箱形图中统一处理条件下`control_1`与其他对照组有明显差异，但是从密度图来看这可能只是低表达数据引起的误导，几个对照组的数据分布可以较好的重叠在一起。
### 2.4 MA图
  MA-plot是一种展示样本间表达倍数差异的图像，因为数据都被取过对数所以被称为对数表达差异倍数（log2 fold change）。两样本间的对数表达差异倍数被称为M-value，而各个转录本表达平均值被称为A-value。其中M作为y值，A作为x值。以`control_2`与`treatment_1`为例，作图代码如下：
```{r}
M <- logdata$control_2 - logdata$treatment_1                        ## 计算对数表达差异倍数M-value
A <- (logdata$control_2 + logdata$treatment_1)/2                    ## 计算表达平均值
judge <- sapply(M, function(x , threshold = 1){                     ## 以±1为界根据表达上调或下调为表达差异分类为：up, down, none
  if(x > threshold){x = 'up'
  } else if (x < -threshold){x = 'down'
  } else {x = 'none'}
})

MA = data.frame(A, M)
ggplot(MA, aes(x = A, y = M)) +                                     
  geom_point(aes(color=judge), size = 1.5, alpha = 1/5) +           ## 利用judge为数据点分类
  geom_hline(aes(yintercept=0), color = 'blue3') +                  ## 添加水平线
  stat_smooth(se = FALSE, method = 'loess', color = 'red3') +       ## 添加趋势线（采用局部加权回归loess）
  scale_color_manual(values=c('forestgreen', 'grey', 'firebrick1')) ## 根据分类为数据点配色
```
  
  MA图中蓝色标明了上调的转录本，红色则是下调，在x轴附近的绿点是没有变化的转录本。这只是一个粗略结果，还不能说具有统计差异。同时平均表达水平5以下的转录本呈射线状，极其分散、整齐和对称，这暗示这些数据主要来自测序系统自身的噪音（泊松噪音），应该过滤掉。所以我们后文将采用5作为阀值初略过滤低表达数据。另外，呈箭头状的图像主体（称为Fan-shape)应以x轴为对称轴近似对称，因为基于假设大部分的转录本都没有显著变化。添加的红色趋势线显示了实际的轴线，与x轴有一定偏移。这再次证明数据归一化是需要的。
** Hint **:
严格的低表达过滤是基于严格的统计学方法在归一化处理后的数据上实施的，而预过滤只是为初步的探索数据形态而基于分析者的经验进行的。后者决不可取代前者，也不是必须的。但是在初步探索阶段预过滤对于提高运算和作图效率的帮助是很明显的。
  
### 2.5 变化倍率散点图
  用散点图来展示对数表达倍数是一种类似于MA图的方法，但是只需要计算M值不需要计算A值。我们将使用上图的数据直接绘制散点图，注意下图在添加趋势线时使用了与上图不同的函数，但是效果是相同的：
```{r}
ggplot(data = logdata, aes(x = treatment_1, y = control_2) ) + 
  geom_point(aes(color=judge), alpha=1/2, size=0.8) + 
  scale_color_manual(values=c('forestgreen', 'grey', 'firebrick1')) + 
  geom_smooth(method = 'loess',se = FALSE, colour = 'red3') +         ## 添加趋势线
  geom_abline(intercept = 0, slope = 1, colour = 'blue')              ## 添加对角线
ggtitle('Scatter plot for log2fc')                                    ## 添加图片标题
```
  
  散点图以两组参与比较的数据分别作为两轴，被减数control为y,减数treatment为x。此图的判读方法与MA图一致。趋势线同样显示数据需要归一化。
### 2.6 聚类
  机器学习算法中的聚类是判断测序质量的重要方法。简言之，聚类依靠计算样本间的相似度（在数据空间内的距离）为样本自动分类（无人为指导，所以属于无监督学习）。如果实验实施得当，我们可以预见接受相同处理条件的样本应当被聚为一类。而聚类的结果可以有多种表达方式，一种是直接使用dendrogram（系统树状图），实例请参考[R语言绘图-常用参数](http://www.bioengx.com/2018/04/04/r%e8%af%ad%e8%a8%80%e7%bb%98%e5%9b%be-%e5%b8%b8%e7%94%a8%e5%8f%82%e6%95%b0/)。
  而另一种方式就是大家喜闻乐见的热图（heat-map），利用色块的颜色深浅来直观显示样本间的相似度。同时热图是可以和树状图结合在一起的。接下来我们将使用`gplots`包中的`heatmap.2`函数进行绘图。吃函数是R自带的`stats`包内的`heatmap`函数的升级版本，可以绘制更加丰富的细节。上文已经反复提到我们需要过滤过多的低表达RNA，他们对聚类和后文的显著性分析都有不利影响。所以在绘图之前先做一个简单的预过滤（pre-filter）是有帮助的。但是后文为了进行更加精确的DE分析我们会基于原始数据用更加严格的方法加以过滤。
```{r}
preindex <- rowSums(logdata >= 5 ) >= 3 ## 找出至少三个样本对数表达值不小于5的转录本在数据中的行号
predata  <- logdata[preindex, ]         ## 提取通过了过滤的转录本
predatat <- t(as.matrix(predata))       ## 将数据转为绘图函数可接受的矩阵且转置。热图是以行为样本与测序数据结构不同，所以需转置 

library(gplots)                                               
heatmap.2(predatat, scale = 'row', key = T, keysize = 2, trace = 'none', dendrogram = 'row', cexRow=0.8, cexCol = 0.8, labCol = F, main = 'Cluster of samples')
##绘图：col设置颜色，scale设定以行为准聚类，key显示标尺, trace关闭基准线, dendrogram为样品添加树状图
```
  
  聚类结果和预想差别较大，这是在实际工作中常常遇到的情况。处理组和对照组无法被各自聚为一类。可能造成这一结果的原因很多，比如实验的处理条件并没有显著效果，也可能是样本文库准备出了问题还可能是测序噪音掩盖了真实情况。当无论如何，我们应当提高警觉。出现此种情况暗示我们最终的结论的是不可靠的。但是放弃还言之过早，我们还有其它工具可以进一步探究聚类失败的成因。
  我们可用如下方法进行分析：
+ 运用统计重新转化数据
  上文我们只是简单运用log转化数据，消除过大的尺度差异，然后基于图像选择了过滤的阈值。这种经验方法通常可以处理大部分的问题，但是当数据集很小时(n<30)就可以使用正则化对数转化（regularized-logarithm）。rlog是基于log的一种转化。当表达量较大时其结果与单纯log接近，否则会对结果施加惩罚将其压缩（shrink），相当于做了过滤。
  著名的DESeq2包提供了rlog函数：
```{r}

```