1. Obtenció dels genomes de protists
2. Obtenció de les querys
3. tBLASTN i selecció dels hits
4. Extracció d'una regió genòmica d'interès
5. Generació d'una anotació: EXONERATE
6. Generació d'una anotació: GENEWISE
7. Alineament de proteïnes: T-COFFEE
8. Cerca dels elements SECIS
9. Cerca de la maquinària de transcripció
10. Cerca de similaritat al NCBI
11. Automatització

Els genomes dels catorze protists a analitzar han estat proporcionats pel Departament de Bioinformàtica de la Universitat Pompeu Fabra. Són accessibles des dels ordinadors del Campus Mar a través de la ruta:

~/cursos/BI/genomes/protists/2012

Les seqüències de selenoproteïnes de la família Thioredoxin Reductase (TR) les hem obtingut a partir de les bases de dades SelenoDB i NCBI. El criteri que hem emprat per seleccionar quines d'aquestes proteïnes utilitzem com a querys en la recerca de selenoproteïnes en els genomes dels protists ha estat el següent:

1. Utilitzar seqüències proteiques d'organismes filogenèticament propers als protists a analitzar, que hem extret de treballs d'anys anteriors versats sobre la família TR. Concretament hem emprat els resultats obtinguts per Antolin, A et al. (2011) amb Saprolegnia Parasitica i els resultats de del Amo, E et al. per Aureococcus Anophagefferens.

2. Utilitzar seqüències proteiques d'altres organismes més distanciats per tal d'estudiar un marc evolutiu més ampli. Hem partit de TR1, TR2 i TR3 d'Homo sapiens TR1, TR2 i TR3 de Mus musculus i TR1 de Caenorhabditis elegans. Hem fet un t-coffee amb totes aquestes seqüències per tal de descartar les querys molt semblants, i així tenir seqüències diferents com a querys per tal d'ampliar la búsqueda i que no s'escapin hits.

Per altra banda, de la família Lmsel1 hem emprat una query que ha estat proporcionada a partir d'un fitxer de la pàgina web de Bioinformàtica de Leishmania major i dues querys obtingudes del treball de Quevedo, M et al. (2010) de Leishamania infantum i de Leishmania mexicana.

Les querys seleccionades han estat emmagatzemades en fitxers FASTA (.fa)

BLAST v.2.2.13 (Basic Local Alignment Search Tool) és un programa bioinformàtic que permet realitzar alineaments de tipus local entre seqüències de DNA o de proteïnes. Dintre de les diferents modalitats de BLAST nosaltres hem emprat tBLASTn, que es basa en comparar una seqüència proteica (query) amb una sèrie de seqüències nucleotídiques que es troben a una base de dades. Per tal de realitzar l'alineament amb la query, el programa tradueix les seqüències de DNA a proteïna en els 6 possibles marcs de lectura (Open Reading Frame). Si algun dels alineaments compleix els paràmetres estadístics establerts, el programa reporta l'alineament com a hit. En aquest treball les bases de dades a emprar són els genomes dels protists i les querys són les seqüències seleccionades de TR i Lmsel1 citades anteriorment. El que hem fet ha estat alinear cadascuna de les querys amb els 14 genomes a estudiar. Si els genomes dels protists estan en format FASTA el programa BLAST no els pot utilitzar com a base de dades. Perquè siguin accessibles els podem formatejar amb la comanda formatdb. L'ordre que s'ha d'introduir al shell és la següent:

$ formatdb -i /cursos/BI/genomes/protists/2012/genome_name/genome.fa -p F -n genoma.fa

On l'argument -i indica la ruta d'arxius per accedir al genoma, el paràmetre -p F (false) ens informa que la base de dades no és un arxiu de proteïna, i on l'argument -n ens permet donar nom a l'arxiu de sortida que ja podrà utilitzar BLAST.

No obstant, no ha sigut necessari utilitzar aquesta comanda perquè els arxius que ens han facilitat ja estaven formatejats. Per extreure el software necessari per la utilització del BLAST cal executar aquestes dues comandes al shell:

$ export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
$ cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

A continuació, ja es pot executar el tBLASTn amb la següent ordre:

On -p indica la modalitat de blast emprada, -i especifica la query, -d indica la base de dades, -e indica el llindar d'E-value i -o dóna nom al fitxer de sortida de BLAST.

