Cada una de les selenoproteïnes s’ha emmagatzemat en un fitxer de text individual (exemple: SPG00000038_1.0.fa). A més, s’ha creat un altre fitxer que conté totes les selenoproteïnes humanes juntes (hs.aa).

Obtenció del genoma de Nomascus leucogenys

El genoma de Nomascus leucogenys que s’ha utilitzat es troba al directori:

cursos/BI/genomes/project_2013/Nomascus_leucogenys/genome.fa

El genoma de Nomascus leucogenys també es pot descarregar de la base de dades Ensembl, accedint a la pàgina web http://www.ensembl.org/Nomascus_leucogenys/Info/Index

Primerament interessa trobar quines regions del genoma de Nomascus leucogenys són susceptibles de contenir selenoproteïnes. Per saber-ho es fa un blast de les selenoproteïnes d'una espècie propera (Homo sapiens) contra el genoma problema (el de Nomascus leucogenys).

Requeriments

1. Tenir instal·lat a l’ordinador el software NCBI Blast
2. Un fitxer amb totes les selenoproteïnes humanes. Aquestes poden estar en forma aminoacídica o nucleotídica. En el nostre cas, tot i que disposàvem de les dues fonts, vam escollir la forma aminoacídica (hs.aa) ja que el programa s’executa a molta més velocitat.
3. Un fitxer amb la seqüència del genoma en el que volem trobar les selenoproteïnes (genome.fa). El tindrem en format FASTA, és a dir, nucleotídica.

Cluster de docència

Per accedir al cluster de docència només cal introduir la següent ordre i a continuació, quan és requerida, la contrassenya.

ssh UXXXX@luke.upf.edu

Caldrà transferir el fitxer amb les selenoproteïnes humanes “hs.aa” al cluster. Per fer-ho cal usar la següent comanda:

$ scp < origen > < destí >

On el destí serà: UXXXX@luke.upf.edu:/homes/users/UXXXX/hs.aa

Sistema de cues

Realitzar el blast és un procés que pot durar entre minuts i hores. Per aquest motiu vam usar el sistema de cues per dur-lo a terme. Aquest sistema funciona amb un fitxer de text que conté la descripció del càlcul que volem fer, que anomenem “tblast.job”. Aquest fitxer s’envia al sistema de cues, que s'encarregarà de reservar un processador per aquest treball i executar-lo. El fitxer del job és en realitat un script del shell amb un seguit de línies al principi que permeten proporcionar-li al sistema de cues la informació necessaria per executar el treball satisfactoriament.

Execució

Creació de tblast.job: emacs tblastn.job

Contingut de tblast.job:

#!/bin/bash

#$ -o hs_blast.stdout
#$ -e hs_blast.stderr
#$ -q llicen.q
#$ -N blastJOB
#$ -cwd

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

blastall -p tblastn -i hs.aa –d /cursos/BI/genomes/project_2013/Nomascus_leucogenys/genome.fa -m9 -e 0.001 -o hs.blast

tblast.job conté els diversos fitxers de sortida que ens interessen, els permisos per executar el programa i la comanda del blast, que consta de les següents parts:

Blastall

• -p tblastn: el tipus de blast que es vol usar (en el nostre cas, en comparar aminoàcids contra nucleotids usem tblastn)
• -i hs.aa: query de selenoproteïnes humanes
• -d /cursos/BI/genomes/project_2013/Nomascus_leucogenys/genome.fa: contra què comparem la query, és a dir, el genoma de Nomascus leucogenys
• -m9: format taula amb línies de comentari
• -e 0.001: constitueix un filtre pels valors de e value superiors a 0.001, que no serán inclosos en el blast final
• -o hs.blast: fitxer de sortida (tBLASTn final)

Finalment es va executar el programa amb la següent instrucció:

$ qsub tblastn.job

Quan el procés acaba s’obtenen 3 documents:

1. El fitxer de sortida amb el blast (tBLASTn). Aquest és el resultat que volem obtenir.
2. El fitxer amb la sortida estàndard “tBLASTn.stdout”, que també hauria de contenir informació.
3. El fitxer amb la sortida d’error estàndard “tBLASTn.stderr” que contindrà els possibles errors que s’hagin pogut produir.

