Rosalind problema 10 – Finding a Shared Motif LCS

E la vamos nos em mais um problema de informática. Dado um n de sequências, achar a maior sub-sequência comum a todas as sequências.

Primeiro eu fui dar uma olhada no google e a primeira parada foi essa pagina da wikipedia, parece que esse problema não é especifico a bioinformática. Ele ensina um esquema muito louco usando matriz. Mas após pensar muito eu imaginei três caminhos a seguir. O esquema com matriz do wikipedia, o esquema com árvore que também é comentado no wikipedia e a boa e velha força bruta.

Mas vamos la.

Então o exemplo do problema fornece 3 sequências:

GATTACA TAGACCA ATACA

Aqui qual a maior sub-sequência comum a todo mundo? Nesse caso a maior so tem 2 pares de base. Mas existe mais de uma.

GAT--TA--CA TA--GACCA A--TA--CA

Todo mundo tem TA, mas outras sequências são comuns a todo mundo também como:

GATTA--CA TAGAC--CA CA ATA--CA GATT--AC--A TAG--AC--CA AC AT--AC--A

Certo, esse problema é fácil de entender, mas solucionar é meio complicado para mim pelo menos, o primeiro esquema com matriz é assim.

Vamos supor outras string so para ver uma sequência comum maior.

"ACACACACACGTAC" "GTGTACACACACAAA"

Então primeiro comparando so essas 2 string a gente faz o seguinte, usa uma como nomes das linhas e outra como nome das colunas de uma matriz.

Dai vamos comparar a primera letra da sequência 2 com a sequência 1, e vemos que o G da sequência 2 ocorre na 11 base da sequência 1

A C A C A C A C A C G T A C G 0 0 0 0 0 0 0 0 0 0 1 0 0 0

Dai vamos para a segunda linha.

A C A C A C A C A C G T A C G 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T 0 0 0 0 0 0 0 0 0 0 0 2 0 0

Agora um T da sequência 2, vemos que na 12 base ele ocorreu, mas opa na 11 ja tinha um 1, so olhar a diagonal ali, então eu preencho com um 2 ja que é a segunda ocorrência em linha. Ai preenchendo toda a matriz eu tenho:

A C A C A C A C A C G T A C G 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T 0 0 0 0 0 0 0 0 0 0 0 2 0 0 G 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T 0 0 0 0 0 0 0 0 0 0 0 2 0 0 A 1 0 1 0 1 0 1 0 1 0 0 0 3 0 C 0 2 0 2 0 2 0 2 0 2 0 0 0 4 A 1 0 3 0 3 0 3 0 3 0 0 0 1 0 C 0 2 0 4 0 4 0 4 0 4 0 0 0 2 A 1 0 3 0 5 0 5 0 5 0 0 0 1 0 C 0 2 0 4 0 6 0 6 0 6 0 0 0 2 A 1 0 3 0 5 0 7 0 7 0 0 0 1 0 C 0 2 0 4 0 6 0 8 0 8 0 0 0 2 A 1 0 3 0 5 0 7 0 9 0 0 0 1 0 A 1 0 1 0 1 0 1 0 1 0 0 0 1 0 A 1 0 1 0 1 0 1 0 1 0 0 0 1 0

Dai quanto eu tirar o maior numero dessa matriz, ele será a maior sequência comum, e a sequência será a diagonal, ou o nome das colunas da diagonal que ele ocorre. Mas como podem ocorrer várias sub-sequências, eu pego e retiro todas as subsequências por ordem em uma lista.

lista.cs [[1]] [1] "AC" "AC" "AC" "AC" "AC" "AC" "AC" "AC" "GT" "GT" "AC" "AC" "AC" [[2]] [1] "ACA" "ACA" "ACA" "ACA" "ACA" "ACA" "ACA" "GTA" [[3]] [1] "ACAC" "ACAC" "ACAC" "ACAC" "ACAC" "ACAC" "GTAC" [[4]] [1] "ACACA" "ACACA" "ACACA" "ACACA" "ACACA" [[5]] [1] "ACACAC" "ACACAC" "ACACAC" "ACACAC" [[6]] [1] "ACACACA" "ACACACA" "ACACACA" [[7]] [1] "ACACACAC" "ACACACAC" [[8]] [1] "ACACACACA"