Un cop fet el tBLASTn, hem obtingut diferents hits per a cada query alineada amb cadascun dels genomes dels protists. Amb l'objectiu de poder treballar amb els resultats del BLAST directament, hem afegit l'argument -m8 a la comanda del BLAST. Amb aquest paràmetre hem obtingut un fitxer de sortida en forma de taula per a cada query on trobem per columnes: la identitat de la query, la regió on es troba el hit, el percentatge d'identitat, la longitud de l'alineament, els mismatches, els gaps openings, on comença l'alineament dins de la query, on acaba l'alineament dins de la query, on comença l'alineament dins del genoma, on acaba l'alineament dins del genoma, l'E-value i la puntuació (bit score) de l'alineament. La comanda a executar al shell en aquest cas és:

Dels diferents hits obtinguts hem seleccionat aquells que són estadísticament i biològicament significatius. BLAST és un programa basat en un algorisme heurístic, pel que no garanteix trobar la solució correcta. No obstant això, proporciona paràmetres que permeten jutjar el grau de significança dels resultats obtinguts. L'E-value és un d'aquests paràmetres i indica el nombre esperat de diferents alineaments (high-scoring pair o HSPs) que s'obtindrien per atzar amb una puntuació major o igual a un valor donat quan s'està buscant en una base de dades d'una determinada mida. Nosaltres hem considerat que un hit és significatiu quan aquest valor és igual o inferior a 10^-4 i a més a més quan compleixen un segon criteri: la selenocisteïna (U) de la query ha d'estar alineada amb un aminoàcid ja sigui una selenocisteïna, indicada com a codó stop, o un altre aminoàcid del genoma a analitzar.

En els passos següents hem treballat amb els hits seleccionats prèviament, els resultats BLAST -m8, extraient la regió genòmica on es troben per tal de limitar la búsqueda de les selenoproteïnes a regions petites dels genomes i d'aquesta manera facilitar la manipulació de la informació. Per fer l'extracció de les regions genòmiques d'interès, les quals tenen potencialment el gen que codifica per les proteïnes a estudiar, hem utilitzat una sèrie de programes que formen part del software Exonerate v.2.2.0. Per tal d'utilitzar-los s'ha d'introduir al shell:

$ export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH

El primer pas per l'extracció de les regions d'interès és indexar els genomes, per després poder utilitzar el programa Fastafetch, mitjançant la comanda Fastaindex:

$ fastaindex /cursos/BI/genomes/protists/2012/genome_name/genome.fa genomes_index/genome_name.index

On genome.fa correspon al genoma de l'organisme i genome_name.index al fitxer de sortida.

Els genomes es troben dividits en regions (contigs), i degut que només volem aquells alineaments que continguin U, hem modificat els resultats del blast -m8, per tal que només mostrin les region d'interès. A continuació, hem relacionat els contigs d'interès, amb la llista dels contigs del genoma resultant del Fastaindex, i els genomes. D'aquesta manera, mitjançant la comanda Fastafetch hem pogut generar arxius que contenen les regions genòmiques d'interès:

On genome.fa correspon al genoma de l'organisme, genome_name.index al genoma indexat, contig_name al contig d'interès, i genome_name.contig_name.RETALLAT.fa al fitxer de sortida amb el contig d'interès.

A continuació, hem delimitat encara més les regions d'interès mitjançant el programa Fastasubseq per tal d'obtenir seqüències encara més curtes. La comanda que cal escriure és:

On genome_name.contig_name.RETALLAT.fa és la regió que conté el contig d'interès, start fa referència a la posició d'inici del genoma a partir del qual s'agafarà la subseqüència, length indica la quantitat de nucleòtids que tindrà la subseqüència i genome_name.cotig_name.SUBSEQ.fa el fitxer de sortida que contindrà la subseqüència delimitada.

Els paràmetres start i length s'extreuen a partir dels fitxers que contenen els resultats del BLAST. El que hem mirat és on comença i on acaba l'alineament dins del genoma del protist de cada hit significatiu. Es considera que la posició més petita de l'alineament dins del genoma, ja sigui l'inici o el final segons si el sentit és forward o reverse , és la posició inicial de l'alineament. Per assegurar-nos que la subseqüència inclou el gen d'interès sencer, hem expandit els marges downstream i upstream de l'alineament, considerant la posició inicial de la subseqüència 10.000 nucleòtids upstream de la posició inicial i la posició final 10.000 nucleòtids downstream de la posició final de l'alineament. D'aquesta manera, ens queda una subseqüència d'uns 20.000 nucleòtids que presenta una mida suficient per no perdre informació.

En alguns casos, diferents hits d'un mateix genoma obtinguts a partir de diferents querys pertanyen un mateixa regió. El que hem fet per tal de treballar d'una manera més eficient, evitant tenir informació repetida, ha estat analitzar tots els hits significatius abans de fer l'extracció de la regió del genoma, per tal de només extraure una mateixa subseqüència una vegada.

