Materials i mètodes



A continuació tenim una explicació detallada del procediment que hem seguit per estudiar la presència de selenoproteïnes en genomes de protistes seqüenciats recentment. El format dels paths de les comandes que trobareu al llarg de la pàgina, està expressat de la mateixa manera que en el script d'automització que podeu trobar a l'apartat corresponent. La pipeline que hem utilitzat és semblant a la de Selenoprofiles.[9][10][11][12][13][14]



1. Obtenció de les queries

Hem estudiat la presència de la seqüència proteica de les selenoproteïnes GPx i FrnE en diferents organismes protistes.



Per obtenir la seqüència proteica de GPx hem utilitzat la base de dades SelenoDB. Hem obtingut la seqüència d'aquesta proteïna d'H. sapiens perquè aquest és l'organisme on s'han descrit més selenoproteïnes de la família Glutathione Peroxidase (GP), tot i que el criteri teòric generlament emprat per comparar seqüències és usar un organisme que sigui filogenèticament proper.



MCAARLAAAAAAAQSVYAFSARPLAGGEPVSLGSLRGKVLLIENVASLXGTTVRDYTQMN ELQRRLGPRGLVVLGFPCNQFGHQENAKNEEILNSLKYVRPGGGFEPNFMLFEKCEVNGA GAHPLFAFLREALPAPSDDATALMTDPKLITWSPVCRNDVAWNFEKFLVGPDGVPLRRYS RRFQTIDIEPDIEALLSQGPSCA






La seqüència de la proteïna FrnE ens l'han proporcionat els professors, i correspon a Ciona intestinalis, un urocordat.



MSKFLHIRVVSDIMXPWCWVGKRNLETAMKASENEYKFKVTWEPFLLRPNMPLEGIAKQG DFGPHTPGAQRLINVGRKVGVEFAFKQPRYPYTLFGHCALEFAIKKDPSGSIQTQLQESLFK SYFTDGEYPDVETVSTVAATCGLNREEVKSFISDESNLAAVKRKAAQWSANGVSGVPYFII NDCPVFSGAQEPAAFQNIFAKVAEKYPMPSQ





2. Obtenció de genomes de protistes

Els genomes a investigar, seqüenciats recentment, van ser facilitats pels professors de l'assignatura de bioinformàtica i pertanyen a les següents espècies: A. albachii, A. rara, C. fasciculata, D. discoideum, D. fasciculatum, F. cylindrus, G. niphandrodes, I. multifiliis, L. donovani, L. tarentolae, P. capsici, P. polycephalum, S. artica, T. congolense.





3. Cerca de similaritat: BLAST i tBLASTn

El BLAST (Basic Local Aligment Search Tool) és un programa que utilitza un algorisme heurístic per tal d'obtenir alineaments locals de seqüències. Hi ha diferents tipus de BLAST, nosaltres hem utilitzat tBLASTn, que permet comparar seqüències proteiques (query) contra bases de dades nucleotídiques. Per fer-ho, el programa tradueix les seqüències nucleotídiques a seqüències proteiques en els sis possibles marcs de lectura (tres forward i tres reverse) i els alinea amb la nostra query. [10] La comanda que hem utilitzar per realitzar el BLAST és:



blastall -p tblastn -i ./query_proteins/$protein.fa -d $genome_path/$genomes/genome.fa -o ./output_blast_detailed/$protein.$genomes.txt -F F




blastall -p tblastn -i ./query_proteins/$protein.fa -d $genome_path/$genomes/genome.fa -o ./output_blast/$protein.$genomes.txt -m 9 -F F




on,

-p: indica el tipus de BLAST que estem utilitzat (en el nostre cas tblastn).

-i: indica quina query estem utilitzant.

-d: indica contra quin genoma o base de dades de genomes comparem la nostra query.

-o: indica l'output.

-F F: serveix per eliminar les regions de baixa complexitat i permet augmentar la qualitat de l'alineament.

-m 9: que permet visualitzar la informació d'una forma alternativa.





4. Anàlisi del BLAST

Hem acceptat com a vàlids hits amb un e-value menor a 0.0001, que indica el nombre de hits igual o més bons que l'obtingut que esperariem trobar per atzar a la nostra base de dades.





5. Extracció de la regió genòmica del hit a estudiar: fastasubseq

Per poder extreure la regió genòmica primer hem d'indexar el fitxer fasta de la base de dades (genoma del protista):

fastaindex $genome_path/$stripped/genome.fa genomes/$stripped/$stripped.index




