Более быстрый способ разделить строку и символы количества с помощью R?
Я ищу более быстрый способ вычислить содержание GC для строкового чтения DNA в из файла FASTA. Это сводится к взятию строки и подсчету количества раз, что буква 'G' или 'C' появляется. Я также хочу указать диапазон символов для рассмотрения.
У меня есть рабочая функция, которая является довольно медленной, и она вызывает узкое место в моем коде. Это похоже на это:
##
## count the number of GCs in the characters between start and stop
##
gcCount <- function(line, st, sp){
chars = strsplit(as.character(line),"")[[1]]
numGC = 0
for(j in st:sp){
##nested ifs faster than an OR (|) construction
if(chars[[j]] == "g"){
numGC <- numGC + 1
}else if(chars[[j]] == "G"){
numGC <- numGC + 1
}else if(chars[[j]] == "c"){
numGC <- numGC + 1
}else if(chars[[j]] == "C"){
numGC <- numGC + 1
}
}
return(numGC)
}
Выполнение Rprof дает мне следующий вывод:
> a = "GCCCAAAATTTTCCGGatttaagcagacataaattcgagg"
> Rprof(filename="Rprof.out")
> for(i in 1:500000){gcCount(a,1,40)};
> Rprof(NULL)
> summaryRprof(filename="Rprof.out")
self.time self.pct total.time total.pct
"gcCount" 77.36 76.8 100.74 100.0
"==" 18.30 18.2 18.30 18.2
"strsplit" 3.58 3.6 3.64 3.6
"+" 1.14 1.1 1.14 1.1
":" 0.30 0.3 0.30 0.3
"as.logical" 0.04 0.0 0.04 0.0
"as.character" 0.02 0.0 0.02 0.0
$by.total
total.time total.pct self.time self.pct
"gcCount" 100.74 100.0 77.36 76.8
"==" 18.30 18.2 18.30 18.2
"strsplit" 3.64 3.6 3.58 3.6
"+" 1.14 1.1 1.14 1.1
":" 0.30 0.3 0.30 0.3
"as.logical" 0.04 0.0 0.04 0.0
"as.character" 0.02 0.0 0.02 0.0
$sampling.time
[1] 100.74
Совет для того, чтобы сделать этот код быстрее?
4 ответа
Лучше вообще не разделять, просто посчитайте совпадения:
gcCount2 <- function(line, st, sp){
sum(gregexpr('[GCgc]', substr(line, st, sp))[[1]] > 0)
}
Это на порядок быстрее.
Небольшая функция C, которая просто выполняет итерацию по символам, будет еще на порядок быстрее.
Здесь нет необходимости использовать цикл.
Попробуйте следующее:
gcCount <- function(line, st, sp){
chars = strsplit(as.character(line),"")[[1]][st:sp]
length(which(tolower(chars) == "g" | tolower(chars) == "c"))
}
Не знаю, быстрее ли это, но вы можете посмотреть пакет R seqinR - http://pbil.univ-lyon1.fr/software/seqinr/home.php?lang= eng . Это отличный общий пакет биоинформатики со множеством методов анализа последовательностей. Это в CRAN (который, кажется, не работает, когда я пишу это).
Содержимое сборщика мусора будет выглядеть следующим образом:
mysequence <- s2c("agtctggggggccccttttaagtagatagatagctagtcgta")
GC(mysequence) # 0.4761905
Это строка, которую вы также можете прочитать в файле fasta, используя "read.fasta ()".