从表达量矩阵求各处理组的平均值

link171 2020-01-20

展开全文

上一篇文章我提到了从表达量矩阵画时序性（或者多处理）的单基因的折线图，

实际上我们的表达量矩阵里面每个处理是有2-3个重复的。所以用来画折线图的矩阵应该是合并了重复，求得平均值以后的表达量矩阵。这一部分我们来尝试从最原始的表达量矩阵合并重复、求平均值。

首先我们来构建测试数据

set.seed(19960203)
library(magrittr) # 只是为了单纯调用%>%这个管道操作而已
expr_df <- data.frame(Control_R1 = rnorm(5,mean = 10),
Control_R2 = rnorm(5,mean = 10),
Control_R3 = rnorm(5,mean = 10),
Treat_R1 = rnorm(5,mean = 20),
Treat_R2 = rnorm(5,mean = 20),
Other_R1 = rnorm(5,mean = 30),
Other_R2 = rnorm(5,mean = 30),
Other_R3 = rnorm(5,mean = 30),
Other_R4 = rnorm(5,mean = 30))
rownames(expr_df) <- paste0('Gene',1:5)
expr_df
> expr_df
Control_R1 Control_R2 Control_R3 Treat_R1 Treat_R2 Other_R1 Other_R2 Other_R3 Other_R4
Gene1 10.764324 9.428356 10.044180 20.54518 19.12997 29.55600 29.70370 31.64527 29.70927
Gene2 9.342808 9.345432 10.048242 21.24643 19.99557 29.66273 31.06946 28.56245 30.76298
Gene3 9.642546 8.777341 10.029212 18.77438 17.73295 29.96400 30.16638 29.81993 31.63248
Gene4 10.743541 8.764679 9.297456 20.41032 20.62758 29.48314 30.68438 31.92357 28.57588
Gene5 10.199939 9.814605 8.767410 20.91134 17.49593 28.59100 31.29823 30.26129 29.95323

这里我以5个基因，处理为Control（三个重复）、Treat（两个重复）、Other（四个重复）的表达量矩阵为例。

合并重复

# 得到对应的组织类型
> tissue <- colnames(expr_df) %>% gsub('_R[0-9]','',.)
> tissue
[1] 'Control' 'Control' 'Control' 'Treat'   'Treat'   'Other'   'Other'   'Other'   'Other'  
# 我们这里有三种处理类型
> tissue_type <- unique(tissue)
> tissue_type
[1] 'Control' 'Treat'   'Other'  
# 然后用sapply函数来合并重复
# 这里我用的是 -> 来赋值，大家不要学我……这个有点不太正规
sapply(tissue_type, function(x){
  rowSums(expr_df[,x == tissue])
}) -> expr_df_merge

> expr_df_merge
Control Treat Other
Gene1 30.23686 39.67515 120.6142
Gene2 28.73648 41.24201 120.0576
Gene3 28.44910 36.50733 121.5828
Gene4 28.80568 41.03790 120.6670
Gene5 28.78195 38.40728 120.1038

手动检查两个下对不对

> rowSums(expr_df[1,1:3])
   Gene1 
30.23686 
> rowSums(expr_df[3,1:3])
  Gene3 
28.4491

这里我稍微解释下这个合并重复的操作，具体的以后可以单独写下sapply、lapply这种操作。sapply其实是一个简化版的lapply，其输出的是向量，所以可以看到我们的表达量矩阵再也不是数据框了，而是matrix类型了。

> class(expr_df)
[1] 'data.frame'
> class(expr_df_merge)
[1] 'matrix'

然后这里传入的是 Control、Treat、Other 三种tissue_type，我们可以看下如果我们传入的是Control会怎么样。

# 假设x等于Control
> x <- 'Control'
# 这里就自动帮我们找到了Control所在的几个重复的位置，即前三个
> x == tissue
[1]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
# 从而提出了Conrol所在的几列
> expr_df[,x == tissue]
      Control_R1 Control_R2 Control_R3
Gene1  10.764324   9.428356  10.044180
Gene2   9.342808   9.345432  10.048242
Gene3   9.642546   8.777341  10.029212
Gene4  10.743541   8.764679   9.297456
Gene5  10.199939   9.814605   8.767410
# 然后用rowSums
# 这样就得到了Control的5个基因合并重复的表达量值了。
> rowSums(expr_df[,x == tissue])
   Gene1    Gene2    Gene3    Gene4    Gene5 