Un cop tenim les regions genòmiques d'interès en fitxers específics el que hem fet ha estat executar el programa Exonerate. Aquest programa és una eina de comparació de seqüències que ha estat dissenyada per ser general i ràpida. Permet alinear seqüències utilitzant la programació dinàmica exhaustiva o bé una varietat de models heurístics. Per defecte, empra els models heurístics.
Exonerate, a més de produir un alinement més precís que BLAST, prediu l'estructura exònica de la seqüència introduïda (seqüència extreta amb el Fastasubseq) a partir de l'alineament d'aquesta amb la query. S'ha de tenir en compte que aquest programa no reconeix el símbol U (selenocisteïna) present a les querys, i per tant s'han de substituir per X.

La comanda que cal escriure al shell per tal d'executar el programa és:

On l'argument -m indica el model de l'alineament a utilitzar, en aquest cas p2g que significa protein to genome, és a dir, comparar una seqüència proteica amb una seqüència de DNA. El paràmetre -q especifica la query i -t la seqüència contra la que comparem la query. Per últim, query_name.CONTRA.genome_name.contig_name.exonerate és l'output que conté els resultats del programa.

El fitxer de sortida que s'obté amb la comanda mostra l'alineament entre la query i les diferents regions del la subseqüència. En algun cas en què hem vist que el tros de genoma alineat amb la query era molt curt, hem tornat a repetir l'exonerate però en la modalitat exhaustiva (E) que és més lenta, però el resultat és més òptim

De tots els fitxers obtinguts,hem seleccionat aquells en què la X de la query està alineada amb un asterisc *** (selenocisteïna) o amb una cisteïna.
Amb una segona comanda indicada a continuació, el que hem fet ha estat extreure les seqüències exòniques d'interès per després passar-les a cDNA:

On egrep selecciona les línies que tinguin el patró definit, -w especifica que el patró ha de ser una paraula sencera, exon és el patró a buscar,--showtargetgff especifica que el fitxer de sortida tingui el format GFF, i query_name.CONTRA.genome_name.contig_name.gff és el fitxer de sortida.

A partir de les seqüències exòniques en format GFF, hem utilitzat el programa FastaseqfromGFF.pl per tal d'obtenir el cDNA en format FASTA. Abans d'executar al programa cal indicar al shell el permís següent:

$ export PATH=/cursos/BI/bin:$PATH

L'ordre per executar el programa és:

A partir del cDNA, hem emprat el programa Fastatranslate per tal de traduir la seqüència de cDNA a proteïna. Aquest programa s'executa introduint la següent comanda al shell:

D'aquesta manera hem obtingut la proteïna amb els sis marcs de lectura possibles. Després d'analitzar els diferents ORFs hem seleccionat aquell que és correcte i ens interessa, és a dir, aquell que no té asteriscos que corresponen a codons stop (sense tenir en compte la selenocisteïna). Per generar un arxiu que només contingui el marc de lectura d'interès emprarem la següent comanda:

$ fastatranslate cDNA.fa -F/R X > fastatranslate_F/RX/query_name.CONTRA.genome_name.contig_name.prot

On -F inidica que és la seqüència forward, i -R indica que és la seqüència reverse, i X és el número de pauta a seguir.

De manera alternativa a l'Exonerate podem utilitzar el programa Genewise 2.2.0. Aquest programa ofereix el mateix que l'Exonerate; alinea seqüències de forma més precisa que el BLAST i prediu la seqüència exònica de la subseqüència donada. Per tant, l'utilitzarem per confirmar i contrastar els resultats obtinguts amb l'Exonerate. La seqüència emprada també és la obtinguda amb el programa Fastasubseq.

Abans d'executar el programa cal introduir al shell les ordres següents:

$ export PATH=/disc8/bin:$PATH
$ export WISECONFIGDIR=/disc8/soft/wise-2.2.0/wisecfg

La diferència principal amb l'Exonerate és que per executar-se correctament s'ha d'indicar la direccionalitat de la seqüencia genòmica; si és forward o reverse. Per seqüències que es troben en la strand positiva utilitzarem la comanda següent:

Per les que es troben en la strand reversa:

Si volem que ens busqui en ambdues cadenes:

On -pep especifica que mostri la seqüència peptídica a l'output, -pretty és que mostri l'alineament, -cdna que es mostri el cDNA, -gff que mostri la informació en format gff, -trev que el programa busqui a la seqüència reversa,-both que ens busqui a la seqüència positiva i la negativa,query_name.fa és la nostra query, genome_name.cotig_name.SUBSEQ.fa és la subseqüència i query_name.CONTRA.genome_name.contig_name.genewise és el fitxer de sortida.

