0

我是 R 编程的新手，我正在尝试为 Reverse and Complementary Base 编写一个程序。目的是设计 DNA 引物。所以我有一个 DNA 序列，其碱基为 ATCG 和 A 与 T 互补；T=A；C=G；G=C。我已经想出了如何反转它，但是对于补码，我只能让它回答 1 个碱基但不能是所有序列，而且我不知道如何组合反转和补码功能。这是我的代码，我完全对它感到困惑。有人可以帮我解决这个问题吗？你会是我的救命恩人！

strReverse <- function(x) 

sapply(lapply(strsplit(x, NULL), rev), paste, collapse="")
strReverse(c("ATCGGTCAATCGA"))

complement.base = function(base){ 
  if(base == 'A' | base ==  'a')   print("T") 
  if(base == 'T' | base == 't')         print("A")
  if(base == 'G' | base == 'g')     print("C")
  if(base == 'C' | base == 'c')     print("G")}
complement.base(base="A")
4

2 回答 2

2

您可以使用 Rcpp 有效地执行操作：

library(Rcpp)
revComp.rcpp <- cppFunction(
"std::string comp(std::string x) {
  const int n = x.length();
  for (int i=0; i < n; ++i) {
    if (x[i] == 'A' || x[i] == 'a')  x[i] = 'T';
    else if (x[i] == 'T' || x[i] == 't')  x[i] = 'A';
    else if (x[i] == 'G' || x[i] == 'g')  x[i] = 'C';
    else  x[i] = 'G';
  }
  std::reverse(x.begin(), x.end());
  return x;
}")
revComp.rcpp("ATCGGTCAATCGA")
# [1] "TCGATTGACCGAT"

这似乎比 Biostrings 包中的相关代码要快一些（在具有 1300 万个碱基的字符串上进行测试）：

library(Biostrings)
x <- "ATCGGTCAATCGA"
big.x <- paste(rep(x, 1000000), collapse="")
big.x2 <- DNAString(big.x)
rev.biostr <- function(x) as.character(reverseComplement(x))
all.equal(revComp.rcpp(big.x), as.character(reverseComplement(big.x2)))
# [1] TRUE

library(microbenchmark)
microbenchmark(revComp.rcpp(big.x), as.character(reverseComplement(big.x2)))
# Unit: milliseconds
#                                     expr       min        lq      mean    median        uq      max neval
#                      revComp.rcpp(big.x)  77.21618  78.44534  84.54397  82.21002  87.49367 123.8166   100
#  as.character(reverseComplement(big.x2)) 144.13900 151.12869 170.73765 156.44300 164.41374 399.2948   100
于 2015-03-09T02:19:06.030 回答
1

我实际上会考虑使用chartrbase R，并在一些帮助stringi下反转结果（或输入）。

myFun <- function(invec) {
  require(stringi)
  invec <- stri_reverse(invec)
  chartr(old = "AaTtGgCc", new = "TTAACCGG", invec)
}

x <- "ATCGGTCAATCGA"
myFun(x)
# [1] "TCGATTGACCGAT"

使用@josilber 的样本数据，它与他的 Rcpp 方法非常相似：

all.equal(myFun(big.x), revComp.rcpp(big.x))
# [1] TRUE

library(microbenchmark)
microbenchmark(myFun(big.x), revComp.rcpp(big.x))
# Unit: milliseconds
#                 expr      min       lq     mean   median       uq      max neval
#         myFun(big.x) 349.5797 352.8197 362.3009 356.4484 362.7197 437.9556   100
#  revComp.rcpp(big.x) 359.5485 363.8615 378.3465 368.3360 386.3734 444.2901   100
于 2015-03-09T03:32:07.133 回答