在没有任何包的情况下查找基因组组合

Question

我想知道在一个序列中有多少种基因组组合。我的意思是二进制组合：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

Answer 1

这是一种缓慢的方法。我确信它在 bioconductor 包中更快。

# 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)

Answer 2

简单的事情就是这样：

# 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

这里有一个使用 "function factory" (https://adv-r.hadley.nz/function-factories.html) 的方法。

2 元素和 3 元素的组合是大小为 2 和 3 的 n-gram。所以我们将这个 n-gram 函数工厂化。

# 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))

在没有任何包的情况下查找基因组组合

Find combinations of genome without any package

r

seq