Indicant tots aquests paràmetres aconseguim obtenir amb una sola comanda l'alineament entre la query i la subseqüència, el cDNA, i la seqüència peptídica predita. Tanmateix, per tal de poder establir comparacions amb els resultats del Genewise, els de l'Exonerate i les querys necessitem que els tres estigui en el mateix format i per això generarem outputs del Genewise que només continguin la seqüència proteica:

On els arxius generats query_name.CONTRA.genome_name.contig_name.prot, contenen les dues seqüències proteiques possibles. D'ambdues ens quedarem amb la correcta pauta de lectura.

El programa T-Coffee v8.99 (Tree-based Consistency Objective Function for alignment Evaluation) permet realitzar alineaments globals múltiples utilitzant un mètode progressiu que aparella seqüències similars i usa el resultat com a punt de partida del següent aparellament. Nosaltres l'hem fet servir per alinear la nostra query amb les selenoproteïnes o els homòlegs de cisteïna obtinguts a partir de l'Exonerate i del Genewise. Per tal d'establir totes les comparatives d'interès:

L'alineament resultant permet veure la similaritat que hi ha entre una proteïna predita i la query. A més a més, es mostra una puntuació de l'alineament que dóna una idea de la fiabilitat del resultat. Si les dues proteïnes s'assemblen molt i les selenocisteïnes (U o X) s'alineen entre elles tenim indicis per pensar que ens trobem davant de una selenoproteïna. Per tal de tenir més proves, hem realitzat dos passos més: buscar elements SECIS i la maquinària de transcripció necessària per la síntesis de selenoproteïnes en els genomes del protists.
En els casos on la selenocisteïna s'alinea amb una cisteïna, hem considerat que tenim indicis per pensar que ens trobem davant d'un homòleg amb cisteïna.

Per tal de buscar elements SECIS en els genomes dels protists hem emprat la versió web del software SECISearch v2.19. Aquest programa busca els possibles elements SECIS i valora la seva estructura tridimensional i les seves característiques termodinàmiques. La funció dels elements SECIS és reclutar la maquinària necessària per tal que el codó UGA sigui reconegut com selenocisteïna i no com a codó stop. Per tant, en principi la seva presència és necessària per a la síntesis de selenoproteïnes.

Els elements SECIS ens els organismes eucariotes estan presents a les regions 3' UTRs dels gens, generalment a uns 4000 nucleotids del codó UGA tot i que la distancia és variable. Per realitzar la seva cerca a la versió web l'input que hem emprat ha estat la seqüència obtinguda amb el Fastasubseq. Aquesta seqüència presenta 10.000 nucleòtids upstream del inici del gen, de tal manera que així ens assegurem que la regió que conté l'element SECIS estigui dins del input.

Al programa web s'ha d'indicar en quina de les dues cadenes es vol buscar l'element segons on es trobi el gen en potència. A més a més, SECISearch presenta una sèrie de patrons que permeten modelar el grau d'exigència en la cerca dels elements SECIS. Del més restrictiu al menys trobem strict, default i loose. El que hem fet és començar utilitzant el més restrictiu i si no hem obtingut resultats hem emprat els patrons menys restrictius progressivament. Tot i això, s'ha de tenir en compte, que al utilitzar sobretot patrons loose s'obtenen molts falsos positius.

L'output del SECISearch proporciona informació sobre la regió nucleotídica on es troba l'element SECIS, una imatge representativa de l'estructura tridimensional de l'element i un valor numèric anomenat COVE score. Aquest valor simbolitza la probabilitat de que la predicció sigui real. Valors de COVE score per sobre de 15 indiquen una alta probabilitat de que l'element SECIS que ha predit el programa existeixi i es trobi situat a la regió que descriu el programa.

La trobada d'un element SECIS a la regió 3' UTR d'una regió genòmica que pot ser potencialment una selenoproteïna proporciona més motius per pensar que el que estem analitzant es tracta realment d'una selenoproteïna.

Existeix una maquinària molecular necessària per la síntesis de selenoproteïnes. La troballa d'aquests elements en els genomes dels protists que potencialment tenen una selenproteïna pot ajudar a corroborar la hipòtesis de que es tracta d'una selenoproteïna.

Per buscar la maquinària als genomes hem realitzat tots els passos explicats fins ara excepte la cerca d'elements SECIS però amb les querys corresponents de cada element buscades a les bases de dades SelenoDB o NCBI:

1. SPS2:
   
   • Anopheles Gambiae
   • Drosophila Melanogaster
   • Pan Troglodytes
2. Pstk:
   
   • Drosophila Melanogaster
   • Homo sapiens
3. Sec43p:
   
   • Homo sapiens
4. SBP2:
   
   • Drosophila Melanogaster
   • Mus Musculus
   • Homo sapiens
