B Bioinformati isk Xubuntu Perl scriptingg

  Opslaag til

Bioiinfform B maatisk k S iptiingg i Scri X bun Xub ntu u 122.100 (V Vol. 2) Henrik D. Dahlerup, D Februar F 20155 kontakt: henrik@ddah h hlerup.com

Bemærrk - Gammelddags bioinformaatiske Perl P 5 scripts, brugt b til at analyysere DNA, RNA og proteiner, p kan være udfordren nde at finde run ndt i blandt b andett fordi frremgangen under u bioinformatisk b scripting har æ ændret sig en smule over tid. Dettee dokument løsser ikke dette, men introducere i til en del af dennee tidlige tankega ang i Perl P 5 og derveed med held needsætter sidstnæ ævnte udfordring. Altt dette bliver heer udført i den n nye Perl P 5.15. Udoveer dette er der m mange steder anv vendt talesprog t hvor der ikke tagges hensyn til den bioinformatiske b termologi. Deenne tilgang ha ar til hensigt at gøre det nemmere at gennemskuee den praktiske p værd di af disse ttil tider abstrrakte bioinformatiske b værktøjer - brrugt til at analy ysere DNA, D RNA og proteiner. p  

Anvendte bøger: "Samss Publishingg"-books.

Copyright ©2014 © Henrik k D. Dahlerup p 

Bioinformatisk Xubuntu Perl scripting

Indhold  Overblik og opsamling med en bioinformatisk vinkel ............................................................................. 232 At analysere DNA, RNA og proteiner ................................................................................... 246 Grundlæggende Perl debugging ............................................................................................. 251 Mere konkret scripting i Xubuntu ............................................................................................................. 263 Hashes eller associative array'er ............................................................................................ 263 At omdanne dit output ........................................................................................................... 286 At arbejde med filsystemet .................................................................................................... 307 En opsamling på simple analyser af DNA, RNA og proteiner .............................................. 334 Funktioner der behandler forskellige forhold ........................................................................ 345 Funktioner som behandler scalar'er og lister ......................................................................... 369

Anvendelse

- Her er brugt simple og meget praktiske forklaringer som gør det nemt at forklare det fremlagte videre til andre, uden at tale forbi dem - hvis du for eksempel er hjælpelærer på et bioinformatisk kursus eller lignende. 

 Copyright ©2014 Henrik D. Dahlerup  231   

Bioinformatisk Xubuntu Perl scripting

Overblik og opsamling med en bioinformatisk vinkel  Her er en hurtig beskrivelse af udvalgte redskaber anvendt under bioinformatisk scripting i Perl. Enkelte nye meget anvendelige metoder er til stede i dette kapitel, men du kan nemt skippe dette kapitel hvis du ønsker at undgå gentagelser. Hashes En hash (også kaldt associative arrays) er en samling af nul eller flere par scalar-værdier kaldt keys og værdier. Værdierne er indekseret via keys. En associativ array-variabel starter med et procenttegn % efterfulgt af et tilladt variabel-navn. For eksempel mulige hash-variable navne er: %hash1 %genes_by_name Du kan assign'e en værdi til en key med en simpel assignment statement. For eksempel lad os sige at du har en hash kaldt %baseball_stadiums og en key Phillies som du ønsker skal blive assign'et med værdien Veterans Stadium. Dette statement gennemfører assignment'et: $baseball_stadiums{'Phillies'} = 'Veterans Stadium'; Bemærk at en enkelt hash værdi er refereret via en $ i stedet for en % ved begyndelsen af hashnavnet; Dette er den samme måde du refererede individuelle array værdier på ved at bruge en $ i stedet for et @. Du kan assign'e flere keys og værdier til en hash ved at placere deres scalar værdier i en liste, separeret af kommaer og omgivet af et par parenteser. Hvert par af scalar'er bliver til henholdsvis en key og en værdi i hash'en. For eksempel, kan du assign'e en hash med en tom liste: %hash = ( ); Du kan også assign'e en eller flere scalar key-værdi-par: %genes_by_name = ('gene1', 'AACCCGGTTGGTT', 'CCTTTCGGAAGGTC');

