Строка соответствие с групповыми символами, пытаясь Biostrings пакет

Question

Учитывая строку patt:

patt = "AGCTTCATGAAGCTGAGTNGGACGCGATGATGCG" 

Мы можем сделать коллекцию коротких подстрок str_col:

str_col = substring(patt,1:(nchar(patt)-9),10:nchar(patt)) 

, которые мы хотим, чтобы соответствовать против subject1:

subject1 = "AGCTTCATGAAGCTGAGTGGGACGCGATGATGCGACTAGGGACCTTAGCAGC" 

обработка "N" в patt в качестве шаблона (соответствует любой букве в subject1), поэтому все подстроки в str_col соответствуют subject1.

Я хочу, чтобы это сопоставление строк в большой базе данных строк, и я нашел Bioconductor пакет Biostrings быть очень эффективным для этого. Но для того, чтобы быть эффективными, Biostrings требует, чтобы вы конвертировали свою коллекцию подстрок (здесь str_col) в словарь класса pdict, используя функцию PDidc(). Вы можете использовать этот «словарь» позже в таких функциях, как countPDict(), чтобы подсчитывать совпадения с целью.

Для использования подстановочных знаков вам необходимо разделить словарь на 3 части: голова (слева), доверенная группа (средний) и хвост (справа). У вас могут быть только подстановочные знаки, такие как «N», в голове или хвосте, но не в доверенной группе, и вы не можете иметь доверенный диапазон ширины = 0. Так, например, str_col[15] не будет соответствовать, если вы используете доверительная полоса с минимальной шириной = 1, как:

> PDict(str_col[1:15],tb.start=5,tb.end=5) 
Error in .Call2("ACtree2_build", tb, pp_exclude, base_codes, nodebuf_ptr, : 
    non base DNA letter found in Trusted Band for pattern 15 

потому что «N» находится справа в доверенной группе. Обратите внимание, что строки здесь являются последовательностями ДНК, поэтому «N» - это код для «соответствия A, C, G или T».

> PDict(str_col[1:14],tb.start=5,tb.end=5) #is OK 
TB_PDict object of length 14 and width 10 (preprocessing algo="ACtree2"): 
    - with a head of width 4 
    - with a Trusted Band of width 1 
    - with a tail of width 5 

Есть ли способ обойти это ограничение Biostrings? Я также попытался выполнить такую ​​задачу, используя базовые функции R, но я ничего не мог придумать.

Answer 1

Я считаю, что вам нужно будет сопоставлять некоторые другие дикие карты с IUPAC ambiguity code в какой-то момент, нет?

Если вам нужны идеальные соответствия и базовые функции вам достаточно, вы можете использовать тот же трюк, что и функция glob2rx(): просто используйте преобразование с gsub(), чтобы построить соответствующие шаблоны. Пример:

IUPACtoRX <- function(x){ 
    p <- gsub("N","\\[ATCG\\]",x) 
    p <- gsub("Y","\\[CT\\]",p) #match any pyrimidine 
    # add the ambiguity codes you want here 
    p 
} 

Очевидно, что вам нужна строка для каждой неоднозначности вы хотите запрограммировать в, но это довольно просто, я бы сказал.

Делая это, вы можете, например, сделать что-то вроде:

> sapply(str_col, function(i) grepl(IUPACtoRX(i),subject1)) 
AGCTTCATGA GCTTCATGAA CTTCATGAAG TTCATGAAGC TCATGAAGCT CATGAAGCTG ATGAAGCTGA 
     TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
TGAAGCTGAG GAAGCTGAGT AAGCTGAGTN AGCTGAGTNG GCTGAGTNGG CTGAGTNGGA TGAGTNGGAC 
     TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
GAGTNGGACG AGTNGGACGC GTNGGACGCG TNGGACGCGA NGGACGCGAT GGACGCGATG GACGCGATGA 
     TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
ACGCGATGAT CGCGATGATG GCGATGATGC CGATGATGCG 
     TRUE  TRUE  TRUE  TRUE 

Чтобы найти число совпадений, вы можете использовать, например, gregexpr():

> sapply(str_col,function(i) sum(gregexpr(IUPACtoRX(i),subject1) > 0)) 
AGCTTCATGA GCTTCATGAA CTTCATGAAG TTCATGAAGC TCATGAAGCT CATGAAGCTG ATGAAGCTGA 
     1   1   1   1   1   1   1 
TGAAGCTGAG GAAGCTGAGT AAGCTGAGTN AGCTGAGTNG GCTGAGTNGG CTGAGTNGGA TGAGTNGGAC 
     1   1   1   1   1   1   1 
GAGTNGGACG AGTNGGACGC GTNGGACGCG TNGGACGCGA NGGACGCGAT GGACGCGATG GACGCGATGA 
     1   1   1   1   1   1   1 
ACGCGATGAT CGCGATGATG GCGATGATGC CGATGATGCG 
     1   1   1   1