5. SBP2:
   
   • Drosophila Melanogaster
   • Mus Musculus
   • Homo sapiens
6. SLA/LP:
   
   • Homo sapiens

Un altre element de la maquinària important per la síntesis de selenoproteïnes són els tRNAsec que els hem buscat mitjançant el servidor tRNAscan-SE Search Server. Aquests tRNA són els que presenten l'anticodó de TGA, que és TCA, i per tant transporten l'aminoàcid selenocisteïna al ribosoma.

Per acabar de corroborar els nostres resultats hem realitzat un BLASTp amb les proteïnes predites contra el conjunt no-redundant de les proteïnes disponibles al NCBI. Degut a que les proteïnes predites les hem obtingut mitjançant l'Exonerate i el Genewise, hem triat la que sortia millor de les dues per fer el BLASTp. Les seqüències de proteïnes que ens permetin trobar selenoproteïnes de la mateixa família en altres espècies (entre elles les querys) amb bons E-values ajudaran a donar encara més pes a la hipòtesis de que les proteïnes trobades són selenoproteïnes.

Degut que el nostre treball requereix de la repetició del procés anteriorment comentat amb moltes querys, hem automatitzat totes les comandes anteriors amb scripts del shell d'UNIX. Per dur a terme aquesta automatització hem creat tres scripts bash: l'automatic1.sh, l'automatic2.sh i l'automatic3.sh.

Aquests sistemes bash estan basats en l'indexat dels noms i les rutes d'accés dels diferents arxius que s'han d'utilitzar com a input en els diversos programes, per tal de facilitar la introducció al shell de la localització dels arxius. Degut a que inicialment partim dels genomes i les querys hem emprat una llista per cada pool d'arxius: genomes_list_and_info.tab i querys_list.tab,respectivament. Cal comentar que la llista de les querys l'hem elaborat manualment a partir de les que han estat seleccionades després de dur a terme els T-Coffees; d'altra banda, la llista dels genomes ens ha sigut proporcionada per l'equip de Bioinformàtica als ordinadors del campus mar a la ruta:

/cursos/BI/genomes/protists/genomes_list_and_info.tab

Automatic 1.sh

---

Així un cop obtingudes les dues llistes anteriors, hem pogut automatizar el procés del BLAST en aquest script. On dels resultats del BLAST i el BLAST -m8, per tal d'agilitzar l'anàlisi de resultats, hem eliminat els hits que no compleixen el criteri d'E.value mitjançant uns processos comentats a l'script, i ens hem quedat amb els hits útils. Els quals hem indexat a hitsbons_list.tab

Dels arxius obtinguts després d'aquest processament, n'hem fet una segona selecció, segons el criteri comentat de quedar-nos només amb els alineaments que continguin U. A més, també hem afegit als arxius del BLAST m8 una columna amb la posició d'inici i la llargada de cada contig, tenint en compte les premises ja esmentades.

#!/bin/bash

#Fem els exports necessaris per executar el blast.

export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

#Generem les carpetes requerides per guardar els resultats del blast.

mkdir blast blastm8 hitsbons blast_hits blastm8_hits

#Definim les variables necessaries per dur a terme el blast: la localització dels genomes i les querys, amb els respectius noms.


for genome_path in `grep 2012 /cursos/BI/genomes/protists/genomes_list_and_info.tab | cut -f 4`; do {
	genome_name=`grep $genome_path /cursos/BI/genomes/protists/genomes_list_and_info.tab | cut -f 1`;
	
	for query_path in `cat querys_list.tab | cut -f 2`; do {<
		query_name=`grep $query_path querys_list.tab | cut -f 1`; 

#Correm el blast.
			
	blastall -p tblastn -i $query_path -d $genome_path -e 0.0001 -o blast/${query_name}.CONTRA.${genome_name}.tblastn 
	
		blastall -p tblastn -i $query_path -d $genome_path -e 0.0001 -o blastm8/${query_name}.CONTRA.${genome_name}.tblastn -m8

		
#Generem uns duplicats dels arxius blast on només tindran contigut els que presentin hits.

	cat blast/${query_name}.CONTRA.${genome_name}.blastp | egrep Value > hitsbons/${query_name}.CONTRA.${genome_name}.tblastn
	} done
} done

#Eliminem tots aquells arxius que no presentin hits

find hitsbons -empty -exec rm '{}' \; 

#Guardem en una llista els noms dels arxius que presenten hits.

ls hitsbons > hitsbons_list.tab 

#Definim una variable que conté els noms dels arxius bons, i en fem una copia d'aquests a una altre carpeta, obtindrem els blasts que presentin hits.