Seguidament obtenim una seqüència del fitxer fasta indexat, on es troba la regió que conté el hit, mitjançant la comanda fastafetch:

fastafetch $genome_path/$genome/genome.fa ./genomes/$genome/$genome.index $region > ./regions/$region.fa




Finalment podem realitzar el fastasubseq que consisteix en l'obtenció d'una subseqüència de la regió que acabem d'extreure.

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


Per defecte, start és l'inici del hit 20kb abans, i la longitud de la seqüència extreta és de 50kb.





6. Predicció de gens: Exonerate i Genewise

Per predir l'estructura genòmica de les seqüències d'interès obtingudes hem utilitzat exonerate i genewise. Ambdos programes serveixen per obtenir la mateixa informació, però utilitzen algorismes diferents, per tant, pot ser útil comparar els resultats obtinguts mitjançant algorismes diferents, ja que en algun cas un programa podria ometre un hit.





6.1. Exonerate



exonerate -m p2g --showtargetgff -q ./query_proteins/$proteins.fa -t ./subseq/${proteins}_${genome}_${hit}_genomic.fa $exhaustive > ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff




-m: especifica el model d'alineament (en aquest cas p2g “protein to gene”).

--showtargetgff: indica que els resultats incloguin el format gff.

-q: especifica la query que estem usant.

-t: especifica el target contra el que comparem la nostra query.



Aleshores extraiem la seqüència a partir de les coordenades que ens dóna el fitxer en format gff amb la comanda següent:

fastaseqfromGFF.pl ./subseq/${proteins}_${genome}_${hit}_genomic.fa exon.temp > ./cDNA/cDNA_${proteins}_${genome}_${hit}.fa




Seguidament traduïm la seqüència obtinguda:

fastatranslate --frame 1 ./cDNA/cDNA_${proteins}_${genome}_${hit}.fa > ./translations/aa_${proteins}_${genome}_${hit}.fa




Com que el resultat de l'exonerate ens dóna la seqüència codificant, hem utilitzat el primer marc de lectura (després d'haver parlat amb els tutors).





6.2. GeneWise



genewise -pep -both query_proteins/$proteins.fa subseq/${proteins}_${genome}_${hit}_genomic.fa > output_genewise/fasta/whole/genewise_fasta_${proteins}_${genome}_${hit}




genewise -pep -both -pretty -gff -cdna query_proteins/$proteins.fa subseq/${proteins}_${genome}_${hit}_genomic.fa > output_genewise/pretty/genewise_pretty_${proteins}_${genome}_${hit}




-pep: mostra la seqüència traduïda

-both: mostra ambdós sentits de traducció (forward i reverse).

-pretty: és una opció de visualització estètica

-gff: que inclogui el format gff

-cdna: que inclogui la seqüència de cDNA.



Com que del Genewise obtenim dues seqüències, hem utilitzat un programa (anomenat llarg.pl) que ens selecciona la seqüència més llarga.





Utilitzarem l'anotació de exonerate perquè, en general, s'obtenen scores més elevats. No obstant, també s'ha utilitzat el GeneWise per tots els genomes, per si en algun dels casos exonerate no donava bons resultats. Per tant, en algun dels casos haurem fet servir el GeneWise i el T-Coffee haurà estat basat en aquest. També és important tenir en compte que el T-Coffee no diferencia els gap dels codons stop, de manera que quan la selenocisteïna de la query estigui alineada amb un gap, el més probable és que es tracti d'una selenocisteïna en el genoma del protista.



7. T-Coffee

El programa T-Coffee permet realitzar alineaments globals entre seqüències similars, en aquest cas s'alinea la query i la seqüència predita.


t_coffee query_proteins/$proteins.fa ./translations/aa_${proteins}_${genome}_${hit}.fa > t_coffee_output_exonerate/gloaln_${proteins}_${genome}_${hit}.txt



8. Cerca de SECIS

Hem utilitzat SECIS Search que és un programa que permet predir els elements SECIS en la regió 3' UTR.



SECISearch.pl -p s -T -H -x ../subseq/${proteins}_${genome}_${hit}_genomic.fa




9. Automatització



blast_exonerate_genewise_tcoffee_secissearch.sh


És un script que automàticament executa tots els passos mencionats a l'enunciat (blast, exonerate, genewise, tcoffe i secisserach.sh). Aquest programa realitza un blast de totes les seqüències fasta /query_proteins/ contra tots els genomes. Aleshores, el script automàticament produeix una estructura de carpetes i guarda tots els resultats intermedis (de cada pas), de manera que cada annotació es pot reconstruir manualment.