Finalment s’exporta al directori actual del nostre ordinador l’arxiu tBLASTn amb la següent instrucció:

$ scp UXXXX@luke.upf.edu:/homes/users/UXXXX/resultat.tblastn .

S’aplica Fastaindex al genoma de Nomascus leucogenys amb la següent comanda:

$ Fastaindex /cursos/BI/genomes/project_2013/Nomascus_leucogenys/genome.fa genome.index

El fitxer resultant (genome.index) serà necessari per poder realitzar els següents passos.

Per tal de poder analitzar de manera més ràpida la gran quantitat de dades amb les quals treballem en aquest projecte, es va automatitzar el procés d’obtenció de dades. Aquesta automatització va ser possible amb la creació de dos programes en llenguatge bash. Primerament es va crear un fitxer de text amb la següent informació obtinguda a partir del blast: nom de la query, scaffold d’interès i punt d’inici i final de l'scaffold. A continuació, utilitzant el llenguatge bash, es va escriure un programa que incorporava tot el procés d’identificació de scaffolds, predicció d’exons (mitjançant exonerate i genewise), t-coffee i cerca d’elements SECIs. Finalment, vam crear un programa d’autoanàlisi que permetia executar automàticament el programa de comandes agafant la informació del fitxer de text.

Per a més informació, podeu consultar els programes i documents aquí.

Programes

autoanalysis.sh

Aquest programa permet l’execució de commands.sh de forma automàtica, extraient la informació necessària de cada scaffold emmagatzemada a l’arxiu blast.ordenat.

Requeriments

1. El programa autoanalysis.sh
2. El fitxer blast.ordenat
3. El programa commands.sh

Funcionament bàsic

autoanalisi.sh relaciona el fitxer de dades obtingudes del blast amb el programa commands.sh, que conté tots els mecanismes necessaris per l’anàlisi dels scaffolds. Quan commands.sh acaba d’analitzar un scaffold, autoanalysis fa que es torni a executar per a la següent línia de blast.ordenat, i així fins que s’arriba a l’última línia del fitxer de dades.

Execució

1. Donar permisos al programa amb la comanda chmod u+x autoanalysis.sh
2. Haver donat prèviament els permisos a commands.sh
3. Executar el programa al terminal amb la comanda ./autoanalysis.sh

Codi

#!/bin/bash
while read line; 
do 
   ./commands.sh $line;
done < blast.ordenat

commands.sh

Aquest programa permet executar per a cada scaffold les següents comandes: Fastafetch, Fastasubseq, Exonerate, Genewise, FastaseqfromGFF, Fastatranslate, T-coffee de les dues prediccions i cerca d’elements SECIs (SECISearch.pl). Ho executa de forma automàtica i guarda tots aquests resultats en una mateixa carpeta per a cada scaffold. D’aquesta manera el procés d’obtenció de dades es fa d’una manera molt més ràpida que manualment.

Requeriments

1. El programa commands.sh
2. El fitxer genome.index, que s’ha obtingut fent el Fastaindex i és necessari per fer el Fastafetch.
3. El fitxer blast.ordenat, que conté la informació seleccionada del tBLASTn que ens interessa. Aquesta informació correspon al nom de la query, nom del scaffold i la posició start i end del scaffold. Aquestes són les dades que ens permeten procedir amb la resta de protocol, tot i que hi ha informació del final.blast que també és interessant i que, per tant, també s’ha utilitzat a l’hora d’analitzar els resultats obtinguts (per exemple, l’e-value).
4. Fitxer fasta de la query ($QUERYNAME.fa), que conté la seqüència d’aminoàcids de la query. Hi ha tants fitxers fasta com selenoproteïnes analitzades. Aquests fitxers són necessaris per executar els programes de predicció de gens Exonerate i Genewise.
5. Tenir instal·lat a l'ordinador el software Exonerate.
6. Tenir instal·lat a l'ordinador el software Genewise.
7. Tenir instal·lat a l'ordinador el software T-coffee.
8. El programa FastaseqfromGFF.pl.

Funcionament bàsic

Primerament el programa assigna la següent informació als següents paràmetres:

Els quatre factors que llegeix commands.sh en un inici són els que componen cada línia del fitxer de dades: el primer ($1) correspon al QUERYNAME, el segon ($2) al SEQNAME (scaffold), el tercer ($3) i el quart ($4) al START i al END del scaffold.

D’altra banda assigna a $GENOMA el següent directori:

• “/cursos/BI/genomes/project_2013/Nomascus_leucogenys/genome.fa ”

Finalment assigna a $DIR (directori) el nom de la query i del scaffold ($QUERYNAME.$SEQNAME) de forma que crearà una carpeta per cada scaffold on s’emmagatzemarà tota la informació.

A) REGIÓ GENÒMICA DEL SCAFFOLD

Per tal de poder extreure primer la seqüència que conté la regió, i després extreure la regió d’interès de la seqüència, es necessita l'ajuda d'aquests dos programes que acompanyen a Exonerate.

- Fastafetch

S’ha utilitzat aquest programa per extreure el contig que conté el scaffold d’interès que quedarà registrat amb el nom del scaffold.fa ($DIR/$SEQNAME.fa). La comanda per a fer-ho és la següent:

Fastafetch $GENOMA genome.index $SEQNAME > $DIR/$SEQNAME.fa

- Fastasubseq

S’ha utilitzat aquest programa per extreure el fragment del contig on s’estima que hi hagi el nostre gen d’interès, que quedarà registrat amb el nom del scaffold.subseq ($DIR/$SEQNAME.subseq). Per assegurar que agafàvem el gen d’interès hem afegit 50Kb al START i al END de cada scaffold.

B) PREDICCIÓ D’EXONS

Per a predir els gens de les seqüències estudiades es van utilitzar els programes Exonerate i Genewise. És interessant executar els dos programes, ja que utilitzen algoritmes diferents, de manera que pot ser que la predicció no sigui exactament igual en els dos. Per tant, s’obté informació més completa executant-los tots dos.

- Exonerate

El programa Exonerate utilitza Fastaindex, Fastafetch, Fastasubseq i Fastatranslate. S’han predit els gens de les seqüències d’interès alineant la seqüència d’aminoàcids de la query inicial ($QUERYNAME.fa) amb la seqüència de DNA que s’ha extret a partir del Fastasubseq ($DIR/$SEQNAME.subseq). S’ha de tenir en compte que el fitxer fasta de la query no pot contenir cap “U” de la selenoproteïna, de manera que s’ha de canviar per una X (o un asterisc, tot i que en el nostre cas ho hem canviat per una X). També s’han d’eliminar els símbols de la seqüència extreta de selenoDB com per exemple @ o #, ja que sinó el programa Exonerate mostrarà errors. Aquests canvis es van realitzar manualment abans d’executar el programa commands.sh.

Cal explicar que es van executar i obtenir dos Exonerates diferents. El primer ($DIR/$SEQNAME.exonerate), conté l’opció |grep “cds”, que selecciona directament la part exònica de la seqüència, que és la que serà utilitzada per fer el T-coffee. D’aquesta manera, el programa T-coffee utilitza aquest fitxer, sense introns, per a alinear les seqüències. El segon Exonerate que vam realitzar ($DIR/$SEQNAME.exoneratesencer) no conté l’opció |grep “cds”, de manera que realitza l’Exonerate sencer (amb introns i exons). Això és necessari per a poder visualitzar l’Exonerate complet a l’hora d’interpretar els resultats i fer, per tant, un anàlisi més exhaustiu. Es pot veure, per exemple, la seqüència sencera d’aminoàcids, el número d’introns de la seqüència, etc.

- Genewise

Aquest programa utilitza un algoritme diferent al del programa anterior, de manera que obtenim una nova anotació del gen, que pot ser igual o no a la adquirida amb Exonerate. Els resultats dels dos programes, per tant, poden variar. Genewise permet comparar dues seqüències: la seqüència proteica de la query ($QUERYNAME.fa) amb la seqüència de DNA genòmic obtinguda del Fastasubseq ($DIR/$SEQNAME.subseq). Com passava amb l'Exonerate, s’ha de tenir en compte que, pel bon funcionament de Genewise, el fitxer fasta no pot contenir cap U (s’ha de substituir per una X) ni cap símbol com @ o # (s’han d’eliminar).