for hitsbons in `cat hitsbons_list.tab | cut -f 1`; do {

	cp blast/$hitsbons blast_hits/$hitsbons

	cp blastm8/$hitsbons blastm8_hits/$hitsbons

} done

Automatic 2.sh

---

Després de la manipulació comentada hem seguit amb l'automatització dels paquets Fasta de l'Exonerate mitjançant aquest script. Aquí el que hem fet és treballar declarant les rutes d'accés dels arxius input i output segons les llistes de genomes_list_and_info.tab i querys_list.tab fins el Fastasubseq. Degut que en aquest programa havíem d'introduir diverses variables hem fet referència al fitxer Fastafetch mitjançant l'indexat automàtic del seu resultat, genomes_retallats_list.tab. Els resultats del Fastasubseq també els hem referenciat a la llista genomes_subseq_list.tab.

Com hem comentat anteriorment, les querys que s'han d'introduir a l'Exonerate no poden contenir Us, les hem canviat per Xs de forma automàtica mitjançant el programa perl canviaU.pl i també hem fet un indexat de la seva localització queys_canviades_list.tab

#!/bin/bash

#Executem els permisos necessaris per utilitzar els components del paquet exonerate i el canviaU.pl.

export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH
chmod u+x canviaU.pl  

#Generem les carpetes requerides per guardar els resultats del fastaindex, fastafetch i fastasubseq. Eliminem la carpeta de hitsbons. Renombrem la carpeta blastm8_hits a genomes_contigs, aquest nom orienta millor al seu contingut actual. 

rm -r hitsbons
mv blastm8_hits genomes_contigs
mkdir genomes_index genomes_retallats genomes_subseq

#Definim les variables necessaries per dur a terme els fastes: la localització dels genomes i les querys, amb els respectius noms.

for query_path in `cat querys_list.tab | cut -f 2`; do {
	query_name=`grep $query_path querys_list.tab | cut -f 1`; 

	for genome_path in `grep 2012 /cursos/BI/genomes/protists/genomes_list_and_info.tab | cut -f 4`; do {
		genome_name=`grep $genome_path /cursos/BI/genomes/protists/genomes_list_and_info.tab | cut -f 1`; 

#Correm el fastaindex, per indexar els genomes.

   		fastaindex $genome_path genomes_index/${genome_name}.index 


#Renombrem els arxius que contenen els contigs, per fer més clar el seu contingut. 

		mv genomes_contigs/${query_name}.CONTRA.${genome_name}.tblastn genomes_contigs/${query_name}.${genome_name}.contigs 

#Definim una variable que ens localitzi els contigs per poder dur a terme el fastafetch. 

		for contig_name in `cut -f 2 genomes_contigs/${query_name}.${genome_name}.contigs`; do { 

#Correm el fastafetch, obtenim els genomes retallats.
   
      			fastafetch $genome_path genomes_index/${genome_name}.index $contig_name > genomes_retallats/${genome_name}.${contig_name}.RETALLAT.fa  

#Fem una llista dels genomes retallats per facilitar-ne el seu accés.

echo "${genome_name}.${contig_name}""	""genomes_retallats/${genome_name}.${contig_name}.RETALLAT.fa""	""OK" >> genomes_retallats_list.tab 
	
#Definim l'inici i la llargada d'interès per cada contig, així com la loclaització i nom dels genomes retallats, per poder córrer el fastasubseq.

			for inici in `egrep $contig_name genomes_contigs/${query_name}.${genome_name}.contigs | cut -f 14`; do {

				for llargada in `egrep $contig_name genomes_contigs/${query_name}.${genome_name}.contigs | cut -f 13`; do {

					for genomes_retallats_path in `egrep OK genomes_retallats_list.tab | cut -f 2`; do {
						genomes_retallats_name= `egrep $genomes_retallats_path genomes_retallats_list.tab | cut -f 1`; 

#Correm el fastasubseq.

						fastasubseq $genomes_retallats_path $inici $llargada > genomes_subseq/${genomes_retallats_name}.SUBSEQ.fa 

#Fem una llista dels genomes processats amb el fastasubseq per facilitar-ne el seu accés. 

						echo "${genomes_retallats_name}""	""genomes_subseq/${genomes_retallats_name}.SUBSEQ.fa""	""OK" >> genomes_subseq_list.tab


					} done
				} done
			} done
		} done
	} done 

#Canviem les U presents a les querys per X, per tal de poder córrer l'exonerate.  

	./canviaU.pl < $query_path > TR_queryscanviades/${query_name}.fa 

#Fem una llista de les querys canviades per facilitar-ne el seu accés. 

	echo "${query_name}""	""TR_queryscanviades/${query_name}.fa" >> querys_canviades_list.tab

} done

