Дополнить последовательность ДНК

Предположим, у меня есть последовательность ДНК. Я хочу получить дополнение этого. Я использовал следующий код, но я не получаю его. Что я делаю неправильно?

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"
}

6 ответов

Использование chartr который построен для этой цели:

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

Просто дайте ему две строки символов одинаковой длины и свою строку. Также векторизация аргумента для перевода:

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

Заметьте, я делаю замену на строковом представлении ДНК, а не вектора. Чтобы преобразовать вектор, я бы создал карту поиска как именованный вектор и индекс, который:

> 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"

Пакет http://bioconductor.org/ Biostrings имеет много полезных функций для такого рода операций. Установите один раз:

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

затем используйте

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

Для дополнения, как в верхнем, так и в нижнем регистре, вы можете использовать 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"

Существует также пакет seqinr

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

Biostrings имеет гораздо меньший профиль памяти, но seqinr хорош также потому, что вы можете выбрать регистр оснований (включая смешанные) и изменить их на что угодно, например, если вы хотите смешать T и U в одной и той же последовательности. Biostrings заставляет вас иметь или T или U.

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"
> 

Здесь ответ с использованием базы р. Написано с ужасным форматированием, чтобы прояснить ситуацию и сохранить его как одну строчку. Поддерживает верхний и нижний регистры.

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

Я обобщил решение rev(comp(seq)) с seqinr пакет:

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

Эта версия совместима со строковыми входами и векторизована для обработки списка или векторного ввода нескольких строк. Выходной класс должен соответствовать входным данным, включая регистры и типы. Это также поддерживает входные данные, содержащие "U" для выходных последовательностей РНК и РНК.

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

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

См. Руководство или репозиторий пакетов github https://github.com/TomKellyGenetics/tktools.

Другие вопросы по тегам