30.23686 28.73648 28.44910 28.80568 28.78195

求平均值

在合并完了重复之后，如果我们每个处理都是3个重复，那么其实我们不求平均值也行。但有时候事与愿违，可能某个处理下的某个重复由于质量不好，被我们剔除了，就变成了2个重复。所以为了应对这种情况，我把三个处理的重复数设置成3,2,4来举例子。

# 首先我们需要得到各个处理重复的频数
# 用的是table函数
> tissue %>% table()
.
Control Other Treat
3 4 2
> tissue_freq <- tissue %>% table() %>% as.numeric()
> tissue_freq
[1] 3 4 2

方法一

> t(t(expr_df_merge) / tissue_freq)
        Control     Treat    Other
Gene1 10.078953  9.918787 60.30712
Gene2  9.578828 10.310502 60.02881
Gene3  9.483033  9.126831 60.79140
Gene4  9.601892 10.259474 60.33349
Gene5  9.593985  9.601819 60.05188

我先来一步步地讲下方法一的原理，因为matrix其本质上还是（原子）向量，其只不过attribute里面多了dim这个属性，使得其变成了matrix这个特殊的数据结构。你可以想象成这里的matrix其实就是15个值，5个一叠，5个一叠，按Gene1……Gene5……Gene1……Gene5……Gene1……Gene5叠成了3列。

# 这里还多了dimnames这个属性
> typeof(expr_df_merge)
[1] 'double'
> attributes(expr_df_merge)
$dim
[1] 5 3
$dimnames
$dimnames[[1]]
[1] 'Gene1' 'Gene2' 'Gene3' 'Gene4' 'Gene5'
$dimnames[[2]]
[1] 'Control' 'Treat' 'Other'
# 双精度类型
> typeof(expr_df_merge)
[1] 'double'

关于原子向量、matrix这种Advanced_R这本书讲的很好，以后有机会也可以给大家讲讲。

然后R对于向量化的操作还具有短向量自动循环补全长向量的骚操作。比如

> 1:10
 [1]  1  2  3  4  5  6  7  8  9 10
> 1:2
[1] 1 2
> 1:10 - 1:2
 [1] 0 0 2 2 4 4 6 6 8 8
> 1:10 / 1:2
 [1] 1 1 3 2 5 3 7 4 9 5

这类1:2就自动补全了1:10的长度，所以你会看到1:10-1:2或者1:10 / 1:2 是这种结果。那么我们就可以利用这个骚操作来求我们的平均值。

# 一步步拆解下
> t(expr_df_merge)
Gene1 Gene2 Gene3 Gene4 Gene5
Control 30.23686 28.73648 28.44910 28.80568 28.78195
Treat 39.67515 41.24201 36.50733 41.03790 38.40728
Other 120.61424 120.05762 121.58279 120.66697 120.10376

这里我们转置了我们的合并重复的矩阵，那么matrix的堆叠顺序就是Control，Treat，Other，Control，Treat……。那么我们就可以根据上面提到的循环补齐操作来除对应的值了。

> t(expr_df_merge) / tissue_freq
            Gene1     Gene2     Gene3     Gene4     Gene5
Control 10.078953  9.578828  9.483033  9.601892  9.593985
Treat    9.918787 10.310502  9.126831 10.259474  9.601819
Other   60.307119 60.028809 60.791396 60.333487 60.051878

这里并不是我们想要的表达量矩阵，所以我们再次转置下。

> t(t(expr_df_merge) / tissue_freq)
Control Treat Other
Gene1 10.078953 9.918787 60.30712
Gene2 9.578828 10.310502 60.02881
Gene3 9.483033 9.126831 60.79140
Gene4 9.601892 10.259474 60.33349
Gene5 9.593985 9.601819 60.05188
# 再次手动检查下
> expr_df_merge[1,1] / 3
[1] 10.07895

方法二

> sweep(expr_df_merge, 2, tissue_freq, FUN = '/')
        Control     Treat    Other
Gene1 10.078953  9.918787 60.30712
Gene2  9.578828 10.310502 60.02881
Gene3  9.483033  9.126831 60.79140
Gene4  9.601892 10.259474 60.33349
Gene5  9.593985  9.601819 60.05188