Automatic 3.sh

---

En aquest script el que hem fet és automatitzar totes les comandes restants fins a l'obtenció de les proteïnes. Degut que l'automatització de l'Exonerate i el Genewise necessita com a input les querys i els genomes, hem utilitzat les respectives llistes ja esmentades. Dels resultats obtinguts hem dut a terme un indexat del l'output dels arxius .gff obtinguts mitjançant l'exonarate, en un document conjunt amb el nom i localització dels trossos genòmics creats amb el Fastasubseq; genomes_exonerate_subseq.tab. Aquesta relació entre documents ens permet automatitzar el fastaseqfromGFF.pl de manera més senzilla, d'aquest programa també hem fet un indexat dels seus resultats cDNA_list.tab. Cal esmentar que aquest darrer programa perl és accessible als ordinadors del Campus Mar, i no l'hem creat.

Així l'input del Fastatranslate l'hem extret del cDNA_list.tab, i de l'output també hem generat una llista, ,Prot_list.tab per tal de poder executar el programa perl canviaSTOP.pl de manera automàtica. Aquest programa l'hem creat per canviar la representació de codons stop en asteriscs a Xs, degut que d'aquesta manera ens sera millor la interpretació de les selenocisteines dels resultats als T-Coffees.

Per últim hem creat un darrer programa perl que ens selecciona de manera automàtica la sequència proteica d'interès dels resultats del Genewise només -pep. Així obtindrem arxius compatibles amb la informació que ens interessa per dur a terme el T-Coffee. Finalment per córrer aquest programa hem referit les rutes d'accés dels diferents inputs a llistes anteriorment generades.

#!/bin/bash

#Donem els permissos necessaris per executar tots els programes.

export PATH=/cursos/BI/bin/exonerate/i386/bin:$PATH
export PATH=/cursos/BI/bin:$PATH 
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg
chmod u+x canviaSTOP.pl trialamesllarga.ppl

#Generem les carpetes necessàries per emmagatzemar els arxius generats. 

mkdir exonerate exonerate_gff genewise genewise_resum cDNA fastatranslate fastatranslate_F1 t_coffee_exonerate_querys  genewise_resum_seleccionat t_coffee_genewise_querys t_coffee_genewise_exonerate fastarevcomp SECISearch SECISearch_rev 


#Definim les variables necessaries per dur a terme l'exonerate i el genewise: la localització dels genomes i les querys, amb els respectius noms.

for query_canviada_name in `egrep OK querys_canviades_list.tab | cut -f 1`;  do {
	query_canviada_path=`grep $query_canviada_name querys_canviades_list.tab | cut -f 2` ; 

	for genomes_subseq_name in `egrep OK genomes_subseq_list.tab | cut -f 1`; do {
		genomes_subseq_path=`grep $genomes_subseq_name genomes_subseq_list.tab | cut -f 2`; 

#Correm l'exonerate complet.

		exonerate -m p2g -q $query_canviada_path -t $genomes_subseq_path > exonerate/${query_canviada_name}.CONTRA.${genomes_subseq_name}.exonerate 

#Correm l'exonerate resumit, obtenim arxius .gff 
	
		exonerate -m p2g --showtargetgff -q $query_canviada_path -t $genomes_subseq_path | egrep -w exon > exonerate_gff/${query_canviada_name}.CONTRA.${genomes_subseq_name}.gff 

#Correm el genewise complet. 

		genewise -pep -pretty -cdna -gff -both $query_canviada_path $genomes_subseq_path > genewise/${query_canviada_name}.CONTRA.${genomes_subseq_name}.genewise 


#Correm el genewise resumit, obtenim arxius que contene les proteines. 

		genewise -pep -both $query_canviada_path $genomes_subseq_path > genewise_resum/${query_canviada_name}.CONTRA.${genomes_subseq_name}.genewise 

#Fem una llista dels arxius .gff de l'exonerate per facilitar-ne el seu accés. 

		echo "${query_canviada_name}.${genomes_subseq_name}""	""exonerate_gff/${query_canviada_name}.CONTRA.${genomes_subseq_name}.gff""	""${genomes_subseq_name}""	""${genomes_subseq_path}""	""OK" >> genomes_exonerate_subseq.tab

		} done
} done 

#Definim la localització i nom dels arixus exonerate .gff i dels genomes subseq, en relació amb aquests, per tal de poder córrer el fastaseqfromGFF.pl 