'gene2',

Der er en anden måde at gøre den samme ting på, som gør key-værdi-par'ene mere synlige. Dette opnår den samme ting som det forrige eksempel: %genes_by_name = ( 'gene1' => 'AACCCGGTTGGTT', 'gene2' => 'CCTTTCGGAAGGTC' ); For at få værdien koblet med en specifik key, kan du sætte $ bag hash-navnet og efterfølge den med et par curly braces { } som indeholder scalar-værdien af key'en: $genes_by_name{'gene1'} Dette returnerer værdien 'AACCCGGTTGGTT', som tidligere var givet værdien assign'et til key 'gene1' i hash'en %genes_by_name. At sammenkæde operatorer Sammenkædende operatorer er brugt til mønster-matchning, substitution og oversættelse og omskrivning af string'e. De er anvendt sammenholdt med regulær ekspressioner som specificer mønsteret. Her er et eksempel: 'ACGTACGTACGTACGT' =~ /CTA/  Copyright ©2014 Henrik D. Dahlerup  232   

Bioinformatisk Xubuntu Perl scripting

Den anden store regulære-ekspression ide er ændringer. Genstande som er opdelt via | metakarakterer match'er enhver af genstandene. For eksempel som i opstillingen nedenfor. Et simpelt program der matcher abc i en scalar variabel. #!/usr/bin/perl $alphabet = 'abcdefghijklmnopqrstuvwxyz'; if( $alphabet =~ /a(b|c|d)c/ ) { print $&; } $ program9_6 abc $ Dette producer outputtet: abc Eksemplet viser også hvordan parenteser grupper ting i en regulær-ekspression. Parenteserne er metakarakterer som ikke er match'et i en string; De grupper ændringen opgivet som b|c|d, som betyder enhver af b, c, eller d ved denne position i mønsteret. Siden b faktisk er i $alphabet ved denne position, er ændringen gennemført og hele mønsteret a(b|c|d)c, matcher i $alphabet. (En mere pointe: ab|cd betyder (ab)|(cd) og ikke a(b|c)d.) Den tredje store regulære-ekspression-ide er gentagelser. Dette er indikeret i et mønster med metakarakteren *. Når * er til stede efter en genstand, betyder den at genstanden kan hænde 0, 1 eller ethvert antal gange på dette sted i string'en. Så, for eksempel, alle af de følgende mønstermatch'ers vil have succes: 'AC' =~ /AB*C/; 'ABC' =~ /AB*C/; 'ABBBBBBBBBBBC' =~ /AB*C/; Kvantum: * + {MIN,} {MIN,MAX} ? Disse metakarakterer indiker gentagelsen af en genstand. * metakarakteren indiker nul, en, eller flere af de efterfølgende genstande. + metakarakterer indiker en eller flere af de efterfølgende genstande. brace { } metakarakterer lader dig specificere præcis antallet af tidligere genstande eller intervaller. For eksempel, {3} betyder præcis 3 af den efterfølgende genstand; {3,7} betyder 3, 4, 5, 6, eller 7 af de efterfølgende genstande; og {3,} betyder 3 eller flere af de efterfølgende genstande. Selve ? matcher ingen eller en af de efterfølgende genstande. At lade kvantummer matche minimalt ved brugen af ? Kvantummer som lige er vist er greedy (eller maksimale) som default, der betyder at de matcher så mange genstande som muligt. Nogle gange ønsker du et minimalt match som vil matche så få genstande som muligt. Du opnår dette via hver af følgende * + {} ? med et ?. Så, for eksempel, *? forsøger at matche så få som muligt, eller ingen, af de efterfølgende genstande før den forsøger at matche en eller flere af de efterfølgende genstande. Her er et maksimalt match: #!/usr/bin/perl 'same old same old same old' =~ /same.*old/; print $&; Dette matcher 'same' efterfulgt af .* (så mange karakterer som muligt), efterfulgt af 'old', og printer: same old same old same old  Copyright ©2014 Henrik D. Dahlerup  240   

Bioinformatisk Xubuntu Perl scripting

