使用 R 计算每百万映射读取的读取答案

【问题标题】：calculate reads per million mapped read using R使用 R 计算每百万映射读取的读取
【发布时间】：2015-09-04 12:31:52
【问题描述】：

df1 <- read.table(text="
   gene_id               A1      A2       A3     A4  length  Total
ENSMUSG00000000028       58      93       48     58   789     200                 
ENSMUSG00000000031       11      7        20     16   364     54                    
ENSMUSG00000000037       3       5         6     98   196     112                                       
ENSMUSG00000000058       66     93        69     71   436     299                                  
ENSMUSG00000000085       55     68        97     67   177     287", header=TRUE)

该表表示不同样本（A1、A2..A4）中基因的读取计数。如何使用 R

计算这些原始读取计数的每百万映射读取 (RPKM) 的读取

RPKM = (一个基因的读取数 * 1e6)/(总*长度)

out_put <-  read.table(text="
   gene_id               A1             A2             A3        A4  
ENSMUSG00000000028       367.5539      589.3536       304.1825    367.5539                   
ENSMUSG00000000031       559.6256      356.1254       1017.5010   814.0008                    
ENSMUSG00000000037       136.6618      227.7697       273.3236    4464.2857                                       
ENSMUSG00000000058       506.2747     713.3871        529.2872    544.6289                               
ENSMUSG00000000085       1082.6985     1338.6090      1909.4864   1318.9236", header=TRUE)

【问题讨论】：

你已经尝试过什么？为什么它不起作用？
@heroka.. 我很困惑如何使用这个函数编写 r 脚本
你是什么功能？你能展示一些预期的输出吗？查看您提供的数据和公式，您唯一缺少的是基因中的读取数（我认为是 A1+A2+A3+A4？）。
@Heroka，请找到预期的输出文件。
这取决于你。查看 dplyr 和 reshape2。

标签： r bioconductor

【解决方案1】：

不写行或循环的一种方法是使用melt和dcast：

library(reshape2)

m_df1 <- melt(df1, measure.vars=c("A1","A2","A3","A4"))
m_df1$RPKM <- with(m_df1, value*1e6 / (Total*length))

output <- dcast(gene_id~variable,value.var="RPKM",data=m_df1)
> output
             gene_id        A1        A2        A3        A4
1 ENSMUSG00000000028  367.5539  589.3536  304.1825  367.5539
2 ENSMUSG00000000031  559.6256  356.1254 1017.5010  814.0008
3 ENSMUSG00000000037  136.6618  227.7697  273.3236 4464.2857
4 ENSMUSG00000000058  506.2747  713.3871  529.2872  544.6289
5 ENSMUSG00000000085 1082.6985 1338.6090 1909.4864 1318.9236

第二种方法是使用 sapply 创建一个估计矩阵，然后您可以将其重命名并添加到您的原始数据中，或者 cbind 到您的gene_ids。

my_cols <- c("A1","A2","A3","A4")
RPKMs <- sapply(my_cols, function(x){
  df1[,x]*1e6/(df1$Total*df1$length)
}
)
output <- cbind(df1$gene_id,RPKMs)

【讨论】：

很好，+1。对于无需重塑的解决方案，请参阅我的答案。

【解决方案2】：

您也可以在不重塑的情况下实现这一目标。使用data.table 包：

library(data.table)
setDT(df1)[,indx:=.I][, lapply(.SD, function(x) (x * 1e6) / (Total * length)),
                      by=.(indx,gene_id,length,Total)]

这给出了：

   indx            gene_id length Total        A1        A2        A3        A4
1:    1 ENSMUSG00000000028    789   200  367.5539  589.3536  304.1825  367.5539
2:    2 ENSMUSG00000000031    364    54  559.6256  356.1254 1017.5010  814.0008
3:    3 ENSMUSG00000000037    196   112  136.6618  227.7697  273.3236 4464.2857
4:    4 ENSMUSG00000000058    436   299  506.2747  713.3871  529.2872  544.6289
5:    5 ENSMUSG00000000085    177   287 1082.6985 1338.6090 1909.4864 1318.9236

解释：

使用setDT(df1) 将数据帧转换为数据表
使用[,indx:=.I] 为每一行创建一个唯一标识符
使用by=.(indx,gene_id,length,Total)，您可以确定要对数据进行分组的列（这些列不会被转换），通过包含indx，您可以确保每一行都是一个唯一的组
使用lapply(.SD, function(x) (x * 1e6) / (Total * length))，您可以将所需的计算应用于by 语句中未指定的每一列

与dplyr 类似的解决方案：

library(dplyr)

func <- function(x,y,z) (x * 1e6) / (y * z)

df1 %>% mutate(indx=seq(1,nrow(.))) %>% 
  group_by(indx,gene_id,length,Total) %>% 
  summarise_each(funs(func(.,Total,length)))

给出：

   indx            gene_id length Total        A1        A2        A3        A4
  (int)             (fctr)  (int) (int)     (dbl)     (dbl)     (dbl)     (dbl)
1     1 ENSMUSG00000000028    789   200  367.5539  589.3536  304.1825  367.5539
2     2 ENSMUSG00000000031    364    54  559.6256  356.1254 1017.5010  814.0008
3     3 ENSMUSG00000000037    196   112  136.6618  227.7697  273.3236 4464.2857
4     4 ENSMUSG00000000058    436   299  506.2747  713.3871  529.2872  544.6289
5     5 ENSMUSG00000000085    177   287 1082.6985 1338.6090 1909.4864 1318.9236

【讨论】：