补充 DNA 序列答案

【问题标题】：Complement a DNA sequence补充 DNA 序列
【发布时间】：2013-12-20 17:34:23
【问题描述】：

假设我有一个 DNA 序列。我想得到它的补充。我使用了以下代码，但我没有得到它。我做错了什么？

s=readline()
ATCTCGGCGCGCATCGCGTACGCTACTAGC
p=unlist(strsplit(s,""))
h=rep("N",nchar(s))
unlist(lapply(p,function(d){
for b in (1:nchar(s)) {    
    if (p[b]=="A") h[b]="T"
    if (p[b]=="T") h[b]="A"
    if (p[b]=="G") h[b]="C"
    if (p[b]=="C") h[b]="G"
}

【问题讨论】：

标签： r replace bioinformatics genetics complement

【解决方案1】：

使用为此目的而构建的chartr：

> s
[1] "ATCTCGGCGCGCATCGCGTACGCTACTAGC"
> chartr("ATGC","TACG",s)
[1] "TAGAGCCGCGCGTAGCGCATGCGATGATCG"

只需给它两个等长的字符串和你的字符串。还对翻译参数进行了矢量化：

> chartr("ATGC","TACG",c("AAAACG","TTTTT"))
[1] "TTTTGC" "AAAAA"

请注意，我正在替换 DNA 的字符串表示，而不是向量。要转换向量，我会创建一个查找图作为命名向量和索引：

> p
 [1] "A" "T" "C" "T" "C" "G" "G" "C" "G" "C" "G" "C" "A" "T" "C" "G" "C" "G" "T"
[20] "A" "C" "G" "C" "T" "A" "C" "T" "A" "G" "C"
> map=c("A"="T", "T"="A","G"="C","C"="G")
> unname(map[p])
 [1] "T" "A" "G" "A" "G" "C" "C" "G" "C" "G" "C" "G" "T" "A" "G" "C" "G" "C" "A"
[20] "T" "G" "C" "G" "A" "T" "G" "A" "T" "C" "G"

【讨论】：

【解决方案2】：

Bioconductor 包Biostrings 为此类操作提供了许多有用的功能。安装一次：

source("http://bioconductor.org/biocLite.R")
biocLite("Biostrings")

然后使用

library(Biostrings)
dna = DNAStringSet(c("ATCTCGGCGCGCATCGCGTACGCTACTAGC", "ACCGCTA"))
complement(dna)

【讨论】：

【解决方案3】：

为了补充，无论大小写，您都可以使用chartr()：

n <- "ACCTGccatGCATC"
chartr("acgtACGT", "tgcaTGCA", n)
# [1] "TGGACggtaCGTAG"

为了更进一步并反向互补核苷酸序列，您可以使用以下功能：

library(stringi)

rc <- function(nucSeq)
  return(stri_reverse(chartr("acgtACGT", "tgcaTGCA", nucSeq)))

rc("AcACGTgtT")
# [1] "AacACGTgT"

【讨论】：

【解决方案4】：

sapply(p, switch,  "A"="T", "T"="A","G"="C","C"="G")
  A   T   C   T   C   G   G   C   G   C   G   C   A   T   C   G   C   G   T 
"T" "A" "G" "A" "G" "C" "C" "G" "C" "G" "C" "G" "T" "A" "G" "C" "G" "C" "A" 
  A   C   G   C   T   A   C   T   A   G   C 
"T" "G" "C" "G" "A" "T" "G" "A" "T" "C" "G"

如果您不想要互补名称，您可以随时使用unname 去除它们。

unname(sapply(p, switch,  "A"="T", "T"="A","G"="C","C"="G") )
 [1] "T" "A" "G" "A" "G" "C" "C" "G" "C" "G" "C" "G" "T" "A" "G" "C" "G" "C"
[19] "A" "T" "G" "C" "G" "A" "T" "G" "A" "T" "C" "G"
>

【讨论】：

太好了，我从来不知道switch函数。
似乎几乎是为此目的而定制的。
也许是由没有审美意识的人定制的！ :)

【解决方案5】：

还有一个包seqinr

library(seqinr)
comp(seq) # gives complement
rev(comp(seq)) # gives the reverse complement

Biostrings 的内存配置要小得多，但 seqinr 也很好，因为您可以选择碱基的大小写（包括混合）并将它们更改为您想要的任何内容，例如，如果您希望在相同的顺序。 Biostrings 强制您使用 T 或 U。

【讨论】：

我不确定这是否完全准确：comp("u") 返回 NA
@TomKelly 您的序列必须采用字符向量的形式，每个字符的长度为 1。JeremyS 忘记提及的是，如果它是单个字符，您需要先剪切序列.重写，这应该工作：comp(s2c(seq))
@YoannPageaud 实际上，我在这里添加了对字符串（任意长度）的支持：github.com/TomKellyGenetics/tktools/blob/master/R/…
@YoannPageaud 不，问题是不支持“U”。与“A”或“T”相同的输入有效。 ``` > library(seqinr) > comp("t") [1] "a" > comp("u") [1] NA > comp("a") [1] "t" > comp(c( "u")) [1] NA > comp(c("U")) [1] NA > comp(c("T")) [1] "a" > comp(s2c("u")) [ 1] 不适用```

【解决方案6】：

这里是使用底数 r 的答案。以一种可怕的格式编写，以使事情清晰并保持单行。支持大小写。

revc = function(s){
       paste0(
           rev(
            unlist(
             strsplit(
                chartr("ATGCatgc","TACGtacg",s)
                      , "")                        # from strsplit
                   )                               # from unlist
               )                                   # from rev
             , collapse='')                        # from paste0
       }

【讨论】：

【解决方案7】：

我已经用seqinr 包概括了解决方案rev(comp(seq))：

install.packages("devtools")
devtools::install_github("TomKellyGenetics/tktools")
tktools::revcomp(seq)

此版本与字符串输入兼容，并且被矢量化以处理多个字符串的列表或向量输入。输出类应与输入匹配，包括案例和类型。这也支持包含“U”的输入，用于 RNA 和 RNA 输出序列。

> seq <- "ATCTCGGCGCGCATCGCGTACGCTACTAGC"
> revcomp(seq)
[1] "GCTAGTAGCGTACGCGATGCGCGCCGAGAT"

> seq <- c("TATAAT", "TTTCGC", "atgcat")
> revcomp(seq)
  TATAAT   TTTCGC   atgcat 
 "ATTATA" "GCGAAA" "atgcat"

查看manual 或TomKellyGenetics/tktools github 包存储库。

【讨论】：