Subset column and compute operations for each subset

Question

Here is a minimal example of dataframe to reproduce.

df <- structure(list(Gene = structure(c(147L, 147L, 148L, 148L, 148L, 
87L, 87L, 87L, 87L, 87L), .Label = c("genome", "k141_1189_101", 
"k141_1189_104", "k141_1189_105", "k141_1189_116", "k141_1189_13", 
"k141_1189_14", "k141_1189_146", "k141_1189_150", "k141_1189_18", 
"k141_1189_190", "k141_1189_194", "k141_1189_215", "k141_1189_248", 
"k141_1189_251", "k141_1189_252", "k141_1189_259", "k141_1189_274", 
"k141_1189_283", "k141_1189_308", "k141_1189_314", "k141_1189_322", 
"k141_1189_353", "k141_1189_356", "k141_1189_372", "k141_1189_373", 
"k141_1189_43", "k141_1189_45", "k141_1189_72", "k141_1597_15", 
"k141_1597_18", "k141_1597_23", "k141_1597_41", "k141_1597_55", 
"k141_1597_66", "k141_1597_67", "k141_1597_68", "k141_1597_69", 
"k141_2409_34", "k141_2409_8", "k141_3390_69", "k141_3390_83", 
"k141_3390_84", "k141_3726_25", "k141_3726_31", "k141_3726_49", 
"k141_3726_50", "k141_3726_62", "k141_3726_8", "k141_3726_80", 
"k141_3790_1", "k141_3993_114", "k141_3993_122", "k141_3993_162", 
"k141_3993_172", "k141_3993_183", "k141_3993_186", "k141_3993_188", 
"k141_3993_24", "k141_3993_25", "k141_3993_28", "k141_3993_32", 
"k141_3993_44", "k141_3993_47", "k141_3993_53", "k141_3993_57", 
"k141_3993_68", "k141_4255_80", "k141_4255_81", "k141_4255_87", 
"k141_5079_107", "k141_5079_110", "k141_5079_130", "k141_5079_14", 
"k141_5079_141", "k141_5079_16", "k141_5079_184", "k141_5079_185", 
"k141_5079_202", "k141_5079_24", "k141_5079_39", "k141_5079_63", 
"k141_5079_65", "k141_5079_70", "k141_5079_77", "k141_5079_87", 
"k141_5079_9", "k141_5313_16", "k141_5313_17", "k141_5313_20", 
"k141_5313_23", "k141_5313_39", "k141_5313_5", "k141_5313_51", 
"k141_5313_52", "k141_5313_78", "k141_5545_101", "k141_5545_103", 
"k141_5545_104", "k141_5545_105", "k141_5545_106", "k141_5545_107", 
"k141_5545_108", "k141_5545_109", "k141_5545_110", "k141_5545_111", 
"k141_5545_112", "k141_5545_113", "k141_5545_114", "k141_5545_119", 
"k141_5545_128", "k141_5545_130", "k141_5545_139", "k141_5545_141", 
"k141_5545_145", "k141_5545_16", "k141_5545_169", "k141_5545_17", 
"k141_5545_172", "k141_5545_6", "k141_5545_60", "k141_5545_62", 
"k141_5545_63", "k141_5545_86", "k141_5545_87", "k141_5545_88", 
"k141_5545_89", "k141_5545_91", "k141_5545_92", "k141_5545_93", 
"k141_5545_94", "k141_5545_96", "k141_5545_97", "k141_5545_98", 
"k141_5545_99", "k141_5734_13", "k141_5734_2", "k141_5734_4", 
"k141_5734_5", "k141_5734_6", "k141_6014_124", "k141_6014_2", 
"k141_6014_34", "k141_6014_75", "k141_6014_96", "k141_908_14", 
"k141_908_2", "k141_908_5", "k141_957_126", "k141_957_135", "k141_957_136", 
"k141_957_14", "k141_957_140", "k141_957_141", "k141_957_148", 
"k141_957_179", "k141_957_191", "k141_957_35", "k141_957_47", 
"k141_957_55", "k141_957_57", "k141_957_59", "k141_957_6", "k141_957_63", 
"k141_957_65", "k141_957_68", "k141_957_77", "k141_957_95"), class = "factor"), 
    depth = c(9L, 10L, 9L, 10L, 11L, 14L, 15L, 16L, 17L, 18L), 
    bases_covered = c(6L, 3L, 4L, 7L, 4L, 59L, 54L, 70L, 34L, 
    17L), gene_length = c(1140L, 1140L, 591L, 591L, 591L, 690L, 
    690L, 690L, 690L, 690L), regioncoverage = c(54L, 30L, 36L, 
    70L, 44L, 826L, 810L, 1120L, 578L, 306L)), .Names = c("Gene", 
"depth", "bases_covered", "gene_length", "regioncoverage"), row.names = c(1L, 
2L, 33L, 34L, 35L, 78L, 79L, 80L, 81L, 82L), class = "data.frame")

