Automatització

Per tal d'automatitzar els programes que havíem de fer servir, hem creat uns sèrie de programes en llenguatge bash i en llenguatge perl per tal de facilitar el processament d'un gran volum de dades. Concretament hem creat 5 programes bash i 3 programes perl:

Per tal que algun usuari desitgi fer-los servir hem escrit un tutorial bàsic per a fer-lo servir:

En primer lloc, hem de crear una carpeta al directori dels programes anomenada Llibreria, a on s'han de colocar cada seqüència que s'utilitzarà com a query en un arxiu diferent en format fasta (ex: creatina.fa). Cal col·locar els genomes que volem fer servir en un directori tots junts de la següent manera. S'ha de posar dins d'una carpeta amb el nom de l'organisme el genoma formatejat(link) i indexat(link) de manera que el nom de la base de dades sigui "genoma.fa". Si la direcció no és la mateixa que hem fet servir nosaltres, s'ha de modificar el programa i escriure-la on es defineix la variable "adresa". També cal tenir instal·lats els programes blast, exonerate, genewise i t_coffee.

1. run1_blast.bash: Aquest programa t'executarà un blast de cada proteïna amb cada genoma que estiguin en les carpetes de Llibreria i en la ubicació de genomes.

Requereix:

Per visualitzar i valorar el resultat obtingut es crearà un arxiu anomenat "resultats_blast.txt" on es veuran només els hits que siguin significatius.

Codi:

#!/bin/bash

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu

### Aqui definim quina és la ruta pels genomes

adresa="/cursos/BI/genomes/protists/2011"

### Això és per a detectar quines proteines hi ha a la carpeta Llibreria