També es van executar i obtenir dos Genewise diferents. El primer ($DIR/$SEQNAME.genewise), conté l’opció |grep “cds”, de manera que es queda només amb els exons del gen. El segon Genewise generat ($DIR/$SEQNAME.genewisesencer) no conté l’opció |grep “cds”. Aquest últim Genewise generat no l’utilitzarem per a fer el T-coffee (utilitzarem $DIR/$SEQNAME.genewise) però el necessitem per a analitzar els resultats.

C) TRADUCCIÓ

- FastaseqfromGFF i Fastatranslate

Un cop obtinguts els fitxers amb el programa Exonerate i Genewise, procedim a traduir els exons amb el programa FastaseqfromGFF.pl i Fastatranslate. Fent anar el programa FastaseqfromGFF.pl s’emmagatzema el cDNA corresponent en un fitxer fasta ($DIR/$SEQNAME.cdnaE o $DIR/$SEQNAME.cdnaG) i, mitjançant Fastatranslate, en un altre fitxer fasta emmagatzemem la proteïna corresponent ($DIR/$SEQNAME.aaE.fa o $DIR/$SEQNAME.aaG.fa).

Els arxius que s’obtenen d’un Fastatranslate contenen sis seqüències fasta, que corresponen a les sis pautes de lectura possibles, tres forward i tres reverse. És important comparar els resultats obtinguts amb Exonerate amb els obtinguts amb Genewise, perquè pot ser que l’anotació d’un no sigui del tot precisa o correcta.

D) ALINEAMENT

- T-coffee

Un cop s’ha obtingut la seqüència proteica predita amb els programes esmentats anteriorment, executem T-coffee. Aquest ens permet alinear múltiples seqüències, de manera que l’utilitzem per alinear les seqüències proteiques obtingudes del genoma de Nomascus leucogenys amb les seqüències proteiques inicials (querys). D’aquesta manera, el programa ens permet veure si els aminoàcids de les dues seqüències comparades es corresponen i, per tant, inferir l’homologia entre aquestes. El resultat en tindrem en dos arxius, un pel t-coffe de la seqüencia obtinguda per Exonerate ($DIR/$SEQNAME.tcoff.E.html) i un per la del Genewise ($DIR/$SEQNAME.tcoff.G.html).

E) CERCA D'ELEMENTS SECIs

- SECISearch

Per saber si les proteïnes del genoma de Nomascus leucogenys predites contenen elements SECIs, hem executat el programa SECIsearch.pl.

Tots aquests passos descrits són els que s'han inclòs al programa commands.sh, de manera que aquest realitza tot el procés de manera automàtica. Per a què el programa commands.sh funcioni correctament cal:

1. 1. Donar els permisos al programa amb la comanda chmod u+x commands.sh
2. 2. Executar el programa
3. 3. Executar el programa al terminal amb la comanda ./commands.sh

Codi del programa commands.sh

#!/bin/bash
GENOMA="/cursos/BI/genomes/project_2013/Nomascus_leucogenys/genome.fa"

QUERYNAME=$1
SEQNAME=$2
START=$3
END=$4

DIR=$QUERYNAME.$SEQNAME

# Creem la carpeta de resultats per a cada query
mkdir -p $DIR 

echo Procesant $QUERYNAME $SEQNAME $START $END
echo
echo "Fastafetch"
Fastafetch $GENOMA genome.index $SEQNAME > $DIR/$SEQNAME.fa

let START=START-50000
let LENGTH=(END-START)+50000
echo "FASTASUBSEQ"
fastasubseq $DIR/$SEQNAME.fa $START $LENGTH > $DIR/$SEQNAME.subseq

# Executem Exonerate
echo "Exonerate"
export PATH=/cursos/BI/soft/Exonerate/Exonerate-2.2.0-i386/bin/:$PATH
Exonerate -m p2g --showtargetgff -q $QUERYNAME.fa -t $DIR/$SEQNAME.subseq |grep "cds" > $DIR/$SEQNAME.Exonerate

