I'm trying to use fuzzy string matching to convert strings to specific ids and perform grouped summarization using dplyr. The basic idea is combining imperfect gene sequences into a single gene name via a dictionary lookup approach and counting how many times the gene is detected. This way, counts for sequences aaaaaa and aaaxaa match to gene1 and get added together.
I can do what I want using for and if statements via a row-by-row comparison of the raw data against the dictionary but I find this will be inefficient when I scale up (raw data files have 15k rows on average, the dictionary has 200 rows). Please see my solution below I'm trying to improve and let me know if you can think of a more efficient and elegant way of doing this.
df <- data.frame(str_var = rep(c("aaaaaa", "aXaaaa", "bbbbbb", "bbbXbb"), 3),
grp_var = rep(c("grp1","grp2"), each=6),
num_var = rep(c(1,2), 6))
df
#> str_var grp_var num_var
#> 1 aaaaaa grp1 1
#> 2 aXaaaa grp1 2
#> 3 bbbbbb grp1 1
#> 4 bbbXbb grp1 2
#> 5 aaaaaa grp1 1
#> 6 aXaaaa grp1 2
#> 7 bbbbbb grp2 1
#> 8 bbbXbb grp2 2
#> 9 aaaaaa grp2 1
#> 10 aXaaaa grp2 2
#> 11 bbbbbb grp2 1
#> 12 bbbXbb grp2 2
dictionary <- data.frame(string = c("aaaaaa","bbbbbb", "cccccc", "dddddd"),
id = c("gene1", "gene2", "gene3", "gene4"))
dictionary
#> string id
#> 1 aaaaaa gene1
#> 2 bbbbbb gene2
#> 3 cccccc gene3
#> 4 dddddd gene4
for(i in 1:nrow(df)){
for(j in 1:nrow(dictionary)){
match_found <- agrepl(dictionary$string[j], df$str_var[i],
max.distance = list(sub=1, ins=0, del=0, all=1-1e-9))
if(match_found == TRUE){
gene = dictionary[j, "id"]
df$gene_id[i] <- gene
break
}
}
}
suppressPackageStartupMessages(library(dplyr))
new_df <- df %>%
group_by(grp_var, gene_id) %>%
summarize(gene_count=sum(num_var))
#> `summarise()` has grouped output by 'grp_var'. You can override using the `.groups` argument.
new_df
#> # A tibble: 4 x 3
#> # Groups: grp_var [2]
#> grp_var gene_id gene_count
#> <chr> <chr> <dbl>
#> 1 grp1 gene1 6
#> 2 grp1 gene2 3
#> 3 grp2 gene1 3
#> 4 grp2 gene2 6
Created on 2021-06-08 by the reprex package (v2.0.0)