Найти комбинации генома без какой-либо упаковки - PullRequest
Новая
       7

Найти комбинации генома без какой-либо упаковки

0 голосов
/ 21 января 2019

Я хочу узнать, сколько комбинаций генома найдено в последовательности.Я имею в виду бинарные комбинации: AA, AT, AG, AC, ... 16 таких комбинаций или для 3-элементных комбинаций ATG, ACG, ... 64 таких комбинаций.Я знаю, как сделать это с пакетом, и я напишу это здесь.Я хочу создать свой собственный код для выполнения этого пакета

seqinr, идеально подходящего для своей работы.Это код, который я использовал для 

install.packages('seqinr')    
library(seqinr)    
m = read.fasta(file='sequence.fasta')     
mseq = m[[1]]     
count(mseq,2)   # gives how many binary combinations are found in the seq     
count(mseq,3)   # gives how many 3-elemented combinations are found in the seq

Ответы [ 2 ]

0 голосов
/ 21 января 2019

Самое простое, что нужно сделать, это примерно так: 

# Generate a test sequence
set.seed(1234)
testSeq <- paste(sample(LETTERS[1:3], 100, replace = T), collapse = "")

# Split string into chunks of size 2 and then count occurrences
testBigram <- substring(testSeq, seq(1, nchar(testSeq), 2), seq(2, nchar(testSeq), 2))
table(testBigram)

AA AB AC BA BB BC CA CB CC 
10 10 14  3  3  2  2  5  1

Вот способ использования "фабрики функций" (https://adv -r.hadley.nz / function-factories.html ).

2-элементные и 3-элементные комбинации - это n-граммы размера 2 и 3. Таким образом, мы делаем эту фабрику n-граммных функций. 

# Generate a function to create a function
ngram <- function(size) {
  function(myvector) {

    substring(myvector, seq(1, nchar(myvector), size), seq(size, nchar(myvector), size))

  }
}

# Assign the functions names (optional)
bigram <- ngram(2)
trigram <- ngram(3)

# 2 element combinations
table(bigram(testSeq))

AA AB AC BA BB BC CA CB CC 
10 10 14  3  3  2  2  5  1 

# count of 2 element combinations
length(unique(bigram(testSeq)))

[1] 9

# counting function
count <- function(mseq, n) length(unique(ngram(n)(mseq)))
count(testSeq, 2)

[1] 9

# and if we wanted to do with with 3 element combinations
table(trigram(testSeq))
0 голосов
/ 21 января 2019

Это медленный способ сделать это. Я уверен, что это быстрее в пакете биокондуктора. 

# some practice data
mseq = paste(sample(c("A", "C", "G", "T"), 1000, rep=T), collapse="")

# define a function called count
count = function(mseq, n){
  # split the sequence into every possible sub sequence of length n
  x = sapply(1:(nchar(mseq) - n + 1), function(i) substr(mseq, i, i+n-1))
  # how many unique sub sequences of length R are there?
  length(table(x))
}

На самом деле только что проверил, и вот как они это сделали: 

function (seq, wordsize, start = 0, by = 1, freq = FALSE, alphabet = s2c("acgt"), 
    frame = start) 
{
    if (!missing(frame)) 
        start = frame
    istarts <- seq(from = 1 + start, to = length(seq), by = by)
    oligos <- seq[istarts]
    oligos.levels <- levels(as.factor(words(wordsize, alphabet = alphabet)))
    if (wordsize >= 2) {
        for (i in 2:wordsize) {
            oligos <- paste(oligos, seq[istarts + i - 1], sep = "")
        }
    }
    counts <- table(factor(oligos, levels = oligos.levels))
    if (freq == TRUE) 
        counts <- counts/sum(counts)
    return(counts)
}

Если вы хотите найти код для функции, используйте getAnywhere () 

getAnywhere(count)
...