El programa realitza un tblastn de totes les queries contra tots genomes en el mode normal així com amb l'opció m9. L'opció m9 s'utilitza per extreure hits significatius (e-value menor a 0,0001), és a dir, el programa guarda cada fila del ouput de m9 que conté un hit significatiu com un arxiu de text propi per al seu processament posterior. Després de indexar el genoma amb fastaindex, els “arxius hit” s'utilitzen per extreure una subregió de 50kb al voltant de la regió on es troba aquest hit mitjançant el fastasubseq. En el cas que el hit es trobi massa a prop de l'inici o del final del contig i per això, no es puguin extreure 50kb, el script automàticament reajusta les coordenades per ser capaç d'extreure una subregió.



Amb la subregió extreta es duen a terme exonerate, genewise i secissearch. Els outputs de genewise apareixen en mode “pretty” i “fasta”. El fitxer sortida de fasta es tradueix en sentit forward i reverse. L'annotació més llarga es selecciona mitjançant el script de perl llarg.pl. Els fitxers de sortida de exonerate venen amb un “gff-dump” inclòs. El fastaseqfromGFF s'executa en aquest arxiu per extreure una seqüència fasta de les coordenaes gff i es tradueixen mitjançant fastatranslate. Aquestes anotacions (exonerate i genewise) després són sotmeses a un alineament global mitjançant T-Coffee.


Els scripts uni.pl, llarg.pl i substitucio.pl els hem obtingut de les pàgines web dels grups de l'any passat. (link)


#!/bin/bash

#
# script to automaticaly blast different querys against diferent database
# script adaptet/modified from original found at http://bioinformatica.upf.edu/2011/projectes11/Bx/web/materialsimetodes.html 
#


# The directory "query_protein" must exist and contain the query proteins in fasta-format
# The perl-scripts uni.pl llarg.pl substitucion.pl must be present


mkdir output_blast
mkdir output_blast_detailed
mkdir genomes
mkdir regions
mkdir subseq
mkdir exonerate_output
mkdir cDNA
mkdir translations
mkdir t_coffee_output_exonerate
mkdir t_coffee_output_genewise
mkdir output_genewise
mkdir output_genewise/pretty
mkdir output_genewise/fasta
mkdir output_genewise/fasta/whole
mkdir output_genewise/fasta/longer/
mkdir output_secissearch



#specify genome paths

genome_path="/cursos/BI/genomes/protists/2012"

#read in proteins, cut away pathname

