Procesando resultados de BLAST

Con el archivo FASTA generado en el ejercicio anterior 1lfu.fas hacemos dos búsquedas BLAST (Altschul et al., 1997) con los siguientes comandos (necesitas tener instalado BLAST localmente 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 en Julio de 2005

# Programa extraeBLAST.pl que lee un archivo generado por BLAST pasado como argumento 
# y extrae la información relevante de cada alineamiento reportado

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 está leyendo Q o S
my ($firstQ,$lastQ,$seqQ);     # primer resíduo, último resíduo 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 y asígnalo a  
}


## abre $archivoBLAST en modo lectura y recórrelo línea a línea ####################

$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 guárdalo 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 (sólo 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 guardó así: [($firstQ,$lastQ,$firstS,$lastS,$name,$id,$ev,$seqQ,$seqS)]
# el primer alineamiento será $hits[0], que a su vez es un arreglo, por eso necesitamos otra
# coordenada para llegar a cada elemento

if($hits[0])
{
	print "Posible homólogo:\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.