上一篇文章我提到了从表达量矩阵画时序性（或者多处理）的单基因的折线图，实际上我们的表达量矩阵里面每个处理是有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
# 这里是错误的……因为table会像ggplot2一样，会把输入字符串当因子，
# 然后根据字母顺序排序输出
> rep(c("a","c","b"),c(2,3,3))
[1] "a" "a" "c" "c" "c" "b" "b" "b"
> rep(c("a","c","b"),c(2,3,3)) %>% table()
.
a b c
2 3 3
# 当然，如果你输入的字符串本来就是因子，那么table并不会强制帮你更改
> test <- rep(c("a","c","b"),c(2,3,3))
> (test <- factor(test,levels = unique(test)))
[1] a a c c c b b b
Levels: a c b
> test %>% table()
.
a c b
2 3 3
# 所以我们这里在导出频数的时候应该多一步
# 这样才能与合并重复矩阵的列名相对应
> tissue %>% table() %>% "["(tissue_type)
.
Control Treat Other
3 2 4
> (tissue_freq <- tissue %>%
+ table() %>%
+ "["(tissue_type) %>%
+ as.numeric)
[1] 3 2 4
# "["这个其实就是[]这个函数
# 所以你也可以这样
> table(tissue)[tissue_type]
tissue
Control Treat Other
3 2 4
> as.numeric(table(tissue)[tissue_type])
[1] 3 2 4
# 或者提前转因子
> tissue <- factor(tissue,levels = tissue_type)
> table(tissue)
tissue
Control Treat Other
3 2 4
之前没有意识到table也会自动转因子，非常抱歉……
table强制转因子的顺序还是通用的数字字母排序，这意味着如果你的表达量矩阵是来自于featureCount求bam的count，在最后指定bam的时候，你直接用的*.bam，那么count的顺序就一般来说会和table强制转因子的顺序相同了。因为featureCount *.bam的顺序就是你ls *.bam的顺序，也是通用的数字字母顺序。
$ touch {a,c,b}.bam $ ll total 0 -rw-r--r--. 1 test bioinfo 0 Jan 20 20:03 a.bam -rw-r--r--. 1 test bioinfo 0 Jan 20 20:03 b.bam -rw-r--r--. 1 test bioinfo 0 Jan 20 20:03 c.bam $ ls *.bam a.bam b.bam c.bam这里是因为我自己构建的测试数据的缘故。但建议大家还是加一步[tissue_type]，或者直接对你的tissue转成因子再放入table，这样频数顺序就一致了。
unique函数似乎没有转因子的问题，大家可以放心用，其输出排序似乎是按照第一个出现的值
> test <- c("a","a","b","b","c","a","b") > test [1] "a" "a" "b" "b" "c" "a" "b" > unique(test) [1] "a" "b" "c" > test <- c("a","a","c","b","b","c","a","b") > unique(test) [1] "a" "c" "b"
方法一
> t(t(expr_df_merge) / tissue_freq)
Control Treat Other
Gene1 10.078953 19.83757 30.15356
Gene2 9.578828 20.62100 30.01440
Gene3 9.483033 18.25366 30.39570
Gene4 9.601892 20.51895 30.16674
Gene5 9.593985 19.20364 30.02594
我先来一步步地讲下方法一的原理，因为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.07895 9.578828 9.483033 9.601892 9.593985
Treat 19.83757 20.621004 18.253663 20.518948 19.203638
Other 30.15356 30.014404 30.395698 30.166743 30.025939
这里并不是我们想要的表达量矩阵，所以我们再次转置下。
> t(t(expr_df_merge) / tissue_freq)
Control Treat Other
Gene1 10.078953 19.83757 30.15356
Gene2 9.578828 20.62100 30.01440
Gene3 9.483033 18.25366 30.39570
Gene4 9.601892 20.51895 30.16674
Gene5 9.593985 19.20364 30.02594
# 再次手动检查下
> expr_df_merge[1,1] / 3
[1] 10.07895
方法二
> sweep(expr_df_merge, 2, tissue_freq, FUN = "/")
Control Treat Other
Gene1 10.078953 19.83757 30.15356
Gene2 9.578828 20.62100 30.01440
Gene3 9.483033 18.25366 30.39570
Gene4 9.601892 20.51895 30.16674
Gene5 9.593985 19.20364 30.02594
这里我们用到的是sweep函数
#Usage
# 这里x是数据，我们用的是合并了重复的矩阵
# MARGIN是维度，我们用2，即为列
# STATS是统计值，我们这里是个组织的重复频数
# FUN我们用除法
sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...)
# 其实就相当于每列除一个统计值，这里就是除我们的重复频数
其实sweep就类似于apply，具体的用法大家可以看这个R函数介绍：sweep()
感谢师兄对于求平均值过程中的指导╰(￣▽￣)╭
在勘误table的问题时候，我突然想到其实不需要先rowSum，然后再除频数。其实只需要RowMeans这个函数就可以了，由此又可以推广到求中位数，四分数等等……这个下一篇文章再说。
网友评论