# Executem Exonerate sencer
echo "Exonerate SENCER"
export PATH=/cursos/BI/soft/Exonerate/Exonerate-2.2.0-i386/bin/:$PATH
Exonerate -m p2g --showtargetgff -q $QUERYNAME.fa -t $DIR/$SEQNAME.subseq > $DIR/$SEQNAME.Exoneratesencer

# Executem Genewise
echo "Genewise"
export PATH=/cursos/BI/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg
Genewise -pep -pretty -cdna -gff $QUERYNAME.fa $DIR/$SEQNAME.subseq|grep "cds" > $DIR/$SEQNAME.Genewise

# Executem Genewise sencer
echo "Genewise SENCER"
export PATH=/cursos/BI/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg
Genewise -pep -pretty -cdna -gff $QUERYNAME.fa $DIR/$SEQNAME.subseq > $DIR/$SEQNAME.Genewisesencer

# Traduim els exons Exonerate
echo "GFF"
FastaseqfromGFF.pl $DIR/$SEQNAME.subseq $DIR/$SEQNAME.Exonerate > $DIR/$SEQNAME.cdnaE 
echo "Fastatranslate"
Fastatranslate -F 1 $DIR/$SEQNAME.cdnaE > $DIR/$SEQNAME.aaE.fa

# Traduim els exons Genewise
echo "GFF2"
FastaseqfromGFF.pl $DIR/$SEQNAME.subseq $DIR/$SEQNAME.Genewise > $DIR/$SEQNAME.cdnaG 
echo "Fastatranslate2"
Fastatranslate -F 1 $DIR/$SEQNAME.cdnaG > $DIR/$SEQNAME.aaG.fa

# Executem T-coffee de les seqüències Exonerate i Genewise
export PATH=/cursos/BI/bin:$PATH
echo "TCOFFEE"
t_coffee $QUERYNAME.fa $DIR/$SEQNAME.aaE.fa -outfile $DIR/$SEQNAME.tcoff.E.aln
t_coffee $QUERYNAME.fa $DIR/$SEQNAME.aaG.fa -outfile $DIR/$SEQNAME.tcoff.G.aln

# Executem SECIsearch
echo "SECIsearch"
SECIsearch.pl $DIR/$SEQNAME.subseq > $DIR/$SEQNAME.SECIs





---

Recerca de selenoproteïnes de Pan troglodytes en el genoma de Nomascus leucogenys

S'ha complementat l'anàlisi buscant les selenoproteïnes de Pan troglodytes en el genoma de Nomascus leucogenys, ja que s'ha considerat que la proximitat filogenètica ho requeria. Aquestes selenoproteïnes han estat GPx1, GPx2, GPx3, GPx4, GPx5, GPx6, GPx7, SPS1, SPS2 i SelP, totes elles també presents en el genoma humà. Aquest estudi ha permès confirmar la presència d'aquestes selenoproteïnes en el genoma de Nomascus leucogenys.













---

Complementació de l'estudi mitjançant Selenoprofiles

S'ha utilitzat el software Selenoprofiles de detecció de selenoproteïnes, que identifica les selenoproteïnes presents en el genoma de Nomascus leucogenys comparant-lo amb tots els eucariotes. Els resultats obtinguts els hem utilitzat per verificar els nostres resultats (obtinguts mitjançant automatització i el posterior anàlisi dels resultats). La comanda usada és la següent:

nohup /cursos/BI/bin/Selenoprofiles /homes/users/U63758/results_folder -t

/cursos/BI/genomes/project_2013/Nomascus_leucogenys/genome.fa -s "Nomascus_leucogenys" -p eukaryotic -temp ~/temp &> /homes/users/U63758/results_folder/log

---

Realització del BLAST d'NCBI amb les nostres selenoproteïnes candidates

S'ha realitzat un BLAST comparant les selenoproteïnes que hem predit amb les seqüencies de proteïnes no redundants disponible a la base de dades NCBI. La query, en aquest cas, han estat els scaffolds destacats en cadascuna de les selenoproteïnes estudiades. D'aquesta manera es poden consolidar els resultats obtinguts, observant com aquests scaffolds es troben també presents en altres espècies.