Как посчитать строку букв в кмер

Question

Я пытаюсь рассчитать долю населения 5-меров в 100 000 нуклеотидов, 5-меров AATAA. Сколько AATAA повторяется в данных.

dog_ch38 <- read.GenBank("NC_006620.3")
dog_ch38 <- dog_ch38$NC_006620.3[1:100000]
dog_ch38 <- c2s(ape::as.character.DNAbin(dog_ch38))
dog_ch38 <- str_to_upper(dog_ch38)
 kmer_to_index <- function(kmer){
+     n  <- str_length(kmer)
+     letter_value  <- c("A" = 0, "C" = 1, "G" = 2, "T" = 3)
+     base  <- 1
+     index  <- 1
+     for( i in n:1){
+         nucleotide <- str_sub(kmer,start = i,end = i)
+         index  <- index + base * letter_value[nucleotide]
+         base  <- base * 4
+     }
+     return(as.numeric(index))
+ }
k <- 5
kmers <- numeric(4^k)
kmers
N <- str_length(dog_ch38)
> for (i in 1:(N - k + 1)) {
+   kmer <- str_sub(dog_ch38, i, i + k - 1)
+   index <- kmer_to_index(kmer)
+   kmers[index] <- kmers[index] + 1
+ }

но получаю ошибку вот так:

Error in kmers[index] <- kmers[index] + 1 : 
  NAs are not allowed in subscripted assignments
In addition: Warning messages:
1: In 1:(N - k + 1) :
  numerical expression has 100000 elements: only the first used
2: In n:1 : numerical expression has 100000 elements: only the first used

seqinr::count(dog_ch38[1,], 5)
Error in dog_ch38[1, ] : incorrect number of dimensions

что я на самом деле ожидал увидеть в результате:

## aaaaa accaa taaag aataa
##   75   75     47    92 

Я определенно новичок в этой функции, и если кто-нибудь может подсказать мне, как ее решить, и некоторые примеры, чтобы посмотреть. Спасибо!

Answer 1

Я думаю oligonucleotideFrequency() функция от Biostrings Пакет может помочь. Вот пример с искусственными данными.

  library(Biostrings) # requires appropriate Bioconductor install
  s1 <- sample(c("A", "C", "G", "T"), 10^5, TRUE)
  s1 <- DNAString(paste(s1, collapse = ""))
  kmers <- oligonucleotideFrequency(s1, width = 5)

Функция вернула именованный числовой вектор всех возможных кмеров. Вы можете использовать эту функцию, чтобы извлечь интересующий вас кмер. Там должно быть около 100 в этом примере.

  kmers["AATAA"] # actual count varies because of random sampling
> AATAA 
>  102 

Проверьте страницу справки для этой функции. С аргументами по умолчанию он вернет перекрывающиеся kmers. Это можно контролировать с помощью step опция, как показано в этом примере:

  s2 <- DNAString("AATAATAATAA")
  kmers1 <- oligonucleotideFrequency(s2, width = 5)
  kmers2 <- oligonucleotideFrequency(s2, width = 5, step = 5)

# See all the 5-mers found with step = 1 (default) versus step = 5
  kmers1[kmers1 != 0]
> AATAA ATAAT TAATA 
>     3     2     2 
  kmers2[kmers2 != 0]
> AATAA TAATA 
>     1     1

РЕДАКТИРОВАТЬ

Я нашел (все еще нахожу) разнообразие форматов для обработки последовательностей ДНК, сбивающее с толку, и кажется, что есть необходимость преобразовать компактную двоичную форму, возвращаемую read.GenBank() на представление персонажа в Biostrings, Они оба удивительно эффективны.

Преобразование может быть сделано в списке двоичных объектов, которые read.GenBank() возвращается или вы можете использовать as.character = TRUE возможность вернуть необработанные символы. Я показываю последний подход здесь.

# Using package ape to read GenBank file, Biostrings for analysis
  library(ape)
  library(Biostrings)

# By default, read.GenBank returns a list of DNA sequences in compact binary form.
# This asks it to return a list of character vectors. 
  dog_ch38 <- read.GenBank("NC_006620.3", as.character = TRUE)
  str(dog_ch38)
> List of 1
>  $ NC_006620.3: chr [1:23914537] "n" "n" "n" "n" ...
>  - attr(*, "species")= chr "Canis_lupus_familiaris"

# Now convert the first (and only) member of the list to a single character string
  txt <- paste(dog_ch38[[1]], collapse = "")
  print(nchar(txt))
> [1] 23914537

# And now convert the character string to a DNAString
  s <- DNAString(txt)

# This is the form that can be handed to oligonucleotideFrequency
  km <- oligonucleotideFrequency(s[1:10^5], 5)
  km["AATAA"]
> AATAA 
>   176 

Answer 2

Вадим прав! Незначительное обновление: вам не нужен этот фрагмент кода:

      # Now convert the first (and only) member of the list to a single character string
  txt <- paste(dog_ch38[[1]], collapse = "")

для вашего кода, потому что вы уже использовали функции c2s() и str_to_upper(). Вот простой метод, который работает:

      #load libraries
library(ape)
library(tidyverse)
library(Biostrings) # installed by 1) install.packages("BiocManager") and 2) BiocManager::install("Biostrings")

dog_ch38 <- read.GenBank("NC_006620.3")

#first 100k nucleotides
dog_ch38_sub <- dog_ch38$NC_006620.3[1:100000] # the subset
print(nchar(dog_ch38_sub) # verify you have 100k

#convert to a form that Biostrings understands
dog_ch38_sub <- c2s(ape::as.character.DNAbin(dog_ch38_sub))
dog_ch38_sub <- str_to_upper(dog_ch38_sub)
dog_ch38_sub <- DNAString(dog_ch38_sub)

km <- oligonucleotideFrequency(dog_ch38_sub[1:100000], 5)
count_AATAA <- km["AATAA"]

#your proportion
count_AATAA / nchar(dog_ch38_sub)

Спасибо.