E pimba, eu tenho de baixo para cima as maiores sequências comuns de acordo com o número de bases. Minha ideia original era dai usar essa lista, e procurar ela na outras sequências, sistematicamente, da maior para a menor. Mas como eu usei muito loop de for no R, sem vetorizar as operações, meu script ficou muitooooo lento. Para desenhar essa matriz com muitas bases demorava muito. O código funcionava perfeito, mas mais de uma 500 bases para cima ficava ruim. Eu pensei em algumas coisas para melhorar mais ai mudei de idéia.

Fui para a força bruta. O que fiz foi o seguinte, imaginemos uma sequência. A primeira do exemplo la em cima, se a gente aninhar 2 loops de for, da para retirar todas as sub-sequências possíveis.

> for(i in 1:nchar(sequencia)){ + for(j in 1:nchar(data[1])) { + print(substr(data[1],i,i+j)) + } + } [1] "GA" [1] "GAT" [1] "GATT" [1] "GATTA" [1] "GATTAC" [1] "GATTACA" [1] "GATTACA" … [1] "CA"

A partir daqui é só usar uma função ajudante para testar uma a uma dessas sequências se são sub-string das outras, no R com expressão regular, usando o grepl para um resultado de TRUE or FALSE, a partir dai é só fazer a sequência, a gente começa com um sub-string com 0 carácter, e sempre que algo for sub-string de todo mundo, e tiver mais caracteres que o sub-string atual, a gente sobrescreve ele, senão não, assim no final a gente vai ter um dos maiores sub-strings possíveis que todo mundo compartilha. Isso ficou um pouco mais rápido devido a minha incompetência para calcular a matriz ali em cima, mas no final fica a mesma coisa.

Então eu usei essa função no R, a is.substr é auxiliar, ela so testa se é sub-string de todo mundo.

#Função Auxiliar
is_substr<-function(find,data) {
if(length(data)<1 & nchar(find)<1) { T } else {
resp<-logical()
for(i in 1:length(data)) {
if(grepl(find,data[i])) {resp<-TRUE} else {resp<-FALSE;break}
}
resp
}
}
 
#Exemplo de como vai estar os dados
data<-c("GATTACA","TAGACCA","ATACA")
 
#Função para achar o menor sub-string comum
long_substr<-function(data) {
sub.m<-""
if(length(data) > 1 & nchar(data[1]) > 0) {
for(i in 1:nchar(data[1])){
for(j in 1:nchar(data[1])) {
if(is_substr(substr(data[1],i,i+j),data) &
nchar(sub.m)<nchar(substr(data[1],i,i+j))) {
sub.m<-substr(data[1],i,i+j)
}
}
}
}
sub.m
}

E ainda sim ficou meio lento, meio eu to sendo gentil, muito lento na verdade. Mas foi preguiça também, pois é possível comparar tudo de forma mais eficiente se compararmos do maior substring para o menor, e parar no primeiro que encontrar dai, mas ja estava triste com o primeiro código e parei.

E em python o código ficou idêntico ao do R. Que na verdade eu tinha visto esse código aqui em python, mas ele é bem simples, tanto que ficou a mesma coisa em R.

def is_substr(find, data):
    if len(data) < 1 and len(find) < 1:
        return False
    for i in range(len(data)):
        if find not in data[i]:
            return False
    return True
 
def long_substr(data):
    substr = ''
    if len(data) < 1 and len(data[0]) < 0:
        for i in range(len(data[0])):
            for j in range(len(data[0])-i+1):
                if j < len(substr) and is_substr(data[0][i:i+j], data):
                    substr = data[0][i:i+j]
    return substr
 
f = open('/home/augusto/Downloads/rosalind_lcs.txt', 'r')
lista = f.read().splitlines()
f.close()
 
print long_substr(lista)

