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