2016-07-15 3 views
2

단백질 (동일한 알파벳의 긴 시퀀스)에 정확한 펩타이드 (26 자리 알파벳의 짧은 문자열 인 A-Z)를 효율적으로 매핑하려고합니다. 이것을 수행하는 가장 효율적인 방법은 Aho-Corasick trie (펩타이드는 키워드 임)입니다. 불행히도 비 뉴클레오티드 알파벳 (Biostrings 'PDict 및 Starr의 match_ac은 모두 DNA에 대해 하드 코딩 됨)로 작동하는 R 버전의 AC를 찾을 수 없습니다.여러 문자열/키워드를 여러 텍스트와 효율적으로 일치 R

버팀목으로 나는 기본적인 grep 접근 방식을 병렬화하려고 노력했습니다. 하지만 상당한 IO 오버 헤드를 발생시키지 않으면 서 그렇게 할 수있는 방법을 찾는 데 어려움을 겪고 있습니다. 각 노동자에서 반환하는 것은 매우 자세한 그래서 ,

  • peptideInstances 여기 밀도 행렬 : 여기에 몇 가지 문제가 있습니다

    peptides = c("FSSSGGGGGGGR","GAHLQGGAK","GGSGGSYGGGGSGGGYGGGSGSR","IISNASCTTNCLAPLAK") 
    if (!exists("proteins")) 
    { 
        biocLite("biomaRt", ask=F, suppressUpdates=T, suppressAutoUpdate=T) 
        library(biomaRt) 
        ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") 
        proteins = getBM(attributes=c('peptide', 'refseq_peptide'), filters='refseq_peptide', values=c("NP_000217", "NP_001276675"), mart=ensembl) 
        row.names(proteins) = proteins$refseq_peptide 
    } 
    
    library(snowfall) 
    library(Biostrings) 
    library(plyr) 
    sfInit(parallel=T, cpus=detectCores()-1) 
    
    allPeptideInstances = NULL 
    i=1 
    increment=100 
    count=nrow(proteins) 
    while(T) 
    { 
        print(paste(i, min(count, i+increment), sep=":")) 
        text_source = proteins[i:min(count, i+increment),] 
        text = text_source$peptide 
    
        #peptideInstances = sapply(peptides, regexpr, text, fixed=T, useBytes=T) 
        peptideInstances = sfSapply(peptides, regexpr, text, fixed=T, useBytes=T) 
        dimnames(peptideInstances) = list(text_source$refseq_peptide, colnames(peptideInstances)) 
    
        sparsePeptideInstances = alply(peptideInstances, 2, .fun = function(x) {x[x > 0]}, .dims = T) 
    
        allPeptideInstances = c(allPeptideInstances, sparsePeptideInstances, recursive=T) 
        if (i==count | nrow(text_source) < increment) 
        break 
        i = i+increment 
    } 
    
    sfStop() 
    

    : 여기에 간단한 예입니다. 나는 블록을 으로 나누어서 40,000 (단백질) x 60,000 (펩타이드) 매트릭스를 다루지 않을 것입니다.

  • 단백질이 더 커지기 때문에 병렬 처리하는 것이 더 효과적 일 때 펩티드를 병렬 처리합니다. 그러나 나는 단백질에 의해 그것을하려고하는 것에 좌절감을 느꼈다 :
  • text_source에 단 하나의 단백질 만 있다면이 코드는 깨진다.

다른 방법으로 R의 더 나은 솔루션을 알고 있다면 기꺼이 사용하겠습니다. 나는 Aho-Corasick을 구현하는 데 더 도움이되었을 것입니다.

일부는 모호성 코드이지만 일부는 무시합니다.

답변

1

나는 Rcpp를 배우고 Aho-Corasick을 직접 구현했습니다. 이제 CRAN은 여러 가지 키워드 검색이 좋은 일반적인 목적을 가지고 있습니다. package.

listEquals = function(a, b) { is.null(unlist(a)) && is.null(unlist(b)) || !is.null(a) && !is.null(b) && all(unlist(a) == unlist(b)) } 

# simple search of multiple keywords in a single text 
keywords = c("Abra", "cadabra", "is", "the", "Magic", "Word") 
oneSearch = AhoCorasickSearch(keywords, "Is Abracadabra the Magic Word?") 
stopifnot(listEquals(oneSearch[[1]][[1]], list(keyword="Abra", offset=4))) 
stopifnot(listEquals(oneSearch[[1]][[2]], list(keyword="cadabra", offset=8))) 
stopifnot(listEquals(oneSearch[[1]][[3]], list(keyword="the", offset=16))) 
stopifnot(listEquals(oneSearch[[1]][[4]], list(keyword="Magic", offset=20))) 
stopifnot(listEquals(oneSearch[[1]][[5]], list(keyword="Word", offset=26))) 

# search a list of lists 
# * sublists are accessed by index 
# * texts are accessed by index 
# * non-matched texts are kept (to preserve index order) 
listSearch = AhoCorasickSearchList(keywords, list(c("What in", "the world"), c("is"), "secret about", "the Magic Word?")) 
stopifnot(listEquals(listSearch[[1]][[1]], list())) 
stopifnot(listEquals(listSearch[[1]][[2]][[1]], list(keyword="the", offset=1))) 
stopifnot(listEquals(listSearch[[2]][[1]][[1]], list(keyword="is", offset=1))) 
stopifnot(listEquals(listSearch[[3]], list())) 
stopifnot(listEquals(listSearch[[4]][[1]][[1]], list(keyword="the", offset=1))) 
stopifnot(listEquals(listSearch[[4]][[1]][[2]], list(keyword="Magic", offset=5))) 
stopifnot(listEquals(listSearch[[4]][[1]][[3]], list(keyword="Word", offset=11))) 

# named search of a list of lists 
# * sublists are accessed by name 
# * matched texts are accessed by name 
# * non-matched texts are dropped 
namedSearch = AhoCorasickSearchList(keywords, list(subject=c(phrase1="What in", phrase2="the world"), 
                verb=c(phrase1="is"), 
                predicate1=c(phrase1="secret about"), 
                predicate2=c(phrase1="the Magic Word?"))) 
stopifnot(listEquals(namedSearch$subject$phrase2[[1]], list(keyword="the", offset=1))) 
stopifnot(listEquals(namedSearch$verb$phrase1[[1]], list(keyword="is", offset=1))) 
stopifnot(listEquals(namedSearch$predicate1, list())) 
stopifnot(listEquals(namedSearch$predicate2$phrase1[[1]], list(keyword="the", offset=1))) 
stopifnot(listEquals(namedSearch$predicate2$phrase1[[2]], list(keyword="Magic", offset=5))) 
stopifnot(listEquals(namedSearch$predicate2$phrase1[[3]], list(keyword="Word", offset=11))) 

# named search of multiple texts in a single list with keyword grouping and aminoacid alphabet 
# * all matches to a keyword are accessed by name 
# * non-matched keywords are dropped 
proteins = c(protein1="PEPTIDEPEPTIDEDADADARARARARAKEKEKEKEPEPTIDE", 
      protein2="DERPADERPAPEWPEWPEEPEERAWRAWWARRAGTAGPEPTIDEKESEQUENCE") 
peptides = c("PEPTIDE", "DERPA", "SEQUENCE", "KEKE", "PEPPIE") 
peptideSearch = AhoCorasickSearch(peptides, proteins, alphabet="aminoacid", groupByKeyword=T) 
stopifnot(listEquals(peptideSearch$PEPTIDE, list(list(keyword="protein1", offset=1), 
               list(keyword="protein1", offset=8), 
               list(keyword="protein1", offset=37), 
               list(keyword="protein2", offset=38)))) 
stopifnot(listEquals(peptideSearch$DERPA, list(list(keyword="protein2", offset=1), 
               list(keyword="protein2", offset=6)))) 
stopifnot(listEquals(peptideSearch$SEQUENCE, list(list(keyword="protein2", offset=47)))) 
stopifnot(listEquals(peptideSearch$KEKE, list(list(keyword="protein1", offset=29), 
               list(keyword="protein1", offset=31), 
               list(keyword="protein1", offset=33)))) 
stopifnot(listEquals(peptideSearch$PEPPIE, NULL)) 

# grouping by keyword without text names: offsets are given without reference to the text 
names(proteins) = NULL 
peptideSearch = AhoCorasickSearch(peptides, proteins, groupByKeyword=T) 
stopifnot(listEquals(peptideSearch$PEPTIDE, list(1, 8, 37, 38))) 
stopifnot(listEquals(peptideSearch$DERPA, list(1, 6))) 
stopifnot(listEquals(peptideSearch$SEQUENCE, list(47))) 
stopifnot(listEquals(peptideSearch$KEKE, list(29, 31, 33))) 
+1

니스 :

는 여기에 몇 가지 사용 예입니다! 패키지를 사용하여 예제 (또는 두 개)를 추가 할 수 있습니다. – Jota