for exonerate_gff_name in `egrep OK genomes_exonerate_subseq.tab | cut -f 1`; do {
	for genomes_subseq_path in `egrep $exonerate_gff_name genomes_exonerate_subseq.tab | cut -f 4`; do {
		for exonerate_gff_path in `egrep $exonerate_gff_name genomes_exonerate_subseq.tab | cut -f 2`; do { 

#Correm el fastaseqfromGFF.pl per obtenir el cDNA. 

		fastaseqfromGFF.pl $genomes_subseq_path $exonerate_gff_path > cDNA/${exonerate_gff_name}.cdna 

#Fem una llista dels arxius de cDNA obtinguts per tal de facilitar-ne la seva manipulació 

		echo "${exonerate_gff_name}""	""cDNA/${exonerate_gff_name}.cdna""	""OK" >> cDNA_list.tab

		} done
	} done
} done 

#Definim la ruta d'accés i nom dels arxius de cDNA per tal d'executar el fastatranslate. 	

for cDNA_path in `egrep OK cDNA_list.tab | cut -f 2`; do {
	for cDNA_name in `grep $cDNA_path cDNA_list.tab | cut -f 1`; do { 

#Correm el fastatranslate, obtenim el cDNA. 

	fastatranslate $cDNA_path > fastatranslate/${cDNA_name}.prot

	fastatranslate $cDNA_path -F 1 > fastatranslate_F1/${cDNA_name}.prot 

#Fem una llista dels arxius de les proteines obtinguts per tal de facilitar-ne la seva manipulació. 
	
	echo "${cDNA_name}""	""fastatranslate_F1/${cDNA_name}.prot""	""OK" >> Prot_list.tab

	} done
} done	 

#Definim la ruta d'accés i el nom dels arxius de les proteines per executar el canviaSTOP.pl 

for prot_path in `egrep OK Prot_list.tab | cut -f 2`; do {
	for prot_name in `egrep $prot_path Prot_list.tab | cut -f 1`; do { 

#Canviem els STOPs (*) presents als arxius de les proteines per X, perquè aquests siguin alineats al t_coffee amb el canviaSTOP.pl 

	./canviaSTOP.pl < $prot_path > prots_canviades/${prot_name}.fa 

	} done
} done 

#Definim les variables necessaries per dur a terme el t_coffee: la localització dels genomes i les querys, amb els respectius noms. 

for genomes_subseq_name in `egrep OK genomes_subseq_list.tab | cut -f 1`; do {
		genomes_subseq_path=`grep $genomes_subseq_name genomes_subseq_list.tab | cut -f 2`;

	for query_canviada_name in `egrep OK querys_canviades_list.tab | cut -f 1` ;  do {
		query_canviada_path=`grep $query_canviada_name querys_canviades_list.tab | cut -f 2` ; 

#Correm el t_coffee que compara les querys canviades amb les proteines canviades, obtingudes a partir de l'exonerate. 

		t_coffee prots_canviades/${query_canviada_name}.${genomes_subseq_name}.prot TR_queryscanviades/${query_canviada_name}.fa -outfile tcoffee_exonerate_querys/${query_canviada_name}.${genomes_subseq_name}.CONTRA.${query_canviada_name}.tcoffee 

#Definim la variable prot que ens obrira els arxius proteïna vàlids de genwise_resum. Sobre aquests correrem el programa trialamesllarga.pl que permetrà quedar-nos amb la cadena vàlida per dur a terme el t_coffee, la més llarga.  

		prot=`cat genewise_resum/${query_canviada_name}.CONTRA.${genomes_subseq_name}.prot`; 
		
		if [ "$prot" != "" ]; then
		
		./trialamesllarga < genewise_resum/${query_canviada_name}.CONTRA.${genomes_subseq_name}.prot > genewise_resum_seleccionat/${query_canviada_name}.CONTRA.${genomes_subseq_name}.prot
			
		fi 

#Correm el t_coffee que compara les querys canviades amb les proteines genewise seleccionades.

		t_coffee genewise_resum_seleccionat/${query_canviada_name}.CONTRA.${genomes_subseq_name}.prot TR_queryscanviades/${query_canviada_name}.fa -outfile t_coffee/${query_canviada_name}.${genomes_subseq_name}.CONTRA.${query_canviada_name}.tcoffee 

#Correm el t_coffee que compara les querys canviades amb les proteines genewise seleccionades amb les proteines canviades de l'exonerate.

		t_coffee genewise_resum_seleccionat/${query_canviada_name}.CONTRA.${genomes_subseq_name}.prot prots_canviades/${query_canviada_name}.${genomes_subseq_name}.prot -outfile t_coffee_genewise_exonerate/${query_canviada_name}.${genomes_subseq_name}.CONTRA.${query_canviada_name}.${genomes_subseq_name}.tcoffee


	} done
} done 

Si et vols descarregar els nostres Automatics...

   Automatic1.sh                        Automatic2.sh                          Automatic3.sh