DB<12> W $receivingcommittment DB<13> c TA Debugged program terminated. Use q to quit or R to restart, use O inhibit_exit to avoid stopping after program termination, h q, h R or h O to get additional info. DB<13> Vent en gang! W kommandoen skulle indiker når $receivingcommittment ændre sin værdi. Men når programmet fortsætter med at køre med c kommandoen, kører den til enden, som fortæller at $receivingcommittment aldrig ændrede sin værdi. Så lad os starte programmet op igen og stop på linjen som ændre dens værdi: DB<13> R Warning: some settings and command-line options may be lost! Default die handler restored. Loading DB routines from perl5db.pl version 1.07 Editor support available. Enter h or 'h h' for help, or 'man perldebug' for more help. main::(program9_13:11): my $dna = 'CGACGTCTTCTAAGGCGA'; DB<13> w 45 42: next; 43 } elsif ($previousbase eq $base1) { 44: if ( /$base2/ ) { 45: print $base1, $base2; 46: $recievingcommitment = 1; 47 } 48 } 49: $previousbase = $_; 50 } 51 DB<14> b 46 DB<15> c TAmain::(program9_13:46): $recievingcommitment = 1; DB<15> n main::(program9_13:49): $previousbase = $_; DB<15> p $receivingcommittment DB<16> Hmm? Koden siger at den assign'er variablen med værdien 1, men efter du execut'ede koden med n og forsøgte at udprinte værdien printer den ingenting. Hvis du ser nærmere på programmet, kan du se at i linje 66 har du stavet $receivingcommittment forkert (Det vil sige: $recievingcommitment). Det forklarer alt; Dette fikses og nu køres den igen: $ perl program9_13 TAAGGCGA Løst! At løse et andet problem  Copyright ©2014 Henrik D. Dahlerup  260   

Bioinformatisk Xubuntu Perl scripting

12: @words = split (/\s+/, $input); 13: $playerlist[$playercount++] = $words[0]; 14: if ($words[1] eq "player") { 15: @fields = @playerfields; 16: } else { 17: @fields = @pitcherfields; 18: } 19: for ($count = 1; $count <= @words; $count++) { 20: $playerbase{$words[0].$fields[$count-1]} = 21: $words[$count-1]; 22: } 23: } 24: 25: # now, print out pitcher win totals and player home run totals 26: foreach $player (@playerlist) { 27: print ("$player: "); 28: if ($playerbase{$player."KEY"} eq "player") { 29: $value = $playerbase{$player."HR"}; 30: print ("$value home runs\n"); 31: } else { 32: $value = $playerbase{$player."W"}; 33: print ("$value wins\n"); 34: } 35: } $ program10_6 Gehrig player 2164 493 1991 .340 Ruth player 2503 714 2217 .342 Grove pitcher 616 300 3.05 Williams player 2292 521 1839 .344 Koufax pitcher 397 165 2.76 ^D Gehrig: 493 home runs Ruth: 714 home runs Grove: 300 wins Williams: 521 home runs Koufax: 165 wins $ Dette program er designet så det er nemt at lægge et nyt felt til databasen. Med dette i baghovedet, definere, linjerne 3 og 4 felterne som bruges når du bygger spiller og pitcher informationer. Linjerne 9-23 bygger databasen. Til at starte med starter linje 9 $playercount som 0; Denne globale variabel holder styr på antallet af spillere i databasen. Linjerne 10-12 læser en linje fra standard input filen, tjekker om denne fil er tom og fjerner ledende og efterfølgende mellemrum fra linjen og opdeler linjen i ord. Linje 13 lægger spillerens navn (Det første ord i input linjen) til listen af spiller-navne lagret i @playerlist. Tælleren $playercount har nu 1 lagt til sig; Dette reflektere det nye totale antal af spillere lagret i databasen. Linjerne 14-18 bestemmer om den nye spiller er en pitcher eller ikke. Hvis spilleren er en pitcher, er navnet på disse felter lagret i denne spillers informationer taget fra @pitcherfields; ellers er navnet taget fra @playerfields. For at simplificere den senere behandling er en anden array variabel, @fields, brugt til at lagre listen af felter som faktisk bliver brugt for denne spiller.

 Copyright ©2014 Henrik D. Dahlerup  280   

