############################################################### WARNING : FEU AIXO SEMPRE QUE OBRIU UN TERMINAL !!!! ############################################################### export LC_ALL='C' ############################################################### AN&AGRAVE;LISI D'UNA SEQUENCIA GENOMICA: APOCLUSTER ############Obtencio de la sequencia i caracteritzacio######### Anem a la pag. UCSC genome Browser cliquem DNA per obtenir sequencia de l'apocluster.Cliquem get DNA i obtenim sequencia en format fasta que guardarem en un fitxer anomenat "apocluster.fa". Primerament comptarem el nombre de nucleotids i la proporcio de cada un d'ells, passant el format fasta a un format tabulat en el qual tindrem en una primera columna el nom de la sequencia i en una segona columna la sequpencia com una unica cadena.Aixo ho farem mitjançant la comanda seguent: awk '{printf $1}' apocluster.fa > apocluster.tbl0 Un cop tenim la seq en forma tabular li introduim una tabulacio entre la primera i la 2a columna. Ara mesurarem l'allargada de la sequenica amb la seguent comanda: awk '{print length($2)}' apocluster.tbl0 El resutlat obtingut ha sigut que la nostra sequencia te una llargada de 500000 pb. Tot seguit mesurarem la frequencia dels diferents nucleotids amb la comanda: awk '{print $2}' apocluster.tbl0 | fold -1 | sort | uniq -c | \ gawk '{print $2, $1/500000}' Resultat obtugut: A 0,299292 C 0,22697 G 0,21088 T 0,262858 Per tal de poder treballar en els diferents programes, dividirem la nostra sequencia en dos trossos, que tindran un solapament de 50.000pb per tal que cap gens se'ns quedi partit. Per poder-la dividir utilitzarem el programa fastachunk, que l'exportarem amb la seguent comanda: export PATH=$PATH:/disc8/bin A continuacio l'apliquem amb la comanda: fastachunk apocluster.fa 0 275000 | fold -60 > apo1.fa fastachunk apocluster.fa 225000 275000 | fold -60 > apo2.fa ######Identificacio d'elements repetitius a la sequpencia######### Per poder buscar les regions repetides utilitzarem el programa EMBL RepeatMasker server.Mitjançant la opcio Browse introduim la nostra sequencia al programa i seguidament seleccionem la opcio fast. La resta d'opcions deixerem les que surten per defecte.Ho farem dues vegades, una per la sequencia apo1 i l'altre per apo2.Al finalitzar l'execucio del programa obtindrem 5 fitxers com a resultat de cada sequencia introduida. Fitxers apo1: apo1.seq.stderr progress report apo1.seq.out Annotation of masked sequence apo1.seq.cat crossmatch information apo1.seq.tbl summary of the repeat content apo1.seq.masked Masked Sequence Fitxers apo2: apo2.seq.stderr progress report apo2.seq.out Annotation of masked sequence apo2.seq.cat crossmatch information apo2.seq.tbl summary of the repeat content apo2.seq.masked Masked Sequence Guardarem cada un dels fitxer mitjançant "save link as" i amb el nom indicat anteriorment en dues carpetesdiferents,seq1 i seq2. Per treure informacio sobre les regions repetides haurem de treballar amb el fitxer apo1.seq.tbl i apo2.seq.tbl. Per apo1: * Contingut de GCs es de 47,07%. * La el nombre de bases de la sequencia emmascarada es de 115753pb (42,09%). El fet de tenir aquest elevat percentatge d'elements repetitius fara mes facil l'analisi de la seq. Per apo2: * Contingut de GCs es de 41,69%. * El nombre de bases de la sequencia emmascarada es de 107737 (39,18%). Veiem que el resultat obtinhut es bastant semblant al de apo1. Per visualitzar la distribucio de les repeticions utilitzarem l'eina gff2ps, la qual requereix que passem el fitxer apo1.seq.out a format gff.Per fer-ho utilitzem la seguent comanda: grep apo1 apo1.seq.out | \ awk 'BEGIN{ OFS="\t" } { print $5, $11, "repeat", $6, $7, ".", ".", "."; }' \ > apo1.seq.out.gff grep apo2 apo2.seq.out | \ awk 'BEGIN{ OFS="\t" } { print $5, $11, "repeat", $6, $7, ".", ".", "."; }'\ > apo2.seq.out.gff Un cop tenim els fitxers amb el format corresponent, l'exportarem el gff2ps,per tal de poder-la utilitzar des de l'emacs, tot i que tambe la prodriem utilitzar des de la seva pagina web.Comanda: export PATH=/disc8/soft/perl/bin/:/disc8/bin/:$PATH Un cop tenim els fitxers amb el format gff, ja podem executar el program,que ens realitzara un altre canvi de format, de gff a ps.La comanda es la seguent: gff2ps apo1.seq.out.gff > apo1.seq.out.ps gff2ps apo2.seq.out.gff > apo2.seq.out.ps Per poder visualitzar el fitxer final obtingut utilitzarem el Kghostview. Comanda: kghostview apo1.seq.out.ps kghostview apo2.seq.out.ps ###########Prediccio de gens Ab initio######################### Per realitzar la prediccio de gens Ab initio, utilitzarem tres programes diferents, el Geneid, el Genscan i finalment el fgenesh. GENEID(IMIM geneid server) Un cop estiguem al servidor del programa Geneid, farem la prediccio dels gens utilitzant els fitxers de lasequencia emmascarada en format fasta: * apo1.seq.masked.fa: El programa el farem correr tant en format gff(apo1.masked.geneid.gff)com en format geneid (apo1.masked.geneid.txt), ens els dos casos incloent CDS sequences. * apo2.seq.masked.fa: El programa el farem correr tant en format gff(apo2.masked.geneid.gff)com en format geneid (apo2.masked.geneid.txt), ens els dos casos incloent CDS sequences. Els resultats els guardarem mitjançant "save page as" i en els dos casos en format "text files". Per tal que posteriorment poguem utilitazar aquests fitxer, es molt important que eliminem qualsevol resta que hagi pogut quedar del format html (ho farem manualment, obrint els fitxers i eliminant aquelles parts del format html). GENSCAN El programa Genscan el farem servir des del shell. Per tal de poder-lo executar, utilitzarem la seguent comanda: genscan /disc8/bin/genscanparams/HumanIso.smat \ seq1/apo1.seq.masked.fa -cds > seq1/apo1.masked.genscan_full genscan /disc8/bin/genscanparams/HumanIso.smat \ seq2/apo2.seq.masked.fa -cds > seq2/apo2.masked.genscan_full Per poder passar els fitxers obtinguts a format gff, utilitzem la seguent comanda: gawk 'BEGIN{OFS="\t"} $2 ~ /Term|Intr|Init/ \ { print "apo1", "genscan", $2, start=($4<$5 ? $4 : $5),\ end=($5<$4 ? $4 : $5), $13, $3, $7, $1; }' apo1.masked.genscan_full | \ sed 's/\.[0-9][0-9]$//' > apo1.masked.genscan_full.gff gawk 'BEGIN{OFS="\t"} $2 ~ /Term|Intr|Init/ \ { print "apo1", "genscan", $2, start=($4<$5 ? $4 : $5),\ end=($5<$4 ? $4 : $5), $13, $3, $7, $1; }' apo1.masked.genscan_full | \ sed 's/\.[0-9][0-9]$//' > apo2.masked.genscan_full.gff FGENESH Per utilitzar aquest programa anirem a la seva pagina web i li introduirem els fitxer seguents: * apo1.seq.masked.fa: Clicant el boto "search" obtindrem un nou fitxer que anomenarem apo1.masked.fgenesh.txt * apo2.seq.masked.fa: Clicant el boto "search" obtindrem un nou fitxer que anomenarem apo2.masked.fgenesh.txt Per canviar el format dels fitxers obtinguts al format gf, utilitzarem la seguent comanda, semblant a l'utilitzada anteriorment pero una mica modificada.Comanda: gawk 'BEGIN{OFS="\t"} $4 ~ /CDSf|CDSi|CDSl/ \ { print "apo1", "fgenesh", $4, start=($5<$7 ? $5 : $7), \ end=($7<$5 ? $5 : $7), $8, $2, $3, $1; }' \ apo1.masked.fgenesh.txt | sed 's/.[0-9][0-9]$//' > apo1.masked.fgenesh.gff gawk 'BEGIN{OFS="\t"} $4 ~ /CDSf|CDSi|CDSl/ \ { print "apo2", "fgenesh", $4, start=($5<$7 ? $5 : $7), \ end=($7<$5 ? $5 : $7), $8, $2, $3, $1; }' apo2.masked.fgenesh.txt | \ sed 's/\.[0-9][0-9]$//' > apo2.masked.fgenesh.gff Un cop tenim les prediccions fetes pels tres programes, el que farem sera agrupar-les totes tres en un fitxer nou que tindra format ps, per tal de poder veure en una grafica com de diferents o semblants son.El programa utilitzat per canviar el format es el gff2ps i la comanda utilitzada es la seguent: gff2ps apo1.masked.geneid.gff apo1.masked.genscan_full.gff \ apo1.masked.fgenesh.gff > apo1.geneprediction.ps gff2ps apo2.masked.geneid.gff apo2.masked.genscan_full.gff \ apo2.masked.fgenesh.gff > apo2.geneprediction.ps Creiem que pot ser interessant d'introduir cadascuna de les sequencies apo1 i apo2 en l'agrupament de les prediccions. Utilitzarem els fitxers apo1.seq.out i apo2.seq.out, pero primer el canviarem a un format gff d'una sola linia.Per aixo utilitzarem la seguent comanda: gawk 'BEGIN{OFS="\t"} {$2="Repeatmasker"; print $0}' \ apo1.seq.out.gff > apo1.seq.out_singleline.gff gawk 'BEGIN{OFS="\t"} {$2="Repeatmasker"; print $0}' \ apo2.seq.out.gff > apo2.seq.out_singleline.gff Ara ja podem fer l'agrupament de les 3 prediccions mes la sequencia apo1 i apo2, respectivament en una sola linia. Comanda: gff2ps apo1.masked.geneid.gff apo1.masked.genscan_full.gff \ apo1.masked.fgenesh.gff apo1.seq.out_singleline.gff > \ apo1.seq.out.geneprediction_total.ps gff2ps apo2.masked.geneid.gff apo2.masked.genscan_full.gff \ apo2.masked.fgenesh.gff apo2.seq.out_singleline.gff > \ apo2.seq.out.geneprediction_total.ps El fitxer resultant de l'agrupament el visualitzarem amb el kghostview. ################Validacio dels gens predits###################### * VALIDACIO AMB ESTs HUMANS: Ara intentarem validar les prediccions dels gens obtingudes i investigarem si hi ha informacio addicional que ens pugui ajudar a inferir l'estructura exonica dels gens. Per fer la validacio,compararem la nostra sequencia contra una base de dades de EST d'huma, fent servir el BLAST del NCBI BLAST server.Com que estem introduint una sequencia força llarga i nomes ens interessen els alineaments identics, utilitzarem el MEGABLAST. Per fer correr el programa triarem les seguents opcions: - "choose database": EST_human. - "alignment view": pairwise. - "format": html. - "show": treuerem el tick de les tres primeres caselles. - "layout": one window. - "formatting": at the bottom. Els fitxers utilizats per correr el megablst seran els seguents: * apo1.seq.masked.fa * apo2.seq.masked.fa Els resultats obintguts els guardarem en dos formats: * apo1.masked.megablast.html * apo1.masked.megablast.txt * apo2.masked.megablast.html * apo2.masked.megablast.txt Com que els resultats son dificils d'analitzar, primerament els visualitzarem utilitzant el gff2ps, per tant haurem de canviar el format dels resultats al format gff. Per tal de visualitzar elresultat mitjantçant el gff2ps primerament hem de canviar el format dels fitxers a gff. Per a fer-ho utilitzem un programa que es diu PARSEMEGABLAST (elaborat per Josep F.Abril), que el guardarem en un fitxer amb el nom de parsemegablast.pl. Per poder-lo executar utilitzarem la seguent comanda: perl parsemegablast.pl -G seq1/apo1.masked.megablast.txt > \ seq1/apo1.masked.megablast.gff perl parsemegablast.pl -G seq2/apo2.masked.megablast.txt > \ seq2/apo2.masked.megablast.gff Per poder visualitzar el fitxer, el passarem a format ps amb la comanda seguent: gff2ps apo1.masked.megablast.gff > apo1.masked.megablast.ps gff2ps apo2.masked.megablast.gff > apo2.masked.megablast.ps Finalment el visualitzem amb el kghostview. Com que el que ens interessa es poder veure quines prediccions tenen evidencies, el que farem sera agruparel fitxer obtingut amb el megablast, que conte els EST, i el fitxer en el qual tenim les tres prediccions en un nou fitxer: * apo1.seq.out.geneprediction_total.ps * apo2.seq.out.geneprediction_total.ps Fins aqui hem pogut visualitzar tots els ESTs, pero els unics que ens donaran evidencies de l'existencia d'un gen son els spliced ESTs.Per aconseguir nomes aquests ESTs, farem servir el programa getsplicedhsp.awk. Per poder-lo aplicar ferem servir la comanda seguent: gawk -f getsplicedhsp.awk apo1.masked.megablast.gff > apo1.megablast.spliced.gff gawk -f getsplicedhsp.awk apo2.masked.megablast.gff > apo2.megablast.spliced.gff Un cop tenim el fitxer amb els spliced ESTs, els agruparem amb les 3 prediccions de gens, per veure quins gens dels predits tenen mes evidencies. Utilitzarem la seguent comanda: gff2ps apo1.megablast.spliced.gff apo1.masked.geneid.gff \ apo1.masked.genscan_full.gff apo1.masked.fgenesh.gff \ > apo1.spliced+genepredictions.ps gff2ps apo2.megablast.spliced.gff apo2.masked.geneid.gff \ apo2.masked.genscan_full.gff apo2.masked.fgenesh.gff \ > apo2.spliced+genepredictions.ps Per poder visualitzar els spliced ESTs nomes en la direccio forward, treurem el frame dels fitxers apo1.megablast.spliced.gff i apo2.megablast.spliced.gff amb la seguent comanda: gawk '{$7="."; print $0}' apo1.megablast.spliced.gff \ > apo1.megablast.spliced.noframe gawk '{$7="."; print $0}' apo2.megablast.spliced.gff \ > apo2.megablast.spliced.noframe.gff A continuacio agrupem el fitxer resultant amb les tres prediccions de gens en un nou fitxer en format ps per tal de poder-lo visualitzar: gff2ps apo1.masked.fgenesh.gff apo1.masked.geneid.gff \ apo1.masked.genscan_full.gff apo1.megablast.spliced.noframe.gff \ > apo1.pspliced+predictions+noframe.ps gff2ps apo2.masked.fgenesh.gff apo2.masked.geneid.gff \ apo2.masked.genscan_full.gff apo2.megablast.spliced.noframe.gff \ > apo2.pspliced+predictions+noframe.ps * VALIDACIO AMB TBLASTX ABANS DE FER TBLASTX farem un BLAT per trobar en quin cromosoma del genoma de Mus musculus podem trobar la regio homologa a la nostra. Per fer el Blat anirem al UCSC genome Browser. Ara agafem apo1.fa i ho dividim en trossos mitjançant el fastachunk. Primer tros: de 25000pb de la regio compresa entre el nucleotid 0 i el 160000 (farem 8 subsequencies). Per apo2.fa agafarem la regio del nucleotid 84000 i 224000 i la dividirem en diferents trossos tambe utilitzant el fastachunk. Aquestes regions son les que els gens predits no tenen suport de ESTs, per tant el que farem sera correr un TBLASTX per tal de trobar proteines homologues en el Mus musculus, que ens puguin donar suport als gens predits. Parametres:genoma de Mus musculus, DNA, Chrm score Per triar els resultats,agafarem aquell que tingui mes score. Resultats: apo1.1 1180 11698 24165 25000 87.1% 9 - 46378153 46388988 10836 apo1.2 395 2039 15142 25000 83.2% 9 - 46362272 46372484 10213 apo1.3 1582 9821 9654 25000 91.5% 9 - 46337755 46358396 20642 apo1.5 1530 6713 24021 25000 89.4% 9 - 46289676 46313650 23975 apo1.6 1006 3892 20174 25000 84.5% 9 - 46276550 46288652 12103 apo1.7 806 11747 24692 25000 86.7% 9 - 46248302 46257022 8721 apo1.8 1592 1155 24469 25000 86.6% 9 - 46233515 46246572 13058 Obtencio de la seq de Mus musculus a la pag.UCSC genome Browser, despres de fer el blat ens posem a browser. A continuacio cliquem DNA, hem canviat les coordenades (46.240.000-46.400.000) per obtenir el tros de la sequencia que correspongui amb la nostra regio d'interes. Seguidament anem a formatting options: - all upper case - mask repeats: to N - reverse complement (nosaltres tenim la sequencia en reverse i el programa la dona en forward, d'aquesta manera al fer el tblastx el mateix programa alineara les dues sequencies en reverse). cliquem get DNA i obtindrem la seq del Mus musculus que guardarem amb el nom de apo1.mouse.fa a la carpeta de treball. Farem un fastachunk de la sequencia apo1.fa, que ens agafi una regio de 200.000pb.Comanda: fastachunk apo1.fa 0 200000 | fold -60 > apo1.200.fa * Per apo1: A continuacio alinearem la nostra sequencia apo1.fa amb la del Mus musculus(apo1.mouse.fa) amb un TBLASTX LOCAL.Per poder-lo executar farem les seguents comandes: export PATH=$PATH:/disc8/bin/:/disc8/bin/ncbiblast/bin/:/disc8/soft/R/bin export BLASTMAT="/disc8/bin/ncbiblast/data" formatdb -p F -i apo1.mouse.fa blastall -p tblastx -d apo1.mouse.fa -i seq1/apo1.200.fa -e 0.1 > \ apo1.200.dbmouse.tblastx.out parseblast -S -G apo1.200.dbmouse.tblastx.out | sort > \ apo1.200.dbmouse.tblastx.gff gff2ps apo1.200.dbmouse.tblastx.gff > apo1.200.dbmouse.tblastx.ps gff2ps seq1/apo1.masked.fgenesh.gff seq1/apo1.masked.geneid.gff \ seq1/apo1.masked.genscan_full.gff apo1.200.dbmouse.tblastx.gff \ > apo1.gensprediction.tblastx.ps * Per apo2: Primer farem el blat per trobar la regio de Mus musculus. Agafarem la regio que va de 56000 pb i 224000pb i la dividirem en trossos de 25000 pb mitjançant el fastachunk per poder fer el blat. fastachunk apo2.fa 56000 25000 | fold -60 > apo2.1.fa - apo2.1 (56000-81000) apo2.1 2725 15 19292 25000 87.7% 9 - 46145062 46162967 17906 - apo2.2 (81000-106000) - apo2.3 (106000-131000) - apo2.4 (131000-156000) - apo2.5 (156000-181000) - apo2.6 (181000-206000)apo2.6 442 1269 8389 25000 81.4% 9 - 46060105 46065630 5526 - apo2.7 (206000-231000) Regio escollida de Mus musculus: 46.030.000-46.170.000. Ho guardem al fitxer apo2.mouse.fa Tallem la regio de la nostra sequencia que farem servir per fer el TBLASTX: fastachunk seq2/apo2.fa 56000 168000 | fold -60 > apo2.168.fa Comandes per fer el tblastx local: export PATH=$PATH:/disc8/bin/:/disc8/bin/ncbiblast/bin/:/disc8/soft/R/bin export BLASTMAT="/disc8/bin/ncbiblast/data" formatdb -p F -i apo2.mouse2.fa blastall -p tblastx -d apo2.mouse2.fa -i \ apo2.168.fa -e 0.1 > apo2.168.dbmouse2.tblastx.out parseblast -S -G apo2.168.dbmouse2.tblastx.out | sort | \ gawk '{ $4+=56000-1; $5+=56000-1; print $0 }' > \ apo2.168.dbmouse2.tblastx.gff gff2ps apo2.168.dbmouse2.tblastx.gff > apo2.168.dbmouse2.tblastx.ps gff2ps seq2/apo2.masked.fgenesh.gff seq2/apo2.masked.geneid.gff \ seq2/apo2.masked.genscan_full.gff apo2.168.dbmouse2.tblastx.gff \ > apo2.gensprediction2.tblastx.ps kghostview apo2.gensprediction2.tblastx.ps ##########################Analisi de proteines################################### * BLASTP Ara farem el BLASTP. Per a aixo hem de triar la seq de proteines amb la que ens quedem. La tria l'hem fet a partir dels resultats de l'agrupament dels tres programes de prediccio de gens junt amb les evidencies dels spliced ESTs.I hem decidit que per l'apo1 tenim: Parametres del programa: Swissprot,homo sapiens, max. alineaments * Apo1 Gen 4 (-) de Geneid (186.575 - 170.757) Gen 5 (-) de Geneid (192.536 - 201.251) * Apo2: Gen 4 (-) de Geneid (149.941 - 24.491) Gen 2 (-) de Genscan (11.754 - 7.865) * INTERPRO Ara farem interpro per trobar dominis funcionals a les proteines predites, de manera que ens confirmara la prediccio feta per BlatP.Per fer aixo anem a la pagina de Interpro. * apo1: - Gen4 de geneid:ens ha sortit un proteina que encara no ha estat caracteritazada. Hem mirat si aquella proteina estava en algun altre organisme i caratcteritzada, i hem trobat el de la Rattus norvegicus. Hem guardat la sequencia en format fasta amb el nom de apo1.gen4.Rattus norvegicus.fa.Hem vist que la funcio que tenia era la de NF kappa activating protein. Seguidament hem anat a la pagina del NCBI i hem buscat a gens aquesta proteina. Hem trobat que en el cromosoma 4 d'humans, hi havia una proteina amb la mateixa funcio, i ens hem guardat la sequencia d'aminoacids en format fasta amb el nom de nfkappa.crom4.fa.A continuaicio el que farem sera un clustalw per veure si les proteines s'assemblen i per tant poder definir la nostra proteina. - Gen 5 del geneid: hem trobat un ZINC-FINGER PROTEIN ZPR1 * apo2: - Gen2 de geneid:apopoliprotein (A1, A4) - Gen4 de geneid:apopoliprotein. Anem a la pagina del EMBL-EBI i anem a clustalw. Enganxem les tres sequencies en format fasta: - apo1.gen4.rata.fa - nfkappa.crom4.fa - gen 4 de la prediccio de geneid que trobarem a fitxer apo1.masked.geneid.txt com que no ens surt amb els 3, alinnearem de dos en dos: * gen4geneid amb rata * gen4geneid amb crom4huma Farem el clustalw local. el cridem amb aquesta comanda /disc8/bin/clustalw1.8/clustalw parametres: pairwise, blosum80 ( /disc8/bin/blast-2.2.13/data/BLOSUM80 ) Primer hem posat en un fitxer anomenat alineament.txt, les tres sequencies en format fasta.Com que els resultats no han sortit massa be, pq la sequencia de Rattus norvegicus es força diferent, farem un alinement nomes amb les sequencies d'huma. per aixo hem fet un nou fitxer anomenat queryhuma.txt amb aquestes dues sequencies i hem tornat a fer clutalw.I finalment hem fet Rattus norvegicus-huma en un fitxer anomenat queryrathum.txt. resultats: 3 sequencies: sequence.aln i sequence.dnd 2 sequencies: alinehuma.aln i alinehuma.dnd 2 sequencies: seqrathum.aln i seqrathum.dnd Per poder visualitzar les imatges i poder-les penjar a la pagina Web, farem la seguent comanda: convert -antialias -rotate 90 apo2.geneprediction.ps apo2.geneprediction.png ############################################################## Nota:sempre que s'ha executat un programa en perl previament s'han donat els permisos amb la comanda chmod u+x nomedelprograma.pl ##########################################################