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

Question

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)

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

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

RPKM =（一个基因的reads数*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)

Answer 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)

Answer 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

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

calculate reads per million mapped read using R

r

bioconductor