Bioinformatisk Xubuntu Perl scripting

At ændre header print formatet For at ændre navnet på print formatet som printer en side header for en specifik fil, ændre da værdien som er lagret i den specielle system variabel $^. Som med $~, er det kun værdien for den aktuelle default fil som kan blive ændret. For eksemel for at bruge print formatet MYHEADER som header filen for filen MYFILE, læg da følgende statement'er til: $oldfile = select(MYFILE); $^ = "MYHEADER"; select($oldfile); . Disse statement'er sætter MYFILE til at være den aktuelle default fil og ændrer header for MYFILE til at være i print formatet MYHEADER, og nul-stiller den aktuelle fil til dens originale værdi. At sætte sidelængden Som default er sidelængden 60 linjer. For at specificere en anden sidelængde kan du ændre værdien lagret i system variablen $=: $= = 66; # set the page length to 66 lines Denne assignment skal ske før den første write statement. Hvis side-længden er ændret i midten af programmet, vil den nye sidelængde ikke blive brugt før en ny side er startet. Opstillingen nedenfor viser hvordan du kan sætte sidelængden og definer et side-header print format for din output fil. Et program som sætter længden og printer formatet for en side. 1: #!/usr/bin/perl 2: 3: open (OUTFILE, ">file1"); 4: select (OUTFILE); 5: $~ = "WRITELINE"; 6: $^ = "TOP_OF_PAGE"; 7: $= = 60; 8: while ($line = ) { 9: write; 10: } 11: close (OUTFILE); 12: 13: format TOP_OF_PAGE = 14: - page @< 15: $% 16: . 17: format WRITELINE = 18: @>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 19: $line 20: . Lad os sige at du tilfører følgende input: $ program11_8 Here is a line of input.  Copyright ©2014 Henrik D. Dahlerup  300   

Bioinformatisk Xubuntu Perl scripting

programmet b ud, og så videre. I dette eksempel, er det alfabetiske bogstav a til z og cifrene er 0 til 9 tastet. Et program som demonstrere brugen af getc. 1: #!/usr/bin/perl 2: 3: &start_hot_keys; 4: while (1) { 5: $char = getc(STDIN); 6: last if ($char eq "\\"); 7: $char =~ tr/a-zA-Z0-9/b-zaB-ZA1-90/; 8: print ($char); 9: } 10: &end_hot_keys; 11: print ("\n"); 12: 13: sub start_hot_keys { 14: system ("stty cbreak"); 15: system ("stty -echo"); 16: } 17: 18: sub end_hot_keys { 19: system ("stty -cbreak"); 20: system ("stty echo"); 21: } $ program12_8 bcdefghijklmnopqrstuvwxyza1234567890 $ Subrutinen start_hot_keys ændre miljøet til at understøtte hot-key input. For at gøre dette, bruger den 2 call's til built-in funktionen system, som tager dens argument og execut'er det. Kommandoen stty cbreak fortæller systemet at det skal behandle inputtet en karakterer af gangen og kommandoen stty -echo fortæller systemet at den ikke skal vise karakterne tastet på keyboardet. Bemærk - Nogle gamle computere supporter ikke hot keys eller bruger ikke forskellige kommandoer til at skabe hot-key-miljøet. Hvis du bruger en computer der anvender forskellige kommanoder til at skabe miljøet, kan du stadig køre dette program; Du skal bare ændre stty kommandoen til det der virker på din computer. Loop'et i linjerne 4-9 læser og skriver en karakter per loop iteration. Linje 5 starter med at læse en karakter fra standard input filen ved at bruge getc. Linje 6 tester om karakteren læst er et backslash. Hvis det er det, terminere loop'et. Hvis karakteren ikke er et backslash, fortsætter programmet med linje 7. Denne linje oversætter alle alfanumeriske karakterer til det næste - bogstav eller nummer; for eksempel, oversætter den g til h, E til F og 7 til 8. Karakterne z, Z og 9 er oversat til a, A og 0. Linje 8 printer ud de oversatte karakterer. Netop fordi karakterne du skrev på keyboardet ikke er vist, får programmet det til at se ud som om dit keyboard er gået i stykker. Subrutinen end_hot_keys gendanner det normale miljø igen ved at fjerne system calls der er udført af start_hot_keys.

 Copyright ©2014 Henrik D. Dahlerup  320   