The dataframe looks like this:

      Gene depth bases_covered gene_length regioncoverage
1   k141_908_2     9             6        1140             54
2   k141_908_2    10             3        1140             30
33  k141_908_5     9             4         591             36
34  k141_908_5    10             7         591             70
35  k141_908_5    11             4         591             44
78 k141_5079_9    14            59         690            826
79 k141_5079_9    15            54         690            810
80 k141_5079_9    16            70         690           1120
81 k141_5079_9    17            34         690            578
82 k141_5079_9    18            17         690            306

What i want is that for each Gene (e.g k141_908_2) i want to sum region coverage and divide by unique(gene length). In fact gene length is always the same value for each gene.

For example for Gene K141_908_2 i would do: (54+30)/1140 = 0.07 For example for Gene K141_908_5 i would do: (36+70+44)/591 = 0.25

The final dataframe should report two columns.

     Gene Newcoverage
1   k141_908_2 0.07
2   k141_908_5 0.25
3   ......

and so on .

Thanks for your help

score 2 · Answer 1 · answered Feb 23 '17 at 15:34

This is straightforward with dplyr:

library(dplyr)
df_final <- df %>%
  group_by(Gene) %>%
  summarize(Newcoverage = sum(regioncoverage) / first(gene_length))

df_final
# # A tibble: 3 × 2
# Gene Newcoverage
# <fctr>       <dbl>
# 1 k141_5079_9  5.27536232
# 2  k141_908_2  0.07368421
# 3  k141_908_5  0.25380711

score 1 · Answer 2 · answered Feb 23 '17 at 15:34

I needed to set the first column to character and others to numeric. But after that you can just split the df by gene and then do the necessary calculations.

df[,2:5] = lapply(df[,2:5], as.numeric)
df$Gene = as.character(df$Gene)
sapply(split(df, df$Gene), function(x) sum(x[,5]/x[1,4]))
#k141_5079_9  k141_908_2  k141_908_5 
# 5.27536232  0.07368421  0.25380711

akrun · Answer 3 · 2017-02-23T15:43:21.190

1

We can use tidyverse

library(tidyverse)
df %>% 
  group_by(Gene) %>%
  summarise(Newcoverage = sum(regioncoverage)/gene_length[1])
# A tibble: 3 × 2
#          Gene Newcoverage
#        <fctr>       <dbl>
#1 k141_5079_9  5.27536232
#2  k141_908_2  0.07368421
#3  k141_908_5  0.25380711

Or a base R option is

by(df[4:5], list(as.character(df[,'Gene'])), FUN= function(x) sum(x[,2])/x[1,1])

edited Feb 23 '17 at 15:43

answered Feb 23 '17 at 15:35

akrun

874,273
37
540
662

user5249203 · Accepted Answer · 2017-02-23T20:45:42.277

0

quick approach is

require(data.table)
DT <- setDT(df)  
#just to output unique rows 
DT[, .(New_Coverage = unique(sum(regioncoverage)/gene_length)), by = .(Gene)]

output

         Gene New_Coverage
1:  k141_908_2   0.07368421
2:  k141_908_5   0.25380711
3: k141_5079_9   5.27536232

edited Feb 23 '17 at 20:45

answered Feb 23 '17 at 15:35

user5249203

4,436
1
19
45

Should be unique Gene (no duplicates in output). Also gene_length is unique. – david Feb 23 '17 at 15:47
updated for unique columns – user5249203 Feb 23 '17 at 15:53

score 0 · Answer 5 · answered Feb 23 '17 at 15:38

0

I use dplyr a lot. So here's one way:

library(dplyr)

df %>% 
  group_by(Gene) %>% 
  mutate(Newcoverage=sum(regioncoverage)/unique(gene_length))

If you want only unique values per Gene:

df %>% 
  group_by(Gene) %>% 
  transmute(Newcoverage=sum(regioncoverage)/unique(gene_length)) %>% 
  unique()

answered Feb 23 '17 at 15:38

Rahul

2,579
1
13
22

Great, thanks work perfect – david Feb 23 '17 at 15:57

Subset column and compute operations for each subset

5 Answers5