for proteins_fa in ./query_proteins/*; do {
basename_fa=`basename $proteins_fa`
protein=${basename_fa%.fa}
proteins=$proteins" $protein"
} done

#Two loops to combine genomes and proteins, extract query header

rm results_blast.txt 2>> error.txt

       
for protein in $proteins ; do {
	for genomes in A.laibachii_Nc14 A.rara C.fasciculata D.discoideum_AX4 D.fasciculatum F.cylindrus G.niphandrodes I.multifiliis_strain_G5 L.donovani_BPK282A1 L.tarentolae P.capsici P.polycephalum S.arctica T.congolense ; do {
		
		echo $genomes
		pattern=`cat ./query_proteins/$protein.fa | egrep ^">"`
		pattern=${pattern#*>}
		echo $pattern > temppattern
		pattern=`cut -f 1 -d " " temppattern`
		rm temppattern
		echo fasta-header is $pattern
                       
	#blasts standard and -m9 option
    
		blastall -p tblastn -i ./query_proteins/$protein.fa -d $genome_path/$genomes/genome.fa -o ./output_blast/$protein.$genomes.txt -m 9 -F F
		blastall -p tblastn -i ./query_proteins/$protein.fa -d $genome_path/$genomes/genome.fa -o ./output_blast_detailed/$protein.$genomes.txt -F F
			             
	#Directory and hit-file structure, extract target header to 
                       
		echo ${protein}_$genomes: >> results_blast.txt
		egrep "$pattern" output_blast/$protein.$genomes.txt | egrep e- >> results_blast.txt
		egrep "$pattern" output_blast/$protein.$genomes.txt | egrep e- > temp
		x=`egrep "$pattern" output_blast/$protein.$genomes.txt | egrep e-`
		if [ "$x" != "" ]; then
			mkdir hits 2>> error.txt
			mkdir hits/$protein 2>> error.txt
			mkdir hits/$protein/$genomes 2>> error.txt
			egrep -v '^#' temp > temp_stripped
			./uni.pl $protein $genomes > temp_stripped
		fi
		echo ${protein}_$genomes ready
	} done
} done

echo "###BLAST READY" >> error.txt

echo "#################"
echo "###BLAST READY###"
echo "#################"


#automaticaly anotate various blast -m9 hits with exonerate

# Genome indexing with fastaindex
	for genome in /cursos/BI/genomes/protists/2012/* ; do {
		
		stripped=`basename $genome`
		mkdir ./genomes/$stripped 
		echo $genome
		fastaindex $genome_path/$stripped/genome.fa genomes/$stripped/$stripped.index
	} done


# select if exonerate should run in exhaustive mode

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

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

		#chop of pathname-part from filenames
		              
			proteins=`basename $proteins`
			genome=`basename $genome`
			hit=`basename $hit`
			hit=${hit%.txt}
			echo ${proteins}_${genome}_${hit} 

		# detection of fasta-header to extract multifasta fragment from -m9 blast

			header_pattern=`cat ./hits/$proteins/$genome/$hit.txt | egrep ^">"`
			header_pattern=${header_pattern#*>}
			echo $header_pattern > tempheader
			header_pattern=`cut -f 1 -d " " tempheader`
			rm tempheader 2>> error.txt       

		# extract header and perform fastafetch
										        
			region=`cat $header_pattern ./hits/$proteins/$genome/$hit.txt | cut -f 2`
			fastafetch $genome_path/$genome/genome.fa ./genomes/$genome/$genome.index $region > ./regions/$region.fa

		# extract beginning and end of hit from -m9 blast

			beg=`cat $header_pattern ./hits/$proteins/$genome/$hit.txt | cut -f 9`
			ending=`cat $header_pattern ./hits/$proteins/$genome/$hit.txt | cut -f 10`

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

		# perform fastasubseq on fetched fragment, various conditions
			
			fastasubseq ./regions/$region.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg.temp

		# variable scenarios in case sequence is not long enough/too short...
													     
			error=`cat error_msg.temp`
			if [ "$error" != "" ]; then
			error_long=`cat error_msg.temp | egrep "before end"`
			error_short=`cat error_msg.temp | egrep "Unknown flag"`
																         
		# if ending is too near chop of lenght to max lenght
				if [ "$error_long" != "" ]; then
					long=$error_long
					long=${long##*(}
					long=${long%)}
					length=$(( $long - $start))
					fastasubseq ./regions/$region.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa
				fi
																	        
		# If beginning is too near set start to 0
				if [ "$error_short" != "" ]; then
					start="0"
					fastasubseq ./regions/$region.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_long=`cat error_msg_2.temp | egrep "before end"`
						if [ "$error_2_long" != "" ]; then
							long=$error_2_long
							long=${long##*(}
							long=${long%)}
							long=$(( $long - $start))
							fastasubseq ./regions/$region.fa $start $length > ./subseq/${proteins}_${genome}_${hit}_genomic.fa 2> error_msg_3.temp
						fi
					fi 
				fi
			fi

		
		#												   
		# perform genprediciton with exonerate
		#												   
			echo "performing: exonerate"
			
			exonerate -m p2g --showtargetgff -q ./query_proteins/$proteins.fa -t ./subseq/${proteins}_${genome}_${hit}_genomic.fa $exhaustive > ./exonerate_output/exonerate_${proteins}_${genome}_${hit}.gff 2> error_de_substitucio.temp
														         
		# run perl program "substitucion.pl" in case there are selenocysteins marked as U in the sequence

			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
					./substitucion.pl < ./query_proteins/$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} : >> results_exonerate.txt
			fi

			echo "performing: GFF-extraction and translation"

		# extract and translate cDNA from predictet gene-structure with fastaseqfromGFF.pl and fastatranslate in frame 1
																        																	   
			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 --frame 1 ./cDNA/cDNA_${proteins}_${genome}_${hit}.fa > ./translations/aa_${proteins}_${genome}_${hit}.fa
			fi
																		    
			echo exonerate ${proteins}_${genome}_${hit} ready
			
		# for each subseq of a hit, realize genewise annotations in "pretty" and fasta mode (check options)
	
		echo "*** performing : genewise "
	
		genewise -pep -both query_proteins/$proteins.fa subseq/${proteins}_${genome}_${hit}_genomic.fa > output_genewise/fasta/whole/genewise_fasta_${proteins}_${genome}_${hit}

		genewise -pep -both -pretty -gff -cdna query_proteins/$proteins.fa subseq/${proteins}_${genome}_${hit}_genomic.fa > output_genewise/pretty/genewise_pretty_${proteins}_${genome}_${hit}
				
																		       
		} done
	} done
} done

rm *.temp
echo "##################################"
echo "###EXONERATE AND GENEWISE READY###"
echo "##################################"

	#select longer genewise annotation

for genewise_out in ./output_genewise/fasta/whole/* ; do {
			
			
			temp_strp=`basename $genewise_out`	
			echo selecting longer hit from $temp_strp			
			./llarg.pl < $genewise_out > ./output_genewise/fasta/longer/$temp_strp.fa
			
} done

	# perform global alignment with t_coffee & secissearch

	#chop of pathname part from filesnames 

for proteins in hits/* ; do {
	for genome in $proteins/* ; do {
		for hit in $genome/*; do {
		
		proteins=`basename $proteins`
        genome=`basename $genome`
        hit=`basename $hit`
		hit=${hit%.txt}
		
		echo "performing: tcoffee"

	#perform global alignment with t_coffee

		if [ -f ./translations/aa_${proteins}_${genome}_${hit}.fa ]; then
			t_coffee query_proteins/$proteins.fa ./translations/aa_${proteins}_${genome}_${hit}.fa > t_coffee_output_exonerate/gloaln_${proteins}_${genome}_${hit}.txt
		fi
		
		if [ -f ./output_genewise/fasta/longer/genewise_fasta_${proteins}_${genome}_${hit}.fa ]; then
			t_coffee query_proteins/$proteins.fa ./output_genewise/fasta/longer/genewise_fasta_${proteins}_${genome}_${hit}.fa > t_coffee_output_genewise/gloaln_${proteins}_${genome}_${hit}.txt
		fi

		echo "*** performing SECISearch of subseq of $hit , from $proteins in $genome"
		cd output_secissearch
		SECISearch.pl -p s -T -H -x ../subseq/${proteins}_${genome}_${hit}_genomic.fa 
		cd ..
		
		} done 
	echo tcoffee of $genome ready
	
	} done
	echo tcoffee of $proteins ready
	
} done

echo "###SECISSEARCH READY###"

echo "##################################"
echo "###T-COFFEE FOR EXONERATE READY###"
echo "##################################"

echo NAME / EXONERATE / GENEWISE >>exonerate_genewise_compare.list

for proteins in hits/* ; do {
	for genome in $proteins/* ; do {
		for hit in $genome/*; do {
			
		proteins=`basename $proteins`
        genome=`basename $genome`
        hit=`basename $hit`
               
        if [[ "$hit" = "*" ]] ; then
        	echo shit ; else
        
        hit=${hit%.txt}
		
		if [ -f ./t_coffee_output_exonerate/gloaln_${proteins}_${genome}_${hit}.txt ]; then
			ali_exo=`head -n 1 ./t_coffee_output_exonerate/gloaln_${proteins}_${genome}_${hit}.txt | cut -d ',' -f 3,5` ; else
			ali_exo="NONE"
		fi
		
		if [ -f ./t_coffee_output_genewise/gloaln_${proteins}_${genome}_${hit}.txt ]; then
			ali_gen=`head -n 1 ./t_coffee_output_genewise/gloaln_${proteins}_${genome}_${hit}.txt | cut -d ',' -f 3,5` ; else
			ali_gen="NONE"
		fi
		
			
				
		echo  $proteins $genome $hit / $ali_exo / $ali_gen >> exonerate_genewise_compare.list
		fi
		
		} done
	} done
} done

echo "lei lei!"





uni.pl


#!/usr/bin/perl
use warnings;
use strict;
my $contador=0;

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


substitution.pl


#!/usr/bin/perl

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;

llarg.pl


#!/usr/bin/perl

#CAL QUE L'INPUT SIGUI EN FORMAT GENEWISE

use strict;

my @v;
my $inici=0;
my $longv=0;
my $c1=0;
my $par1=0;
my $par2=0;
my $temp1="";
my $temp2="";


while ()
{
    @v=split(//,$_);
    $longv=scalar(@v);
    $c1=0;
   
	while ($c1 < $longv){
		if ($v[$c1] eq ">"){
	    	$inici++;
		}
	
		if ($inici==2){
			$temp2.=$v[$c1];
			$par2++;
		}
		else{
	    $temp1.=$v[$c1];
	    $par1++;
		}
		
	$c1++;
	
	}

}

if ($par1 > $par2){
	print $temp1;
}
else{
	print $temp2;
}