Подсчитать элементы списка и сохранить только те, которые появляются x раз

Question

Это немного похоже на этот вопрос о подсчете вхождений строк, но мне нужно решение, которое работает с dplyr и способно применять тесты на основе результатов.

В геномике естьидея k-мер, которая является каждой подстрокой длины k внутри более длинной строки.Мне нужно взять несколько строк (геномных последовательностей) и найти (потенциально перекрывающиеся) k-мер длины 5, которые встречаются внутри строки ровно 4 раза.

Я могу подсчитать количество вхождений каждого k-мер с помощью следующего:

sequence1 <- "CGGACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAA"
sequence2 <- "GGACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAAT"
sequence3 <- "GACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAATG"
sequences <- list(sequence1, sequence2, sequence3)

#Generate all k-mers of length 5 within each sequence
k <- 5
kmers <- map(sequences, function(x) {
    map_chr(seq_len(nchar(x) - k + 1), function(y) str_sub(x, y, y + k - 1))}) %>% 
    set_names(sequences)

kmers

Дает k-мер:

#> $CGGACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAA
#>  [1] "CGGAC" "GGACT" "GACTC" "ACTCG" "CTCGA" "TCGAC" "CGACA" "GACAG"
#>  [9] "ACAGA" "CAGAT" "AGATG" "GATGT" "ATGTG" "TGTGA" "GTGAA" "TGAAG"
#> [17] "GAAGA" "AAGAA" "AGAAA" "GAAAT" "AAATG" "AATGT" "ATGTG" "TGTGA"
#> [25] "GTGAA" "TGAAG" "GAAGA" "AAGAC" "AGACT" "GACTG" "ACTGA" "CTGAG"
#> [33] "TGAGT" "GAGTG" "AGTGA" "GTGAA" "TGAAG" "GAAGA" "AAGAG" "AGAGA"
#> [41] "GAGAA" "AGAAG" "GAAGA" "AAGAG" "AGAGG" "GAGGA" "AGGAA" "GGAAA"
#> [49] "GAAAC" "AAACA" "AACAC" "ACACG" "CACGA" "ACGAC" "CGACA" "GACAC"
#> [57] "ACACG" "CACGA" "ACGAC" "CGACA" "GACAT" "ACATT" "CATTG" "ATTGC"
#> [65] "TTGCG" "TGCGA" "GCGAC" "CGACA" "GACAT" "ACATA" "CATAA"
#> 
#> $GGACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAAT
#>  [1] "GGACT" "GACTC" "ACTCG" "CTCGA" "TCGAC" "CGACA" "GACAG" "ACAGA"
#>  [9] "CAGAT" "AGATG" "GATGT" "ATGTG" "TGTGA" "GTGAA" "TGAAG" "GAAGA"
#> [17] "AAGAA" "AGAAA" "GAAAT" "AAATG" "AATGT" "ATGTG" "TGTGA" "GTGAA"
#> [25] "TGAAG" "GAAGA" "AAGAC" "AGACT" "GACTG" "ACTGA" "CTGAG" "TGAGT"
#> [33] "GAGTG" "AGTGA" "GTGAA" "TGAAG" "GAAGA" "AAGAG" "AGAGA" "GAGAA"
#> [41] "AGAAG" "GAAGA" "AAGAG" "AGAGG" "GAGGA" "AGGAA" "GGAAA" "GAAAC"
#> [49] "AAACA" "AACAC" "ACACG" "CACGA" "ACGAC" "CGACA" "GACAC" "ACACG"
#> [57] "CACGA" "ACGAC" "CGACA" "GACAT" "ACATT" "CATTG" "ATTGC" "TTGCG"
#> [65] "TGCGA" "GCGAC" "CGACA" "GACAT" "ACATA" "CATAA" "ATAAT"
#> 
#> $GACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAATG
#>  [1] "GACTC" "ACTCG" "CTCGA" "TCGAC" "CGACA" "GACAG" "ACAGA" "CAGAT"
#>  [9] "AGATG" "GATGT" "ATGTG" "TGTGA" "GTGAA" "TGAAG" "GAAGA" "AAGAA"
#> [17] "AGAAA" "GAAAT" "AAATG" "AATGT" "ATGTG" "TGTGA" "GTGAA" "TGAAG"
#> [25] "GAAGA" "AAGAC" "AGACT" "GACTG" "ACTGA" "CTGAG" "TGAGT" "GAGTG"
#> [33] "AGTGA" "GTGAA" "TGAAG" "GAAGA" "AAGAG" "AGAGA" "GAGAA" "AGAAG"
#> [41] "GAAGA" "AAGAG" "AGAGG" "GAGGA" "AGGAA" "GGAAA" "GAAAC" "AAACA"
#> [49] "AACAC" "ACACG" "CACGA" "ACGAC" "CGACA" "GACAC" "ACACG" "CACGA"
#> [57] "ACGAC" "CGACA" "GACAT" "ACATT" "CATTG" "ATTGC" "TTGCG" "TGCGA"
#> [65] "GCGAC" "CGACA" "GACAT" "ACATA" "CATAA" "ATAAT" "TAATG"

ИЯ могу найти счет с помощью

kmers %>%
    imap(~ str_count(.y, .x))

, который возвращает