So olhando que agora eu aprendi a abrir arquivos em python, pelo menos para tirar só as sequências. Eu abro com open, dai leio com read e separo em uma lista com splitlines e depois fecho o arquivo, que ele para de ocupar memoria. Mas a solução em python é bem mais rapida, ja que nele não precisa vectoriza as coisas.

E segue o código usado, esse problema foi complicado, principalmente por conta de fazer um método eficiente. No fórum de soluções as pessoas colocam o tempo que seus scripts levam, mas pelo menos deu para encontrar uma solução. Otimização de código ainda não esta em pauta por enquanto, segue todo o código que usei caso alguém queria dar uma olhada.

########
#Função para checar se algo é comum a todos os strings de um vetor
is_substr<-function(find,data) {
if(length(data)<1 & nchar(find)<1) { T } else {
resp<-logical()
for(i in 1:length(data)) {
if(grepl(find,data[i])) {resp<-TRUE} else {resp<-FALSE;break}
}
resp
}
}
 
#testando
is_substr("A",c("GATTACA","TAGACCA","ATACA"))
is_substr("G",c("GATTACA","TAGACCA","ATACA"))
 
#Agora fornecendo varios strings em um vetor
data<-c("GATTACA","TAGACCA","ATACA")
 
#função para encontrar o maior sub-string dos strings de um vetor
long_substr<-function(data) {
sub.m<-""
if(length(data) > 1 & nchar(data[1]) > 0) {
for(i in 1:nchar(data[1])){
for(j in 1:nchar(data[1])) {
if(is_substr(substr(data[1],i,i+j),data) &
nchar(sub.m)<nchar(substr(data[1],i,i+j))) {
sub.m<-substr(data[1],i,i+j)
}
}
}
}
sub.m
}
 
#Testando, mas fica bem lento conforme se aumenta
#o número e tamanho dos strings
long_substr(data)
 
#Abrindo os dados do Rosalind e transformando num Vetor
strings<-readLines("/home/augusto/Downloads/rosalind_lcs.txt")
strings<-as.vector(strings)
 
#Usando a função
long_substr(strings)
 
#Exemplo com matriz do wikipedia
string1<-"ACACACACACGTAC"
string2<-"GTGTACACACACAAA"
string3<-"ATACA"
 
#picando os 2 primeiros strings
str1<-strsplit(string1,"")[[1]]
str2<-strsplit(string2,"")[[1]]
 
#olhando o tamanho de cada um
m<-length(str1)
n<-length(str2)
 
#fazendo uma matriz em que os strings são os nomes da linhas e colunas
matriz<-matrix(0,ncol=m,nrow=n,dimnames=list(str2,str1))
matriz
 
#Preenchendo com 1 onde ocorre a mesma letra
for(i in 1:m) {
for(j in 1:n) {
if(str1[i]==str2[j]) {matriz[j,i]<-1}
}
}
 
matriz
 
#Colocando a sequencia, aqui precisa de melhoria para
#ficar mais rapido
menor<-min(c(m,n))
 
for(k in 1:(menor-1)) {
for(i in 1:(m-1)){
for(j in 1:(n-1)) {
if(matriz[j,i]==k & matriz[j+1,i+1]>=1) {matriz[j+1,i+1]<-k+1}
}
}
}
 
matriz
 
#Extraindo uma lista com os sub-strings
lcs.n<-max(matriz)
 
lista.cs<-vector("list",lcs.n-1)
lista.cs
 
for(j in (lcs.n-1):1) {
posi<-which(matriz==j+1,arr.ind=T)
for (i in 1:nrow(posi)) {
lista.cs[[j]][i]<-paste(str1[(posi[i,2]-(j+1)+1):posi[i,2]],collapse="")
}
}
 
#Conferindo com o resto dos strings, não cheguei a investir muito tempo
#nessa parte
for (i in 1:length(unlist(lista.cs))) {
if(grepl(rev(unlist(lista.cs))[i],string3)) {
cat(rev(unlist(lista.cs))[i],"presente","n")
break
} else {cat(rev(unlist(lista.cs))[i],"ausente","n")}
}

Leave a Reply

Your email address will not be published. Required fields are marked *