Bioinformatisk Xubuntu Perl scripting

CTATGCGGTA ATCTGACCTC Til at starte med er @sequences array fyldt med listen af korte DNA sekvenser. $index er startet med et nul og $last er startet med den største array index værdi koblet med @sequences array, som er 4. Loop forholdet er sandt for så lang tid værdien af $index er mindre end eller ligmed (<=) værdien af $last index'et. Her er nul mindre end eller ligmed 4, så statement'erne inde i loop'et execut'er. Værdien af elementet for array index'et $index er printet og $index scalar'en er øget ved at bruge ++ operatoren. $index scalar'en har nu værdien 1 og loop-forholdet er tjekket igen for at se om en anden iteration er ved at køre. Derfor er hvert element i denne array behandlet indtil værdien af $index overstiger værdien af $last index'et. Når dette sker stopper loopet. Perl bidrager også med en anden type loop som er foreach loop'et. Her vises hvordan du kan omskrive ovenstående program til at bruge foreach i stedet for while: 1: #! /usr/bin/perl -w 2: 3: # The program wich iterate over an entire array 4: # with 'foreach'. 5: @sequences = ( 'TTATTATGTT', 'GCTCAGTTCT', 'GACCTCTTAA', 6: 'CTATGCGGTA', 'ATCTGACCTC' ); 7: foreach $value ( @sequences ) { 8: print "$value\n"; 9: } $ program12_19 TTATTATGTT GCTCAGTTCT GACCTCTTAA CTATGCGGTA ATCTGACCTC $ Som det ses, er resultaterne produceret af program12_19 er præcis de samme som dem produceret af program12_18 programmet. Lister er nemme at arbejde med Husk at lister er scalar'er opdelt af kommaer. Da @sequences array først var defineret, var det foreslået at benytte (') således: @sequences = ( 'TTATTATGTT', 'GCTCAGTTCT', 'GACCTCTTAA', 'CTATGCGGTA', 'ATCTGACCTC' ); Brugen af enkelt quotes er ikke fuldstændig nødvendig. Hvis en værdi i listen er en string som ikke indeholder white-space, er enkelt quotes valgfrie. Siden korte DNAsekvenser ikke indeholder white-space, kan listen også blive skrevet således: @sequences = ( TTATTATGTT, GCTCAGTTCT, GACCTCTTAA, CTATGCGGTA, ATCTGACCTC ); Det er også tilladt at fjerne kommaerne (,) som separere liste elementerne. Lad os sige at qw operatoren bliver brugt, hvor qw er en forkortelse for "quote words". Følgende er en anden måde at specificere listen med genstande brugt til at starte @sequences array: @sequences = qw( TTATTATGTT GCTCAGTTCT GACCTCTTAA  Copyright ©2014 Henrik D. Dahlerup  340   

Bioinformatisk Xubuntu Perl scripting

Dette call'er til index som så hopper 5 karakter over før den starter med at søge efter foo i string'en lagret i $line. Lige som før, hvis index finder under-string'en, returner den det totale antal af karakterer der skal hoppes over (inklusiv antallet specificeret af det tredje argument til index). Hvis index ikke finder under-string'en på positionen af string'en som den søger i, returner den -1. Denne egenskab for index gør det muligt for dig at finde alle hændelser af en under-string i en string. Nedenstående opstilling er en modificeret version af opstillingen ovenfor som søger efter alle hændelser af the i en input linje. Et program som bruger index til at søge efter en linje flere gange. 1: #!/usr/bin/perl 2: 3: $input = ; 4: $position = $found = 0; 5: while (1) { 6: $position = index($input, "the", $position); 7: last if ($position == -1); 8: if ($found == 0) { 9: $found = 1; 10: print ("pattern found - characters skipped:"); 11: } 12: print (" $position"); 13: $position++; 14: } 15: if ($found == 0) { 16: print ("pattern not found\n"); 17: } else { 18: print ("\n"); 19: } $ program13_8 Here is the test line containing the words. pattern found - characters skipped: 8 33 $ Linje 6 af dette program call'er index. Netop fordi den startende værdi af $position er 0, starter det første call til index en søgning fra begyndelsen af string'en. Otte karakterer hoppes over før den første hændelse af the er fundet; Dette betyder at $position er assign'et med 8. Linje 7 tester om et match er fundet ved at sammenligne $position med -1, som er værdien index returner når den ikke finder string'en som den leder efter. Netop fordi et match er fundet, fortsætter loop'et med at execute. Når loop'et iterate'er igen, call'er linje 6 index igen. Denne gang hopper index ni karakterer over før den begynder at søge igen, som sikre at den tidligere fundet hændelse af the er sprunget over. Et total af 33 bytes er hoppet over før the er fundet igen. Igen fortsætter loop'et fordi forholds-ekspressionen i linje 7 er falsk. I den sidste iteration af loop'et, call'er linje 6 index og hopper 34 karakterer over før den starter søgningen. Denne gang er the ikke fundet og index returner derfor -1, og forholdsekspressionen i linje 7 er sand. Dette medfører at loop'et terminer.

 Copyright ©2014 Henrik D. Dahlerup  360   

Bioinformatisk Xubuntu Perl scripting

Dette program bruger chomp til at fjerne input linje-separatoren fra både en scalar variabel og en array variabel. Call'et til chomp i linje 5 konverter værdien af $scalar fra testing:: til testing. Antallet af karakterer fjernet, 2, er returneret af chomp og assign'et til $num. Call'et til chomp i linje 8 tjekker hvert element på @list. Det første element er konverteret fra test1:: til test1 og det sidste element er konverteret fra test3:: til test3. (Det andet element er ignoreret, netop fordi det ikke er termineret af end-of-line specifier'en). Det totale antal af karakterer fjernet, 4 (som er 2 fra det første element og 2 fra det sidste), er returneret af chomp og assign'et til $num. crypt funktionen crypt funktionen enkrypter en string ved at bruge NBS Data Enkryptering Standard (DES) algoritme. Den bioinformatiske syntaks for crypt funktionen er: result = crypt (original, salt); original er string'en som bliver enkrypteret, og salt er en karakter-string af to karakterer som definer hvordan du ændre DES algoritmen (for at gøre det mere svært at decode). Disse to karakterer kan være ethvert bogstav eller ciffer, eller en af . og / karakterne. Efter algoritmen er ændret, er string'en enkrypteret ved at bruge den resulterende key. result er den enkrypteret string. De første to karakterer af result er de to karakterer specificeret i salt. Du kan bruge crypt til at sætte op en password tjekker ens med dem brugt af UNIX login. Opstillingen nedenfor er et eksempel på et program som prompter brugeren for et password og sammenligner det med et password lagret i en speciel fil. Et program som spørger efter og sammenligner et password. 1: #!/usr/bin/perl 2: 3: open (PASSWD, "/u/jqpublic/passwd") || 4: die ("Can't open password file"); 5: $passwd = ; 6: chop ($passwd); 7: close (PASSWD); 8: print ("Enter the password for this program:\n"); 9: system ("stty -echo"); 10: $mypasswd = ; 11: system ("stty echo"); 12: chop ($mypasswd); 13: if (crypt ($mypasswd, substr($passwd, 0, 2)) eq $passwd) { 14: print ("Correct! Carry on!\n"); 15: } else { 16: die ("Incorrect password: goodbye!\n"); 17: } $ program14_2 Enter the password for this program: bluejays Correct! Carry on! $ Bemærk at password'et du skriver er ikke vist på skærmen. Linjerne 3-7 ekstraher et korrekte password fra filen /u/jqpublic/passwd. Dette password  Copyright ©2014 Henrik D. Dahlerup  372