#> $CGGACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAA
#>  [1] 1 1 1 1 1 1 4 1 1 1 1 1 2 2 3 3 4 1 1 1 1 1 2 2 3 3 4 1 1 1 1 1 1 1 1
#> [36] 3 3 4 2 1 1 1 4 2 1 1 1 1 1 1 1 2 2 2 4 1 2 2 2 4 2 1 1 1 1 1 1 4 2 1
#> [71] 1
#> 
#> $GGACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAAT
#>  [1] 1 1 1 1 1 4 1 1 1 1 1 2 2 3 3 4 1 1 1 1 1 2 2 3 3 4 1 1 1 1 1 1 1 1 3
#> [36] 3 4 2 1 1 1 4 2 1 1 1 1 1 1 1 2 2 2 4 1 2 2 2 4 2 1 1 1 1 1 1 4 2 1 1
#> [71] 1
#> 
#> $GACTCGACAGATGTGAAGAAATGTGAAGACTGAGTGAAGAGAAGAGGAAACACGACACGACATTGCGACATAATG
#>  [1] 1 1 1 1 4 1 1 1 1 1 2 2 3 3 4 1 1 1 1 1 2 2 3 3 4 1 1 1 1 1 1 1 1 3 3
#> [36] 4 2 1 1 1 4 2 1 1 1 1 1 1 1 2 2 2 4 1 2 2 2 4 2 1 1 1 1 1 1 4 2 1 1 1
#> [71] 1

Но теперь мне нужно вернуть каждый уникальный k-мер, который имеет счет 4. Мое единственное решение покавзять индекс каждого k-mer, равный 4, а затем регенерировать k-mer с помощью substr следующим образом:

kmers >%>
    imap(~ str_count(.y, .x)) %>%
#test for k-mers that appear 4 times
    map(function(y) {
        map_lgl(y, function(x) x == 4)}) %>%
#Get the indexes of the matches
    map(which) %>%
#Recreate the k-mers from each sequence
    imap(function(a,b) {
        map_chr(a, ~ substr(b, .x, .x + k -1))}) %>%
    unlist %>%
    unique

Что дает мне желаемый результат

#>  'CGACA' 'GAAGA' 

Но это неэффективновыбросить k-mers и затем воссоздать их.Как я могу получить счетчики, но затем использовать их для фильтрации исходного списка k-мер?table() работает, но я не могу понять, как работать с table(), приводит к конвейеру dplyr, чтобы получить результат с простым списком строк, удовлетворяющих условию.

Answer 1

Я принял ответ Ронака Шаха, потому что он обладал замечательным очевидным свойством сохранять счетчики в отдельном списке и затем использовать этот список для проверки исходного списка k-мер. Однако я продолжал пытаться обеспечить поток данных через один конвейер без необходимости сохранять промежуточные значения. Вот где я получил

map(sequences, function(x) {
    #Generate all k-mers for each sequence
    str_sub(x, a <- seq_len(nchar(x) - k + 1), a + k - 1) %>%
    #Count how many times each k-mer appears in the sequence and keep if equal to 4
    keep(str_count(x, .) == 4)}) %>%
unlist %>%
unique

Answer 2

Вот полностью базовый способ.

Изменить: База kmer вызов был изменен.Вы должны рассмотреть возможность использования другого подхода к кмер.Подход base r примерно в 100 раз быстрее.Если вы удалите set_names, база r будет еще в 20 раз быстрее.Я также упростил вызов lapply(kmers, function (x) ...).

kmers <- lapply(sequences, function(x) substring(x, seq_len(nchar(x) - k + 1), seq_len(nchar(x) - k + 1)+ k - 1))
names(kmers) = sequences

unique(
  unlist(
    lapply(kmers, function(x) names(Filter(function(z) z == 4, table(x))))
    )
  )

[1] "CGACA" "GAAGA"

#Or with no intermediate variables (it doesn't look pretty)
  unique(
    unlist(
      lapply(sequences, function(x) names(Filter(function(z) z == 4, table(substring(x, seq_len(nchar(x) - k + 1), seq_len(nchar(x) - k + 1)+ k - 1)))))
    )
  )

# Or same thing with chains (you can use map instead of lapply):

  lapply(sequences, function(x) {
    table(substring(x, seq_len(nchar(x) - k+1), seq_len(nchar(x) - k+1) + k-1))%>%
      Filter(function(z) z == 4, .)%>%
      names()
    }
  )%>%
    unlist()%>%
    unique()

Производительность альтернативного строкового метода:

Unit: microseconds
               expr    min      lq     mean  median      uq     max neval
         base_kmers   50.0   52.50  129.598   58.15   71.85  6836.7   100
       base_w_names   52.5   55.85  135.946   62.10   77.65  6995.0   100
 purr_w_names_kmers 1651.1 7259.55 7683.990 7569.30 7898.35 10476.2   100
         purr_kmers 1260.8 1294.65 1424.776 1322.50 1364.00  7395.9   100

Производительность фильтра - мое решение примерно в два раза быстрее, чем у @ Ronak.Решение OP примерно в 20 раз медленнее, хотя оно и пропускает создание переменной kmers.

Unit: microseconds
                     expr     min       lq      mean   median       uq     max neval
          base_everything   953.9   991.80  1201.355  1030.50  1077.40 10953.1   100
     base_no_intermediate   946.2   972.20  1196.071  1017.50  1107.00 11104.0   100
 base_no_inter_plus_chain 11437.7 11970.35 13003.189 12341.15 12835.15 43925.7   100
     base_kmers_and_ronak  1964.8  2053.60  2426.282  2120.90  2342.00 10817.2   100
                OP_answer  7204.6 18481.50 19349.587 19095.95 20426.25 22423.6   100
             Ony_solution   721.2   748.55   823.104   766.15   799.90  4803.8   100

Answer 3

Вы можете flatten kmers и cnts и подмножество kmers, где cnts равно 4. Вы можете сгладить значения здесь, так как в конечном итоге вам нужны только unique.

library(tidyverse)

cnts <- kmers %>% imap(~ str_count(.y, .x))
(kmers %>% flatten_chr)[cnts %>% flatten_int == 4] %>% unique
#[1] "CGACA" "GAAGA"

, что совпадает с

unique(unlist(kmers)[unlist(cnts) == 4])