for proteines_fa in ./Llibreria/*; do {
basename_fa=`basename $proteines_fa`
proteina=${basename_fa%.fa}
proteines=$proteines" $proteina"
} done

### Aquí comença el programa pròpiament dit: dos bucles for per a fer les combinaciones entre proteïnes i genomes

rm resultats_blast.txt 2> error.temp

	for protein in $proteines ; do {
		for genomes in P.tricornutum B.hominis H.arabidopsidis P.ultimum S.parasitica E.siliculosus C.muris C.parvum A.taiwanensis C.owczarzaki P.pallidum C.merolae T.trahens N.gruberi; do {

			patroalfa=`cat ./Llibreria/$protein.fa | egrep ^">"`
			patrobeta=${patroalfa#*>}
			echo $patrobeta > temppatro
			patro=`cut -f 1 -d " " temppatro`
			rm temppatro
			
			blastall -p tblastn -i ./Llibreria/$protein.fa -d /cursos/BI/genomes/protists/2011/$genomes/genome.fa -o output_blast/$protein.$genomes.txt -m 9
			echo ">${protein}_$genomes:" >> resultats_blast.txt
			egrep $patro output_blast/$protein.$genomes.txt | egrep e- >> resultats_blast.txt
			
			x=`egrep $patro ./output_blast/$protein.$genomes.txt | egrep e-`
			if [ "$x" != "" ]; then
				egrep $patro output_blast/$protein.$genomes.txt | egrep e- > arxiu.temp
				mkdir ./blast_selec/$protein 2> error.temp
				mkdir ./blast_selec/$protein/$genomes 2> error.temp
				./dosificador.pl $protein $genomes < arxiu.temp

			fi
		
		
			echo ${protein}_$genomes fet
		} done
	} done

rm ./*.temp

Codi "dosificador.pl":

#!/usr/bin/perl

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu

use strict;
my $contador=0;

while ()
{
	open (fitxersortida, ">./blast_selec/$ARGV[0]/$ARGV[1]/hit$contador.txt");
	print fitxersortida $_;
	close fitxersortida;
	$contador++;
}

2. blast_extend.bash (opcional): Aquest programa t'executarà un blast de cada proteïna amb cada genoma que estiguin en les carpetes de Llibreria i en la ubicació de genomes i et mostrarà el resultat de manera extensa (inclou els alineaments).

Requereix:

Codi:

#!/bin/bash

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu


### Aqui definim quina és la ruta pels genomes 

adresa="/cursos/BI/genomes/protists/2011"

### Això és per a detectar quines proteines hi ha a la carpeta Llibreria

for proteines_fa in ./Llibreria/*; do {
basename_fa=`basename $proteines_fa`
proteina=${basename_fa%.fa}
proteines=$proteines" $proteina"
} done

### Aquí comença el programa pròpiament dit: dos bucles for per a fer les combinaciones entre proteïnes i genomes

	for protein in $proteines ; do {
		for genomes in P.tricornutum H.arabidopsidis P.ultimum S.parasitica E.siliculosus C.muris C.parvum A.taiwanensis C.owczarzaki P.pallidum C.merolae T.trahens N.gruberi; do {

			blastall -p tblastn -i ./Llibreria/$protein.fa -d /cursos/BI/genomes/protists/2011/$genomes/genome.fa -o output_blast_extend/$protein.$genomes.txt
			echo ${protein}_$genomes fet
		} done
	} done

3. run2_exonerate.bash: Aquest programa executarà l'exonerate per als hits obtinguts amb el programa "run1_blast.bash".

Requereix:

Per saber quins exonerates han estat capaços de predir gens es crea un arxiu anomenat "resultats_exonerate.txt" on només es mostren els exonerates que han predit gens.

Codi:

#!/bin/bash

rm resultats_exonerate.txt 2> error.temp

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu

adresa="/cursos/BI/genomes/protists/2011"

### Això és per poder triar si volem exonerate exhaustiu o no

if [ "$1" = "-exhaustive" ]; then
exhaustive="--exhaustive yes"
fi

for proteins in blast_selec/* ; do {
       for genome in $proteins/* ; do {
               for hit in $genome/*; do {

		### Assignem a les variables els noms sense rutes ni extensions (sense "/exemple/" ni ".etc")
		
		proteins=`basename $proteins`
        	genome=`basename $genome`
        	hit=`basename $hit`
		hit=${hit%.txt}

		echo ${proteins}_${genome}_${hit}:

		### Això és per a detectar l'identificador que apareix al fitxer fasta a la primera línea

		patroalfa=`cat ./blast_selec/$proteins/$genome/$hit.txt | egrep ^">"`
		patrobeta=${patroalfa#*>}
		echo $patrobeta > temppatro
		patro=`cut -f 1 -d " " temppatro`
		rm temppatro       

		### Per emagatzemar la regió on apareix el hit

		regio=`cat $patro ./blast_selec/$proteins/$genome/$hit.txt | cut -f 2`
		fastafetch $adresa/$genome/genome.fa ./genomes/$genome/genome.index $regio > ./regions/$regio.fa

		### Aquí definim l'inici i el final del hit en la regió

		inici=`cat $patro ./blast_selec/$proteins/$genome/$hit.txt | cut -f 9`
		final=`cat $patro ./blast_selec/$proteins/$genome/$hit.txt | cut -f 10`

		if [ "$inici" -gt "$final" ]; then
			inici="$final"
		fi
		
		start=$(( $inici - 20000 ))
		length="50000"

		### FASTASUBSEQ : Retalla una regió al voltant del gen que volem fer servir		

		fastasubseq ./regions/$regio.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg.temp

		### Això ens servirà per modificar variables en cas que no hi hagi prou bases en la regio o que l'inici estigui massa aprop
					
		error=`cat error_msg.temp`
		if [ "$error" != "" ]; then
			error_llarg=`cat error_msg.temp | egrep "before end"`
			error_curt=`cat error_msg.temp | egrep "Unknown flag"`
						
			### Si es passa de llarg, establim la llargada perquè arribi al final
			if [ "$error_llarg" != "" ]; then
				llargada=$error_llarg						
				llargada=${llargada##*(}
				llargada=${llargada%)}
				length=$(( $llargada - $start))
				fastasubseq ./regions/$regio.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa
			fi
						
			### Si l'inici està massa aprop establim l'start en 0
			if [ "$error_curt" != "" ]; then
				start="0"
				fastasubseq ./regions/$regio.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg_2.temp
				error_2=`cat error_msg_2.temp`
				if [ "$error_2" != "" ]; then
					error_2_llarg=`cat error_msg_2.temp | egrep "before end"`
					if [ "$error_2_llarg" != "" ]; then
						llargada=$error_2_llarg						
						llargada=${llargada##*(}
						llargada=${llargada%)}
						length=$(( $llargada - $start))
						fastasubseq ./regions/$regio.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg_3.temp
					fi
				fi 
			fi
		fi

		
		### EXONERATE : Predicció de gens
		
		echo "	Iniciant exonerate:"

		exonerate -m p2g --showtargetgff -q ./Llibreria/$proteins.fa -t ./subseq/${proteins}_${genome}_${hit}_genomic.fa $exhaustive > ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff 2> error_de_substitucio.temp
	
		### Per si al fasta de sequencia hi ha Us, cal substituir-les per X, però només un cop!				

		error_substitucio=`cat error_de_substitucio.temp | egrep "Unknown amino acid"`
		if [ "$error_substitucio" = "** FATAL ERROR **: Unknown amino acid [U]" ]; then
			if [ "$int_UX" = "0" ]; then
				./substitucio.pl < ./Llibreria/$proteins.fa > ${protein}_temporal.txt
				int_UX="1"
			fi
			exonerate -m p2g --showtargetgff -q ${proteins}_temporal.txt -t ./subseq/${proteins}_${genome}_${hit}_genomic.fa $exhaustive > ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff 2> error_de_substitucio2.temp				
		fi
		gene=`grep gene ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff`
		if [ "$gene" != "" ]; then
			echo ${proteins}_${genome}_${hit} : >> resultats_exonerate.txt
		fi

		echo "	Exonerate acabat"	

		### Per extreure el cDNA i traduir-lo a proteïnes				
		
		echo "	Preparant i executant per a t_coffee:"
				
		exon=`grep -w exon ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff`
		grep -w exon ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff > exon.temp
		if [ "$exon" != "" ]; then
			fastaseqfromGFF.pl ./subseq/${proteins}_${genome}_${hit}_genomic.fa exon.temp > ./cDNA/cDNA_${proteins}_${genome}_${hit}.fa
			fastatranslate ./cDNA/cDNA_${proteins}_${genome}_${hit}.fa > ./translations/aa_${proteins}_${genome}_${hit}.fa
		fi
	
		echo ${proteins}_${genome}_${hit} fet
		
		} done
	} done
} done

rm *.temp

Codi "substitucio.pl":

#!/usr/bin/perl

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu

use strict;

my @v;
my $string;
my $counter = 0;

while ()
{
	if ($counter == 0)
	{
		print $_;
	}
	if ($counter > 0)
	{

		@v=split(//,$_);
	
		while (<@v>)
		{
			if ($_ eq "U")
			{		
				$string.= "X";
			}		
			else
			{
				$string.=$_;
			}
		}
	}
	$counter++;
	$string.="\n";
}

print $string;

4. run3_genewise.bash: Aquest programa executarà el genewise per als hits obtinguts amb el programa "run1_blast.bash".

Requereix:

Codi:

#!/bin/bash

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu

for proteins in blast_selec/* ; do {
       for genome in $proteins/* ; do {
               for hit in $genome/*; do {
		
		proteins=`basename $proteins`
        	genome=`basename $genome`
        	hit=`basename $hit`
		hit=${hit%.txt}

		genewise -pep -pretty -cdna -gff -both Llibreria/$proteins.fa ./subseq/${proteins}_${genome}_${hit}_genomic.fa > genewise/genewise_${proteins}_${genome}_${hit}.txt		

		./getpepgene.pl < genewise/genewise_${proteins}_${genome}_${hit}.txt > genewise/PEP_g_${proteins}_${genome}_${hit}.fa

		echo ${proteins}_${genome}_${hit} fet

		} done
	} done
} done

Codi "getpepgene.pl":

#!/usr/bin/perl -w

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu

use strict;

my $interruptor=0;
my $contador=0;
my @v;
my $cont1=0;
my $cont2=0;
my $string1;
my $string2;
my $int2=0;



while ()
{
	if ($_ =~ /\//)
	{
		if ($interruptor == 1)
		{
			$int2=1;
		}
		$interruptor = 0;
	}
	while ($interruptor == 1 && $contador == 0)
	{
		
		if ($int2 == 0 )
		{
			$string1.= $_;
			@v=split(//,$_);
			$cont1+=scalar(@v);
		}
		if ($int2 == 1 )
		{
			$string2.= $_;
			@v=split(//,$_);
			$cont2+=scalar(@v);
		}
		$contador=1;
	}
	if ($_ =~ /.pep/)
	{
		$interruptor=1;
		if ($int2 == 0 )
		{
			$string1.=$_;
		}
		if ($int2 == 1 )
		{
			$string2.=$_;
		}
	}
	$contador=0
}

if ($cont1 > $cont2)
{
	print $string1;
}
else
{
	print $string2;
}

5. run4_t_coffee.bash: Aquest programa executarà un alineament amb T-Coffee amb els resultats de l'exonerate i el genewise amb la seva respectiva query.

Requereix:

Codi:

#!/bin/bash

# Programa fet pel grup By de 4t de biologia del curs 2010-2011 de la Pompeu Fabra
# Contacte: marcel.costa01@estudiant.upf.edu

for file in exonerate_output/PEP_e_* ; do {
name=`basename $file`
name=${name%.fa}
name=${name#PEP_e_}

protein=${name%%_*}
echo $name
echo $protein
t_coffee Llibreria/$protein.fa $file > ./t_coffee_e/t_coffee_e_$name.txt 2> error.temp
mv $protein.aln ./t_coffee_e/t_coffee_e_$name.aln
mv $protein.html ./t_coffee_e/t_coffee_e_$name.html
echo $file done
} done

for file in genewise/PEP_g_* ; do {
name=`basename $file`
name=${name%.fa}
name=${name#PEP_g_}

protein=${name%%_*}
echo $name
echo $protein
t_coffee Llibreria/$protein.fa $file > ./t_coffee_g/t_coffee_g_$name.txt 2> error.temp
mv $protein.aln ./t_coffee_g/t_coffee_g_$name.aln
mv $protein.html ./t_coffee_g/t_coffee_g_$name.html
echo $file done
} done

rm *.temp