Procesando resultados de BLAST

Con el archivo FASTA generado en el ejercicio anterior 1lfu.fas hacemos dos búsquedas BLAST con los siguientes comandos (necesitas tener instalado BLASTlocalmente y bajarte la biblioteca de secuencias nr):

blastall -p blastp -i 1lfu.fas -d nr -o 1lfu.blast -v 20 -b 20           # -o es el archivo de resultados
blastpgp -i 1lfu.fas -j 2 -Q 1lfu.pssm -d nr -o 1lfu.blast -v 20 -b 20   # -o y -Q son archivos de resultados

Los archivos de resultados generados son 1lfu.blast y 1lfu.blastpgp, donde se encuentra la información relativa a las secuencias potencialmente homólogas encontradas en la biblioteca nr ( non redundant ), que puedes obtener en el NCBI.

El último archivo de resultados generado es 1lfu.pssm, que contiene información también muy valiosa, entre la que destacamos una matriz de sustitución de aminoácidos específica para la familia de proteínas 1lfu y el contenido en información evolutiva de cada aminoácido de la secuencia de entrada 1lfu.fas.

En el siguiente ejemplo os muestro un trozo de código que extrae de la salida de BLAST (indistintamente 1lfu.blast o 1lfu.blastpgp) la información relevante de cada uno de los alineamientos reportados.

#!/usr/bin/perl -w 
# Ejemplo escrito por Bruno Contreras

# Programa extraeBLAST.pl que lee un archivo generado por BLAST pasado como argumento 
# y extrae la informacion relevante de cada alineamiento reportado
# Ejemplo de como usar variables bandera (flags) para leer textos complejos

# Ojo: se puede hacer con Bioperl:
# Ojo2: si ejecutas blast con -outfmt 6 es trivial leer la salida de BLAST, si no te importa perder el alineamiento

use strict; 

## variables ######################################################################

my ($archivoBLAST,@hits);

# otras variables manejadas, donde Q significa 'query' y S ' subject',
# las dos secuencias de un alineamiento en BLAST
my ($aux0,$aux1,$aux2,$aux3);  # variables auxiliares
my ($Q,$S);                    # variables booleanas que indican si se esta leyendo Q o S
my ($firstQ,$lastQ,$seqQ);     # primer residuo, ultimo residuo y secuencia de 'query'
my ($firstS,$lastS,$seqS);
my ($longQ,$name,$id,$ev);     # longitud de 'query', nombre se secuencia, % de identidad y e-value


## revisa los argumentos pasados al invocar el programa ############################
if(scalar(@ARGV) < 1)
{	
	die "uso: extraeBLAST.pl <archivoBLAST>\n";	
}
else 
{	
	$archivoBLAST = $ARGV[0];	# lee el nombre del argumento
}


## abre $archivoBLAST en modo lectura y recorrelo linea a linea

$Q = $S = 0;
open(BLAST,$archivoBLAST) || die "# no puedo abrir $archivoBLAST\n";
while(<BLAST>)
{
	if(/Query:/)
	{
		($aux0, $aux1, $aux2, $aux3) = split;  # split separa $_ en 4 escalares separados por espacios en blanco
		
		if(!$Q){	$firstQ = $aux1;	}
		
		$lastQ = $aux3;
		$seqQ .= $aux2;
		$Q++;  # estamos leyendo 'query'        
	}
	elsif(/Sbjct:/)
	{
		($aux0, $aux1, $aux2, $aux3) = split;
	
		if(!$S){	$firstS = $aux1;	}

		$lastS = $aux3;
		$seqS .= $aux2;
		$S++;  # estamos leyendo 'subject'
	}  
	elsif((/^>/)||(/^  Database:/))	## procesa el alineamiento previo
	{ 
		if($Q)
		{
			# crea un arreglo con los datos de este alineamiento y guardalo en @hits
			# @hits es por tanto un arreglo de arreglos
			push(@hits, [($firstQ,$lastQ,$firstS,$lastS,$name,$id,$ev,$seqQ,$seqS)] );
		}
	
		# reiniciamos valores para el siguiente alineamiento
		$Q = $S = 0;
		$seqQ = $seqS = $name = $id = $ev = "";
		
		# nombre del siguiente alineamiento
		$name = $_;
		chomp($name);
	}
	elsif(/^ Score =/)	
	{	
		if($Q) ## secuencias que tienen varios alineamientos 
		{		
			push(@hits, [($firstQ,$lastQ,$firstS,$lastS,$name,$id,$ev,$seqQ,$seqS)] );
					
			$Q = $S = 0;
			$seqQ = $seqS = $id = ""; 
		}
			
		$ev = (split)[7];
	}
	elsif(/^ Identities/)
	{	
		# toma nota del % de identidad del alineamiento
		$id = (split)[3]; 
		$id =~ s/(\(|\)|\%|,)//g;	
	}
	elsif(/^Sequences/)	# reinicia todo para la siguiente ronda de resultados (solo blastpgp)
	{
		undef @hits;
		$S = $Q = 0;
		$seqQ = $seqS = $name = $id = $ev = "";
	}
	elsif(/letters\)/)
	{	
		# toma nota de la longitud de la secuencia 'query'
		$longQ = (split)[0];	
		$longQ =~ s/\(//;	
	}
}
close(BLAST);

## ahora como ejemplo muestra el primer alineamiento guardado #################################

# cada alineamiento se guarda asi: [($firstQ,$lastQ,$firstS,$lastS,$name,$id,$ev,$seqQ,$seqS)]
# el primer alineamiento es $hits[0], que a su vez es un arreglo, por eso necesitamos otra
# coordenada para llegar a cada elemento

if($hits[0])
{
	print "Posible homologo:\n$hits[0][4]\n";          
	print "Fragmento de $hits[0][0] a $hits[0][1], e-value $hits[0][6] y \%ID $hits[0][5]\n";    
	print "Alineamiento\n$hits[0][7]\n$hits[0][8]\n";
}

Vuestra tarea es ahora escribir un programita que lea 1lfu.pssm y genere una lista de los aminoácidos cuyo contenido en información evolutiva sea elevado. Os sugiero que uséis la función split para leer el archivo línea a línea.