2014-09-02 38 views
1

給出字符串patt字符串使用通配符匹配,試圖Biostrings包

patt = "AGCTTCATGAAGCTGAGTNGGACGCGATGATGCG" 

我們可以更短的子str_col的集合:

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

我們要匹配一個subject1

subject1 = "AGCTTCATGAAGCTGAGTGGGACGCGATGATGCGACTAGGGACCTTAGCAGC" 

在中處理「N」作爲通配符(匹配subject1中的任何字母),因此str_col中的所有子字符串都與subject1匹配。

我想在字符串的大型數據庫中做這種字符串匹配,我發現BiostringsBiostrings是非常有效的做到這一點。但是,爲了高效,Biostrings要求您使用函數PDidc()將您的子串集合(此處爲str_col)轉換爲類pdict的字典。您可以稍後在countPDict()等函數中使用此'詞典'來計算與目標的匹配。

爲了使用通配符,您必須將字典分爲三部分:頭(左),可信帶(中)和尾(右)。您只能在頭部或尾部使用通配符,例如「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」在可信頻段中是正確的。注意這裏的字符串是DNA序列,所以「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基函數來執行這樣的任務,但我無法想出任何東西。

回答

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 
+0

不錯!非常感謝! – vitor 2014-09-02 14:14:12

+0

不確定如果您可以使用長度(gregexpr()),因爲它也會返回1以表示不匹配。例如。如果你省略函數IUPACtoRX,它也會爲每個人返回1。 – vitor 2014-09-02 14:49:16

+0

@agui不錯,我忘記了gregexpr在不匹配的情況下返回-1。我已經改編 – 2014-09-02 16:16:48