Thursday, August 9, 2012

SRS'de Coklu Arama Yapmak

Inceleme yapan scriptin en son hali, oncekilere gore daha fazla okuma inceliyor oldugu icin her okuma icin SRS uzerinde isim aramak oldukca zaman alan bir islemdi. Oyle ki, son inceleme 4 gun surdu.

Bunu azaltmak icin inceleme scriptini tamamen degistirdim. Oncelikle her zaman oldugu gibi esik degerini gecenleri aliyor ama direkt bunlarin ID numaralarini bir dizide (array) listeliyorum. Daha sonra bu listenin herbir elemanini boru karakteri ile ayirarak bir string haline getiriyorum. Son olarak bu stringi direkt getz komutuyla SRS'de aratip, tek seferde tum sayfada esik degerini gecen organizmalarin isimlerini alip, bunlari tek tek okuyup ayirarak hash icinde sakliyorum.

while (@params = $blast_obj->next_alignment(fields => ['name', 'identity', 'overlap'])) {

        $overlap_seen = $params[2];
        $mismatch_seen = $params[2] - $params[1];

        if ($overlap_seen >= $overlap_threshold and $mismatch_seen <= $mismatch_threshold) {

            if ($database eq "refseq_dna") {
                $params[0] =~ /:([A-Z|0-9]*\_[A-Z|0-9]*)/;
                push(@names, $1);
            }

        }

    }

    $ids = join("|", @names);

    if ( !($ids eq "") ) {

        open (PIPE, "getz '[$database:$ids]' -vf 'ORGANISM' |") or die $!;
        while(my $line = <PIPE>) {
            $line =~ /\t(.*)/;
            $organism_names{$1}++;
        }
        close PIPE;


Gene her ID numarasini duzenli ifade ile elde edip, bunu @names dizisine push komutu ile ekliyorum. Bu tumu icin bittiginde liste elemanlarini join komutu ile aralarinda "|" karakteri olacak sekilde birlestirip $ids stringini elde ediyorum. Sonra eger bu string bos degilse getz komutuyla arama yapip, elde ettigim cok satirli ciktiyi ozel bir sekilde okuyup her organizma ismini %organism_names hashine key olarak ekliyorum.

MegaBLAST Sonuclarini Incelemek - Parsing

Pipeline'da son asama, aranan dizilerin urettigi ciktilari baska bir script ile incelemek. Bu islemle herbir megablast dosyasi okunuyor, ve dizilerin name, identity, overlapping length gibi parametrelerinin degerleri saklanarak amaca yonelik sekilde ekrana yazdiriliyor.

Projemde HUSAR paketinde bulunan ve yukarida bahsettigim alanlari bana dizi olarak donduren Inslink adinda bir parser kullaniyorum. Bu parserin yaptigi tek sey, dosyayi okumak ve dosyadaki istenen alanlarin degerlerini saklamak.

Daha sonra ben bu saklanan degerleri, koda eklemeler yaparak gosteriyorum ve birkac ek kod ile de ihtiyacim olan anlamli sonuclar gosteriyorum.

Gene bu scripti de Unix uzerinde Emacs yazilimini kullanarak Perl dilinde yazdim.

Script parserdan gelen bilgileri kullanarak once esik degeri kontrolu yapiyor ve gecen tum okumalarin organizma ismini getz komutuyla alarak bir hash icinde saklayip sayarak bunlari daha sonra ekrana yazdiriyor. Boylece esik degerini gecen organizma sayisi uzerinden yorum yapabiliyorum.

Bu scripti projenin gelisme asamasinda cokca degistirdim. Farkli genomlari arastirirken farkli ihtiyaclar ciktigi icin her seferinde bir yenilik ile daha guvenilir hale geliyor. Daha oncekilero yazmayacagim ama su an son olarak kodladigim scripti paylasacagim.

Kalite Satirinin Degerlendirilmesi - Quality Filter

Kirleten organizma (konaminant) analizi yapacak olan pipeline'i daha fazla gelistirmek, daha anlamli sonuclar elde etmek icin ilk adimlara (henuz fastq dosyasini isliyorken) kalite filtresi eklemeyi dusunduk. Boylece belirli bir esik degerinden dusuk okumalari daha o asamadan filtreleyerek daha guvenilir sonuclar elde elebilecegiz.

Bu kalite kontrolunu fastq dosyasinda her okumanin 4. satirini anlayarak yapacagiz. Bu 4. satir (aslinda okumanin dizileme kalite skoru), cesitli dizileme cihazlari tarafindan cesitli sekillerde yaziliyor (kodlaniyor) ve bu kodlamadan tekrar kalite skorunu elde ederek filtreleme uygulanmasi gerekiyor. Bu yuzden oncelikle dizilme yapan cihazin kodlama (encoding) formatini bilmek ve bu skoru tekrar elde etmek gerekiyor. Henuz bu konuda bir fikrim yok. Dizileri aldigim bolume bir e-posta attim ve yakinda cihaz ile ilgili tum bilgileri edinecegim.

Bunu kendim de yapabilirim ancak yerine bir arac kullanmayi dusunuyorum. Simdilik bu araclar FASTX Toolkit icinde bulunan FASTQ Quality Filter araci, PRINSEQ ve FASTQC. Bu araclari tek tek arastirip hangisinin daha uygun oldugunu belirleyip onu pipeline'da kullanacagim.

Bu araclar temel olarak elimdeki fastq dosyasindan kalitesi dusuk okumalari cikarip atacak. Yani herhangi bir kesme yapmayacak, tamamen o okumayi siliyor olacak. Boylece elimde daha az ama guvenilir okuma kalacak.

Daha sonra bu okumalar arasindan megablast aramasi icin kullanilacaklari rastgele secmeyi dusunuyorum. Simdilik secimi dosyanin belirli noktalarindan baslayip 1000 okuma alarak yapiyorum.

Friday, August 3, 2012

Dorduncu Deneme Veriseti: Mus Musculus Genomu

Simdiye kadar ilk uc veriseti de insan genomuna aitti. Pipeline'i bu genomlarla deneyip, yer yer iyilestirmeler yaptim. Simdi ise baska organizmalarla da deneyip, daha fazla sonuc alip bunlari inceleyecegim ve gene gerekli iyilestirmeleri yapacagim.

Bu ilk farkli veriseti fareden geliyor. Mus Musculus tur adina ve ev faresi olarak yaygin isme sahip bu organizma da model organizma olarak calismalarda kullanildigi icin dizisi daha siklikla cikarilan diger bir organizma.

Bi dizilemeyi yapan, birlikte calistigim laboratuvardan cesitli BAM formatinda dizi dosyalari aldim. Bunlardan birini bam2fastq araciyla FASTQ formatina cevirerek dorduncu analizime baslayacagim. Bu dosyada toplam 37994473 dizi var ve bunlarin 6440475 tanesi eslesmeyen diziler.

Thursday, August 2, 2012

Inceleme Sonuclarini "Ambiguous" Olarak Ayirmak

Cesitli veritabanlarina karsi yaptigim aramalardan aldigim sonuclari incelerken, bunlari cesitli esik degerleri ile degerlendirmek ile beraber belirlenen esik degerlerinin uzerinde ya da altinda olan hitleri "Ambiguous" (belirsiz, cok anlamli) ya da "Unique" (essiz, tek) olarak ayirarak daha da anlamli hale getirmeye calisiyorum.

"Ambiguous" olarak, her bir megablast dosyasinda esik degerlerine uygun ancak birden fazla farkli organizmayi iceren hitleri etiketliyorum. Eger her esik degerine uygun hit, tek bir dosya icinde her zaman ayni organizmaya ait ise bu durumda yaptigim sey onu "unique" olarak etiketlemek.

Perl ile "Ambiguous" ve "Unique" icin birer hash tanimliyorum ve organizma isimlerine gore bu hashlerdeki anahtarlarin degerini artiritorum. Sonunda da bu iki etikete ait toplamlari liste olarak gosteriyorum.

Bunun bir ornegi aslinda ikinci veriseti incelemesinde var, ancak bunda "Ambiguous" hitleri sadece bir organizmaymis gibi kabul ediyorum ve hangi organizmalarin oldugunu gostermiyorum. Inceleme scriptinde yaptigim diger bir desigisiklikle daha fazla bilgi edinmek, "Ambiguous" hitlerin icerini ve organizma dagilimini gormek icin tek tek bunda bulunan organizmalari sayiyorum. Direkt buna ornek ucuncu verisetinin incelemesiyle gelecek.

Ikinci Veriseti Inceleme Sonuclari

Daha az eslenemeyen okumalara sahip ikinci verisetinin incelemesini tamamladim. Bu oncekine gore daha iyi bir dizileme ornegi oldugu icin aldigim sonuclar da oldukca tutarliydi. Insan genomuna ait bir diziden inceleme sonra asagidaki sonuclari elde ettim.

LIST OF ORGANISMS AND THEIR NUMBER OF OCCURENCES
Ambiguous hit   1323
Homo sapiens    312
Pan troglodytes 25
Pongo abelii    18
Nomascus leucogenys     17
Halomonas sp. GFAJ-1    7
Callithrix jacchus      4
Macaca mulatta  3
Oryctolagus cuniculus   2
Loxodonta africana      1
Cavia porcellus 1


"Ambiguous hit" tanimini baska bir yazida aciklayacagim. Simdilik sadece onlarin da bir hit oldugunu soylemem ve bu hitlerin de esik degerlerini gecmis oldugunu belirtmem yeterli olacak.

Yukaridaki sonuc RefSeq genomik veritabani aramasina dayali bilgiyi iceriyor ve birkac insana uzak okaryot ve bakteri disinda cogunlukla insan ve diger maymunlara ait sonuclar almam sonucun tutarli oldugunu gosteriyor. Aslinda daha onceki yazilara gozattiysaniz, sonuclarda "Homo sapiens" gormememiz gerektigi dusunebilirsiniz. Cunku bana gelen verisetinin insan genomuna karsi eslestirildigini ve eslesebilenlerin cikarildigini soylemistim, eger bana gelen veri her ikisini de iceriyorsa ben onlari baska bir araclar cikariyorum ki daha sonra karsima cikmasinlar ve olasi kirletici organizmalari direkt gorebileyim. Ancak buna ragmen sonuclarda "Homo sapiens" de goruyoruz. Aciklamayi eslestirme icin kullandigimiz araci inceleyerek bulmaya calistik ve su ki bwa araci eslesmeyen bazlara (mismatch) yeterince duyarli olmadigi bu sonucu aldik. Bu direkt olarak bwa ile ilgili ve onun kullandigi algoritmanin yarattigi bir durum. Elbette bizim icin bir sorun yaratmiyor, bunu bilerek zaten birtakim "Homo sapiens" hitleri bekliyoruz.

Bu verisetiyle birlikte inceleme yapan Perl scriptinde -- "Ambiguous hit" eklentisini de iceren -- bazi degisiklikler yaptim, onu ayrica aciklayacagim.

Thursday, July 26, 2012

Yeni Verisetinin Incelenmesi

Pipeline'i tasarlama asamasinda deneme amacli kullandigim onceki verinin cok kotu olmasi sebebiyle yeni bir veriseti aldim. Elbette deneme asamasinda birden fazla, farkli karakterlerde verisetleri kullanmak yararlidir. Ancak onceki veriseti anlamli birkac sonuc veremeyecek kadar kotuydu diyebilirim. Ayrintilarina buradan gozatabilirsiniz.

Yeni veriseti, gene bir insan genomu verisi ve BAM dosyasinin boyutu 1.8 GB ve icinde eslenebilen ve eslenemeyen okumalari bulunduruyordu. Ben bam2fastq araciyla hem bu BAM dosyasini FASTQ dosyasina cevirirken hem de eslenebilen okumalardan ayiklayarak 0.2 GB'lik FASTQ dosyasi olusturdum. BAM dosyasinda 26759102 adet okuma (dizi) bulunuyordu ve ayiklamadan sonra 735741 adet eslenemeyen okuma FASTQ dosyasina yazildi. Bu, tum dizinin % 2.7'sinin eslenemedigini gosteriyor ki oldukca normal bir durum.

Bu yeni veriseti ile birlikte scriptleri degistirmeye devam ediyorum. Ozellikle inceleme yapmak icin kullandigim scriptte bircok degisiklik yaptim. Bunu daha sonra yazacagim.

Birden Fazla Dizi Dosyalarindan MegaBLAST'i Calistirmak

Asagidaki scripti, pipeline'in MegaBLAST aramasini daha hizli yapabilmek icin dusundugumuz bir teknige uygun olabilmesi icin yazdim. Yaptigi sey, her okuma icin olusturulmus ve formatlanmis dizi dosyalarini kullanarak veritabanlarinda belirtilen baslangic noktasi ve okuma sayisi ile arama yapmak.

#!user/local/bin/perl

$database = $ARGV[0];
$dir = $ARGV[1]; #directory for sequences
$sp = $ARGV[2]; #starting point
$n = $ARGV[3] + $sp;

while (1){

    system("blastplus -programname=megablast $dir/read_$sp.seq $database -OUTFILE=read_$sp.megablast -nobatch -d");
    $sp++;
    last if ($sp == $n);

}


Burada her sey gercekten cok basit bir programlama ile isliyor. Tek yapilan sonsuz bir dongu olusturup bu dongu icinde tek tek dizileri system komutuyla MegaBLAST aracini kullanarak aramak. Scripte arguman olarak veritabani ismi, baslangic noktasi ve okuma sayisi veriliyor. Sonsuz donguden de, istenen okumalar arandiktan sonra cikiliyor.

Ayrica yeni scriptte megablast aracini artik farkli sekilde calistiriyorum. Bu sadece paketin en guncel halini kullanabilmek amaciyla degistirildi. Bununla birlikte "-nofilter" secenegini de kaldirdim. Boylece tekrar eden dizilerden (CTCTCT....) kurtulmus oluyorum.

Monday, July 23, 2012

Tek FASTA Dosyasindan MegaBLAST'i Calistirmak - Duzenli Ifadeler

Asagida MegaBLAST'i FASTA dosyasi okuyarak calistirmak ve sonuclari bir dizinde toplayabilmek amaciyla yazdigim Perl scripti ve onun aciklamasi var. Bu script tasarlamakta oldugum pipeline'in onemli bir parcasi. Bu script ilk yazdigim olan ve sadece bir FASTA dosyasi uzerinden tum okumalara ulasabilen script.

#!user/local/bin/perl
$database = $ARGV[0];
$fasta = $ARGV[1]; #input file
$sp = $ARGV[2]; #starting point
$n = $ARGV[3] + $sp;

if(!defined($n)){$n=12;} #set default number

open FASTA, $fasta or die $!;

while (my $line = <FASTA>){

    chomp $line;

    if ($line =~ /^>(\w+$sp)\s/){

        $name = $1;
        $input = $fasta."[".$name."]";

        system("megablast $input $database -OUTFILE=megablast/$name.megablast -nobatch -d");

        $sp++;
        last if ($sp == $n);
    }
}

close FASTA;


En son MegaBLAST yapabilmek icin sadece eslesemeyen okumalara sahip FASTQ dosyasini FASTA formatina cevirmistim. Simdi, bu FASTA dosyasiyla cesitli veritabanlarinda okumalari arayacagim ve cikti dosyalarini saklayacagim. Daha sonra da bu dosyalari inceleyecegim (parsing).

Bu projedeki FASTA dosyasi bir insan genomuna ait, bu yuzden oldukca buyuk bir dosya. Daha onceki yazilarimda toplamda yaklasik 5.5 milyon okumanin (ya da dizinin) oldugunu soylemistim. Ayrica dizilemeyi yapan makinenin ve dizileme islemleri genel olarak goz onunde bulunduruldugunda, deneme icin tum bu okumalarin MegaBLAST ile aratilmasi pratik degil. Bu yuzden bir kismini, ama gene de buyuk bir kismini ele alacagim. Ise yarar diyebilecegimiz kismi, dosyanin ortalarindaki okumalar olabilir. Bu okumalar digerlerine gore, daha kaliteli bir sekilde dizileniyor ve kaydediliyor. Bu yuzden ben de aramaya ortalardan baslayacak bir script yazdim. Scripte parametre olarak kacinci okumadan baslamasi gerektigini ve o okumadan itibaren o arama icin kac okuma aramasini istedigimizi ekledim. Elbette veritabani ismi, FASTA dosyasi ismi, her durumda olmasi gereken parametreler.

MegaBLAST arama scriptinde okuma adlarini duzenli ifadeler (regular expressions) kullanarak aliyorum, cunku eger gercekten herkes tarafindan kullanilabilecek bir pipeline tasarlamak istiyorsam, duzenli ifadeleri kullanmam sart. Ayrica, ismi aratirken, parametreden aldigim sayiyi da kullaniyorum. Kullanici scripte parametre olarak hangi okumadan baslamasi gerektigini belirttigi icin, okuma adlarini alirken o sayili okumayi bulup, devam etmem gerekiyor.

Ornegin eger okuma isimleri ">read_1", ">read_345", ">read_3000000" diye gidiyorsa, buyuktur isareti, birkac karekter veya alt cizgi ve bir degiskende saklanan okuma sayisini eslestirerek okuma adini buluyorum. Perl'de ise bunu, bir kontrol yapisi icinde soyle gosteriyoruz.

if ($line =~ /^>(\w+$sp)\s/){

}


Yukaridaki blokta $line degiskeni, while dongusu ile okuyor oldugum satira denk geliyor ve "=~" operatoruyle bu satirdaki baslagictan itibaren ">" karakterini direkt, sonra "read_" ifadesini "(\w+)" ile, okuma sayisini gene direkt "$sp" degiskeni ile karsilastiriyorum ve sonra da "\s" ile sonunda bir bosluk (space) karakteri oldugunu belirtiyorum. Hatirlayacak olursaniz, FASTA dosyasini olustururken, ismin sonunda bir bosluk birakip, daha sonra FASTQ dosyasindan gelen basligi ayni satira eklemistim. Bu ifade aciklamam gereken diger karakter ise "^". Bu karakterle, satirin basindan itibaren karsilastirmak istediginizi belirtiyorsunuz. Elbette duzenli ifadelerde olmasi gereken "/ ifade /" kurali gecerli ve buna uyuyorum.

Duzenli ifadelerle bu tarz karsilastirma yapmanin bircok yolu var. Bunlari baska bir yaziyla ya da yazi dizisiyle paylasmak istiyorum.

system("blastplus -programname=megablast $input $database -OUTFILE=megablast/$name.megablast -nofilter -nobatch -d");

Devaminda system komutu ile megablasti calistiriyorum. Ek olarak kullandigim parametreler, "-nofilter" (tekrar eden dizileri maskelemesin diye) ve "-nobatch".

Bu script bize her okumanin arama sonuclarini ayri .megablast dosyasi olarak veriyor ve onlar uzerinden inceleme yapabiliyoruz.

Unix'te Perl Ile Bir Komut Ciktisini Okumak ve Duzenli Ifadeler

Daha once organizma isimlerini duzenli ifadelerle nasil cikardigimi anlatmistim. Burada, gene benzer bir seyden bahsedecegim ancak bu biraz daha fazla, ozel bir teknikle Perl'de yapilan, veri tabanindan bilgileri birden fazla satir halinde cikti olarak aldigim icin gerek duydugum cok yararli bir yontem. Mutlaka benzerini baska amaclarla da kullanabilir, yararlanabilirsiniz.

Bu ihtiyac, HUSAR gurubu tarafindan olusturulan honest veritabaninin organizma isimlerini direkt sunmamasi ancak birkac satir halinde gostermesi sebebiyle dogdu. Asagida bunun ornegini gorebilirsiniz.

HUSAR %getz "[emblall:EU025714]" -f "ORGANISM"
OS   Salmo salar (Atlantic salmon)
 OC   Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
OC   Actinopterygii; Neopterygii; Teleostei; Euteleostei; Protacanthopterygii;
OC   Salmoniformes; Salmonidae; Salmoninae; Salmo.


Ilk satirda, organizmanin tur ismini ve parantez icinde de yaygin olarak kullanilan ismini goruyoruz. Ayrica TAB, "\t", ile birlikte basinda "OS" karakterleri bulunuyor. Iste bu ciktiyi Perl'de sanki bir dosya okuyormus gibi bir dongu icinde okuyarak, ilk satiri ve ilk satirdan da sadece tur ismini alacagiz.

elsif ($database eq "emblall") {
       $params[0] =~ /:([A-Z|0-9]*)\_?/;
       $params[0] = $1;
       open (PIPE, "getz '[$database:$params[0]]' -f 'ORGANISM' |") or die $!;
       while(my $line = <PIPE>) {
             $line =~ /^OS\s+(\w+)\s+(\w+)/;
             $params[0] = "$1 $2";
             }
       close PIPE;
      }


Yukaridaki elsif yapisi ile, oncelikle her zaman oldugu gibi getz komutunda kullanmak uzere okuma (ya da dizi) isminden erisim numarasini (anahtarini) bir duzenli ifade ile aliyorum. Bu duzenli ifade, veritabalarinin dizi isimlerini nasil sunduguna gore degistigi icin, her veritabani icin ayri bir elsif yapisinda, farkli duzenli ifadeler yaziyorum. Daha sonra bu erisim numarasini degiskende sakladiktan sonra open fonksiyonu ile "getz"in verdigi ciktiyi "aciyoruz" (sanki bir dosya okuyormus gibi...). Burada tek fark, getz komutunun sonuna boru (pipe), |, karakterini de ekliyoruz. Perl'de bu hile tam olarak boyle yapiliyor. Daha sonra gene her zaman oldugu gibi bir while dongusu ile (aslinda sadece bir cikti olan) "dosyamizi" satir satir okuyoruz. Sadece bir satirla ilgilendigimiz icin, bu satiri baska bir duzenli ifade ile ayiriyor ve sadece organizme ismini elde ediyoruz. Daha once de belrttigim gibi organizma isminin bulundugu satirda "OS" ve "\t" karakterleri ile birlikte, bir de bosluk, "\s", ve parantez icin organizmanin yaygin kullanilan ismi var. Elbette tur ismi cins ismi, bosluk ve tur adi seklinde yazildigi icin "(\w+)\s+(\w+)" ifadesiyle $1 degiskeninde orgizmanin cinsini, $2 ifadesinde de tur adini sakliyoruz. Duzenli ifadedeki "\w", bosluk disinda herhangi bir karakter anlamina geliyor ve sonuna arti, "+" ekleyerek bu karakterlerden birden fazla olabilecegini belirtiyoruz. Son olarak da "actigimiz" PIPE "dosyasini" close fonksiyonu ile kapatiyoruz.

Boylece, birkac satirdan olusan bir ciktidan, sadece organizmanin tur ismini elde ederek, raporda sunabiliyorum.

Duzenli Ifadeler ile Tur Ismini Elde Etmek

Projemin sonunda kullaniciya olasi kirleten organizmalarin adlarini (Latince tur isimleri) gosterecegim icin, MegaBLAST sonuclarindaki erisim numaralarini (accession number) kullanarak her dizi icin organizma adlarini elde etmem gerekiyor. Sequence Retrival System (SRS) adinda, HUSAR sunucularinda bulunan baska bir sistem ile bunu yapabiliyorum.

SRS'ten organizma adini ogrenebilmem icin Unix komut satirinda "getz" komutuyla birlikte veritabani ismi, erisim numarasi ve ogrenmek istedigim alani yazmam yetiyor. Asagida, bu isi yapabilen ornek bir kod bulabilirsiniz.

getz "[refseq_dna:NW_001840913]" -vf "ORGANISM"

Yukaridaki kod parcasinda "getz" komutumuz, "refseq_dna" NCBI'in Reference Sequence veritabani, "NW_001840913" erisim numarasi, "-vf" alanlari gostermemizi saglayan parametre ve son olarak "ORGANISM" de -vf parametresine bagli, bize organizma ismini verecek arguman.

Bu komut satiri, asagidaki ciktiyi veriyor.

REFSEQ_DNA:NW_001840913 Homo sapiens

Ve bu ciktiyi kullanarak, "REFSEQ_DNA:NW_001840913" ve hemen sonundaki TAB karakterinden kurtulmak ve kalan kismini bir degiskene atayip, gostermek gerekiyor. Bunu da asagidaki duzenli ifade ile yapiyoruz.

$params[0] =~ /\t(.*)/;
$params[0] = $1;


Bu kodda "$params[0]" getz komutunun ciktisini saklayan degisken. Bu degiskende saklanan tur ismini "\t" (TAB) karkaterinden sonra bulunan herhangi birden fazla karakter ifadesiyle "(.*)" alip, tekrar ayni degiskene, bunu atiyoruz.

Asagidaki kod, inceleme yapan scriptin bir bolumunden alinmis, refseq veritabaninda organizma ismini almak icin yazdigim kisma ait.

elsif ($database eq "refseq_dna") {
       $params[0] =~ /:([A-Z|0-9]*\_[A-Z|0-9]*)/;
       $params[0] = $1;
       $params[0] = `getz "[$database:$params[0]]" -vf "ORGANISM"`;
       $params[0] =~ /\t(.*)/;
       $params[0] = $1;
       }


Bu kod parcasini her megablast dosyasindan okudugumuz dizi ismi icin dondurdugumuzde, tur isimlerinin bulundugu bir listeyi elde edebiliriz. Bu elsif yapisi icinde gorunen diger duzenli ifade de erisim numarasini, getz komutunda kullanmak uzere elde etmemizi sagliyor.

Ayrica bu yapida getz komutunun ciktisini, direkt olarak bir degiskene backtick, "`", karakteri kullanarak atadigimi vurgulamak istiyorum. system komutundan sonra, bu Perl ile dis komutlari calistirmanin baska bir yolu.

Elbette bu yazida inceleme yapan scriptin diger kisimlari yok, onlari da ekledigimde tamamina bakabilirsiniz.

Thursday, July 19, 2012

Bir MegaBLAST Ciktisi Icerigi - RefSeq Veritabani

Asagida, deneme FASTA dosyasini refseq_genomic veritabaninda arayarak elde ettigim dosyadan, bir hitin ayrintilarini goruyoruz.

>>>>refseq_genomic_complete3:AC_000033_0310 Continuation (311 of 1357) of AC_000033 from base 31000001 (AC_000033
             Mus musculus strain mixed chromosome 11, alternate
             assembly Mm_Celera, whole genome shotgun sequence. 2/2012)
          Length = 110000

 Score =  115 bits (58), Expect = 4e-22
 Identities = 74/79 (93%), Gaps = 2/79 (2%)
 Strand = Plus / Minus

                        
                                          
Query: 1     ctctctctgtct-tctctctctctctgtctctctctctttctctctcttctctctctctc 59
             |||||||||||| ||| ||||||||| ||||||||||| |||||||||||||||||||||
Sbjct: 89773 ctctctctgtctgtctttctctctctctctctctctctctctctctcttctctctctctc 89714

                               
Query: 60    tttctctctgccctctctc 78
             ||||||||| |||||||||
Sbjct: 89713 tttctctct-ccctctctc 89696



Ayrintilarda, ilk olarak ">>>>" karakterleriyle hit ile ilgili baslik bilgisi veriyor. Bunda, kullanilan verutabani, ilgili organizme ve genin karsilastirildigi kromozomu ve dizileme bilgisi, tarihiyle birlikte veriliyor. Ayrica, uzunlugunu da gorebiliyoruz.  Bu alani, organizmanin ismini elde etmek icin kullanabilirim, bunun disinda projem icin baska amacla kullanmayacagim.

Daha sonra ayri maddeler halinda score, expectation value, identities, gaps ve strand bilgileri sunuluyor..

Son olarak da aranan dizi ile veritabaninda karsilastirilan dizinin karsilastirma sonucu, acik olarak baz eslestirilmesi gosterilerek veriliyor. Burada "Query" aradigimiz gen dizisi, "Sbjct" ise veritabaninda karsilastirilan dizi.

Monday, July 16, 2012

MegaBLAST Aramasini Hizlandirma

Son zamanlarda sadece farkli veritabanlarinda, MegaBLAST'i en cabuk ve etkili bir sekilde calistirmanin yolunu ariyorum ve FASTA dosyasi olusturma asamasinda, gercekten cokca ise yarayan bir yontem danismanim tarafindan geldi.

Daha once tum dizilerin bulundugu tek bir FASTA dosyasindan arama yapiyordum ve bu zaman kaybina yol aciyordu. Her ne kadar dosya bir sefer acilsa da her seferinde dosya icinde satirlara gidip onu okuman, zaman alan bir islem. Bunu, dosyadaki her okumayi, ayri bir FASTA dosyasi haline getirerek cozduk.

Su an elimizde 5.5 milyon ayri dosya olsa da, eskisine gore cok daha hizli arama yapabiliyoruz.

Buna gore MegaBLAST calistirma scriptim de asagidaki gibi degisti.

#!user/local/bin/perl

$database = $ARGV[0];
$dir = $ARGV[1]; #directory for sequences
$sp = $ARGV[2]; #starting point
$n = $ARGV[3] + $sp;

while (1){

    system("megablast $dir/read_$sp.seq $database -OUTFILE=read_$sp.megablast -nofilter -nobatch -d");
    $sp++;
    last if ($sp == $n);

}

Veritabanina Gore Bir Komutun Calisma Suresi - CPU Runtime

Calisilan dosyalar, veritabanları buyuk olunca ve yeterince bilgisayar gucune sahip olmayınca, her seyden once olcmemiz gereken nasil en etkili ve kisa surede sonucu alabiliyor olmamizdir.

Özellikle projemde, farkli veritabanları ve farkli parametreler kullanarak, bunları arastiriyorum.

Şimdilik dort veritabani deniyorum, bunlar: nrnuc, ensembl_cdna, honest ve refseq_genomic. Ayrica, bunu farkli iki kelime uzunluğuna gore de yapacagim. Kelime uzunluğu (word size) MegaBLAST'in ararken tam olarak eslestirecegi baz cifti sayisi. Yani elimde 151 baz ciftine sahip bir dizilim varsa, ve eger kelime uzunluğu 50 olarak belirlenmişse, bu 151 baz cifti icinden herhangi bir yerden baslayan ama arka arkaya en az 50 bazin dizilendiği kisimlar aranacak.

Bu veritabanlarinin farkli sonuclar vermesi elbette sahip oldugu dizi sayisina ve bunlari buyuklugune bagli. Projemde ben DNA dizisinden, herhangi bir tahminim olmadan, direkt olarak veritabanindan olasi organizmalari aradigim icin bana yeterince genis ve iyi sonuclar verebilecek bir veritabani gerekiyor. Belki iki veritabanini birlikte de kullanabilirim, ancak tum bunlari belirlemeden once veritabanlarini tek tek arastirmam gerekiyor.

Bu veritabanlari ile ilgili arastirmalarimi da gene birer yazi olarak yayinlayacagim.

Friday, July 13, 2012

Veritabani Secimi

Bu projedeki amacim olasi kirleten organizmalari (kontaminantlari) bulmak. Dolayisiyla genis bir veritabanina ihtiyacim var. Ancak veritabanini genis tutmak boyle bir avantaj sagliyorken, her dizi icin o veritabaninda arama yapmak oldukca fazla bilgisayar gucu ve zaman gerektiriyor. Bu yuzden projemi gelistirirken, cesitli veritabanlarini da inceliyorum. Ve ayrica bunlari nasil kisitlayarak, amacim icin en uygun hale getirebilecegimi arastiriyorum.

Ilk olarak NCBI'in Reference Sequence (Kaynak Dizi ya da Referans Sekans) -- RefSeq -- veritabaniyla basladim. RefSeq, cok genis, iyi adlandirilmis, plazmid organel, virus, arke, bakteri ve ökaryotlarin DNA, RNA ve protein dizilerini barindiriyor.1 Ben RefSeq'in genomik kismiyla calisiyorum ve o da gene oldukca genis bir veritabani. Komut satirinda yaptigim olcume gore 151 baz ciftinden olusan her okumanin RefSeq'te aranmasi ortalama 5 dakika aliyor ve CPU %50 oraninda kullaniliyor.  Elbette ilk aramada sunucu veritabanini sakladigi icin biraz daha fazla zaman geciyor (yaklasik 8 dakika) ancak birden fazla yapilan aramalarda ortalamada bu sayi dusuyor. Ayrica eger CPU'nun tamamini kullanabilirse, herbir okuma, tahmini 5 dakikadan daha az aliyor diyebiliyoruz. Eger boylece tum kaynaklari yeterli kullanabilirsek ve islemleri paralel calisacak sekilde tasarlarsak, bir hafta sonu gunu (tahmini sunucunun en az kullanilacagi zaman dilimi) yaklasik 2500 okuma aranabilir.

2500 okuma ne yazik ki toplam sayi olan 5.5 milyonu dusundugumuzde cok kucuk bir sayi. Bu yuzden okumalari paralel bicimde aratarak daha fazla hiz kazanmayi hedefliyorum. Yani, 1 milyonuncu okumadan baslayip 2500 okuma kadar arama yapmasini istiyorum ve bu islemi arkaplana atiyorum, daha sonra bir islem daha baslatiyorum ve bu sefere 1,002,500. okumadan baslayip 2500 okuma kadar aratiyorum. Boylece toplamda 8 paralel islem ile, isleri mumkun oldugunca hizlandirabiliyoruz.

Ancak gene de yetmiyor, hala bu genis veritabanina bir alternatif bulmak gerekli. Eger, yerine gecebilecek, ise yarar ama daha kucuk bir veritabani bulabilirsem, her sey daha cabuk yapilabilir.

RefSeq ile birlikte HUSAR tarafindan olusturulan honest ve Ensembl'in cDNA veritabani var. Simdilik dizileri bunlarda da arayip, karsilastirmalar yapiyorum.


1. The Reference Sequence (RefSeq) Database, http://www.ncbi.nlm.nih.gov/books/NBK21091/

FASTQ'dan FASTA'ya Donusturme Perl Scripti

FASTQ ve FASTA formatlari aslinda ayni bilgiyi iceren ancak birinde sadece herbir dizi icin iki satir daha az bilginin bulundugu dosya formatlari. Projemde onemli olan diger bir farklari ise FASTA formatinin direkt olarak MegaBLAST arama yapilabilmesi. Iste bu yuzden, genetik dizilim yapan makinelerin olusturdugu FASTQ formatini FASTA'ya cevirmem gerekiyor. Ve bu script pipeline'in ilk adimi.

Aslinda deneme amacli aldigim genetk dizilimin, bana bunu ulastiran tarafindan eslestirmesinin yapilmadigi icin, bir on adim olarak bu eslestirmeyi yapmistim. Ancak, pipeline'a boyle bir adim eklemeyecegim. Pipeline, icerisinde sadece eslesemeyen okumalarin (unmapped reads) bulundugu FASTQ dosyasiyla baslayacak.

Asagidaki detaylandirdigim script, create_fasta.pl adini verdigim Perl scripti.

Her zaman oldugu gibi scriptle su satirla basliyoruz:

#!/usr/local/bin/perl

Daha sonra kullanicidan scripti calistirirken istedigimiz iki parametreyi @ARGV dizisinden (array) okuyoruz. $input, FASTQ dosyasinin yolu, $n ise FASTA formatina cevirmek istedigimiz okuma sayisi. Eger tum dosyayi cevirmek istiyorsak, $n degeri girmeyebiliriz. Bu durumda asagidaki kontrol yapisindan anlayacaginiz gibi $n tanimlanmadiginda $n'i cok buyuk bir sayiya (1012) atiyor.

$input=$ARGV[0];
if(defined($n)){$n=$ARGV[1];}else{$n=1000000000000;}


Ardindan girdi dosyamizi (FASTQ dosyasi) aciyoruz.

open INPUT, $input or die $!;

Sonra $i, $j ve $k degiskenlerine sirasiyla 1, 2, 1 degerlerini atiyoruz. Bu degiskenlerden $i her okumanin baslik satirinin, satir numarasi, boylece her $i'inci satiri, baslik olarak gorup, o sekilde FASTA dosyamiza yazdirabiliriz. Ayni sekilde $j ise her okumanin dizilim satirinin satir numarasi, yani her $j'inci satirda bir diziyle karsilasiyoruz. Bu satirlar her 4 satirda bir karisimiza ciktigi icin (cunku FASTQ dosyasinda her okuma 4 satira sahip) bunlari daha sonra kontrol yapilarinin icinde 4'er artiracagiz. $k ise, sadece bir sayici. Her okumada, kontrol yapisi icinde, 1 artirilacak ve onu okuma adindaki okuma numarasini yazdirmak icin kullanacagim.

$i=1;
$j=2;
$k=1;


Simdi tum isi yapan while dongusu var. Bu dongu dosyayi satir satir okuyor ve satiri baslik ise ayri, dizi ise ayri degerlendirip ona gore cikti veriyor. Ayrica eger istedigimiz okuma sayisina ulasmissa donguyu sonlandiriyor.

while dongusunu asagidaki gibi kuruyoruz ve dosyamizdaki her okunan satiri $line degiskeninde sakliyoruz.

while(my $line = <INPUT>){

}

Daha sonra chomp komutu ile bu satirdan, satir sonundaki yeni satir karakterini siliyoruz.

chomp($line);

Son olarak asagida, kodun tumunde goreceginiz while dongusu icindeki kontrol yapilari ile de dosyamizi satir satir olusturuyoruz.


#!/usr/local/bin/perl
$input=$ARGV[0];
if(defined($n)){$n=$ARGV[1];}else{$n=1000000000000;}

open INPUT, $input or die $!;
$i=1;
$j=2;
$k=1;
                                                                                  
    while(my $line = <INPUT>){

        chomp($line);

        if ($i eq $.){
        print ">read_".$k." ".$line."\n";
        $i+=4;
        }
        elsif ($j eq $.){
        print $line."\n";
        $j+=4;
        $k++;
        }

        last if ($k>$n);
               
    }

close INPUT;


Simdilik script bu sekilde, neredeyse tamam ancak okuma isimlerini sabit bir sekilde "read_X" seklinde adlandirmaktansa, dinamik bir sekilde, girdi dosyasina gore belirlemek istiyorum. Bunu ileride degistirecegim.

Thursday, July 5, 2012

Eşleştirme ve Eşleşmeyen Okumaları Çıkarma Sonuçları

Daha önce verinin sadece bir kısmı ile çalışıyordum ancak artık tamamıyla çalışacağım. Bu yüzden bana sıkıştırılmış halde gelen veriyi direkt çalışma klasörüme çıkardım ve onun üzerinden işlemler yaptım.

Başlangıç (FASTQ) dosyamın boyutu 2153988289 bayt (2 GB). Ve bwa aracılığıyla eşleştirmeden sonra toplamda 6004193 dizilim, ya da okuma, (sequences ya da reads) ortaya çıktı. Daha sonra eşleşmeyen okumaları çıkarmam sonrasında toplam okuma sayısı 551065 kadar azaldı ve 5493128 oldu. Yani verinin %9.2'si eşleşeşen okumalara, kalan kısmı ise eşleşmeyen okumalara ait. Eğer her şey yolundaysa, bu sayı elimdeki eşleşmeyen okuma sayısı. Ve bu sadece eşleşmeyen okumaları içeren FASTQ dosyasının boyutu 1932192710 bayt (1.8 GB), yani orijinal dosyanın boyutu 221795579 bayt (0.2 GB) kadar azaldı.

Aslında bu kadar fazla eşleşemeyen okuma beklemiyordum. Elbette elimdeki veri böyle bir sonuç vermiş olabilir, çünkü "kötü" bir veri olduğu bana söylenmişti. Gene de kullandığım yöntemi deneyeceğim ve daha fazla araştırmayla bu sonucu yorumlamaya çalışacağım.

BWA İle Eşleştirme (Mapping - Alignment)

Bunu daha önce yazmayı unutmuşum. Aslında bahsetmiştim ancak nasıl yapıldığına dair bir şeyler yazmamışım ayrıca örnek komutlar da eklememişim.

BWA elimizdeki (FASTQ formatındaki) DNA dizilimini, referans genomunu (projemde bu insan genomu) alarak bir .sai dosyası oluşturuyor. Bu dosya dizinin ve referans genomunun eşleşmesi ile ilgili bilgiler taşiyor ve bu bilgileri kullanarak eşleşmeyenleri ayırabiliyorum.

İlk olarak aşağıdaki komut ile .sai dosyamızı oluşturuyoruz.

bwa aln $NGSDATAROOT/bwa/human_genome37 ChIP_NoIndex_L001_R1_complete_filtered.fastq > complete_alignment.sai

Oluşturduğumuz .sai dosyası çok da kullanışlı bir dosya değil, bu yüzden onu SAM dosyasına çevirerek, işlemlere devam ediyoruz. Elimizdeki veri tek-sonlu okuma (single-end read) olarak kabul edildiği için "samse" ile bu değişimi yapıyoruz, eğer çift-sonlu okuma (paired-end read) olsaydı "sampe" kullanılacaktı.

Bunu da aşağıdaki kod ile gerçekleştiriyoruz.

bwa samse $NGSDATAROOT/bwa/human_genome37 complete_alignment.sai ChIP_NoIndex_L001_R1_complete_filtered.fastq > complete_alignment.sam

Bundan sonra elde ettiğimiz SAM dosyasıyla çalışmamıza devam ediyoruz. Devamı aşağıdaki yazımda bahsettiğim BAM dosyasına dönüştürme ve ardından FASTQ dosyasına çevirme var.

İlgili yazı: SAM Dosyası - BAM Dosyası - samtools

Tuesday, July 3, 2012

SAM Dosyası - BAM Dosyası - samtools

Aslında programlamam gereken pipeline direkt olarak eşleşmeyen okumalar üzerinden analizler yapacak. Ancak böyle bir veri bulamadiığım için, elimdeki tek veri eşleşen ve eşleşmeyen okumaları içerdiği için önce  eşleşenlerden kurtulmam gerekti.

Bunu daha önce de belirttiğim gibi bwa eşleştiricisi (aligner - mapper) ile yapıyorum. bwa bir dizi işlemden sonra SAM dosyası oluşturuyor ancak benim FASTQ dosyasına ihtiyacım var. Bunun için SAM dosyasını samtools1 ile benzer bir format olan BAM dosyasına çevirip, daha sonra da bam2fastq2 aracı ile FASTQ dosyamı elde edeceğim.

SAM dosyasında, tahmin edebileceğiniz gibi sıralamalar (alignment) var. Ve her sıralama 11 satıra sahip ve bu satırlar sıralama, sıralayıcı (eşleştirici - aligner) ilgili bilgiler içeriyor. Bu satırlardan üçüncüsü sıralamanın ismi, onuncusu sıralamanin dizilimi (ya da sekansı) ve on birincisi ise sıralamalanın kalitesi. İşte bu alanlar özellikle FASTQ dosyası için gerekli, elbette diğer bilgiler de başka yerlerde öneme sahipler. BAM dosyası ise SAM dosyasının ikilik sistemde (binary) sürümü.

Aşağıdaki kod ile samtools komutunun view seçeneği üzerinden -S (girdisi SAM dosyası) ve -b (çıktısı BAM dosyası) parametreleriyle dosyanızı çevirebiliyorsunuz.

samtools view -Sb dosya.sam > dosya.bam

Daha sonra HUSAR paketinde bulunan, ücretsiz olarak elde edilebilecek bam2fastq komutu ile BAM dosyasını FASTQ dosyasına çevirerek işlemi tamamlıyoruz. Burada --output parametresi ile FASTQ dosyasının adını, --unaligned ve --no-aligned parametreleri ile de sadece eşleşmeyen okumaları almasını sağlıyoruz.

bam2fastq --output dosya_%#.fastq --unaligned --no-aligned dosya.bam

Sırada bu FASTQ dosyasını FASTA dosyasına çevirip, MegaBLAST ile aradıktan sonra çıktılar incelemek (parsing) var.

1. The Sequence Alignment/Map format and SAMtools, http://www.ncbi.nlm.nih.gov/pubmed/19505943
2. bam2fastq, http://www.hudsonalpha.org/gsl/software/bam2fastq.php

İlk Adım: Eşleşmeyen Okumaları Elde Etmek

Projemin ilk kısmı daha önce bahsettiğim gibi eşleşmeyen okumaları (unmapped reads) FASTQ dosyasından çıkarmak. Böylece, daha sonraki analizler için elimdeki ihtiyacım olmayan dizileri çıkarmış ve bu analizlerdeki iş yükünü azaltmış oluyorum.

Başından beri hedefim, tüm projeyi adım adım gerçekleştiren bir pipeline tasarlamak olduğu için bu işlemi bir Perl scripti ile yapacağım. Bu script pipeline'in ilk scripti ve laboratuvardan gelecek ham (raw) FASTQ formatındaki verinin girdi (input) olarak kullanılacağı yer. Aslında bu scripte ihtiyacım olmayacak, sadece elimdeki verinin eşlenebilen verileri de içermesi sebebiyle bu adımı ekledim.

Scripti yazdığım platform bir Unix bilgisayarı, bu yüzden pipeline komut satırından birkaç parametre ile birlikte çalıştırılıyor olacak, yani bazı parametreleri komut satırından gönderecegim. Bu parametreler, okuma türü, FASTQ dosyası ya da dosyaları.

Okumalar iki tür olabiliyor: tek-sonlu okuma (single-end read), çift-sonlu okuma (paired-end read). Tek-sonlu okuma DNA'nın sadece bir sonundan başlanarak tamamının dizilenmesi demek ve çift-sonlu okuma ise DNA'nın her iki sonundan da dizilemenin gerçekleştirilmesi anlamına geliyor. Bu yüzden tek-sonlu okuma tek bir FASTQ dosyasına sahipken, çift-sonlu okumada ise iki FASTQ dosyasıyla çalışıyorum. İşte bu yüzden kullanıcıdan (tek-sonlu okuma için) "se" ya da (çift-sonlu okuma için) "pe" olarak iki türden birini ilk parametre olarak belirtmesini istiyorum. Ardından da kullanıcıdan FASTQ dosyalarının adreslerini eklemesini bekliyorum. Elbette bunu kontrol etmenin birçok yolu var. Şimdilik böyle başladım, ileride daha da kolaylaştırıcı şekilde değiştirebilirim.

Kullanıcı bu parametrelerde Perl scriptini çalıştırdıktan sonra, program aldığı parametrelerle bir dizi işlem sonunda sadece eşleşmeyen okumalariı içeren FASTQ dosyasını ekrana standart çıktı (standard output) olarak yazdırıyor. Aşağıda örnek komutu görüyorsunuz.

perl extract_unmapped_reads.pl pe hs_sequencing_1.fastq hs_sequencing_2.fastq

İşte böylece Perl scripti çalıştırılmış oluyor.

Bir sonraki yazımda scriptin ayrıntılarına, bwa, samtools gibi kullandığım diğer komutlara değineceğim.

Monday, July 2, 2012

Blog Yazılarını Facebook Twitter ve LinkedIn'e Yönlendirmek

İlgilendiğim bir konu üzerine bir blog açıp, bilgilendirici yazılar yazmak uzun süredir aklımda olan bir şeydi. Sonunda ufak ufak yazılarıma başladım. Umarım şu ana kadar güzel gitmiştir.

Bu yazıda blog başlığıyla çok alakalı olmayan "konu-dışı" bir konudan bahsedeceğim.

Yazılarımı kolay bir şekilde geniş kitleye ulaştırmak için sosyal medyayı kullanmak istiyordum ama her seferinde yazının bağlantısını kopyala-yapıştır yapmak hiç de basit bir iş değil.

Aramalarım sonunda bunu, blogumuzu Facebook, Twitter ve LinkedIn hesaplarımıza bağlayarak aynı anda yeni yazıları yönlendirebildiğimiz bir araç buldum. Bit.ly adresine sahip hepimizin tanıdığı URL kısaltma servisi twitterfeed adında bir başka servisle bunu yapıyor.

Şunu ek olarak söylemeliyim ki, eğer yazınızı Facebook profilinize değil de sadece Facebook sayfa(lar)ınıza yönlendirmek istiyorsanız da, bu aracı kullanabilirsiniz.

Twitterfeed'e ister bilgilerinizi girerek, ister Open-ID ile Blogger, Facebook; Twitter gibi hesaplarınız aracılığıyla kayıt olduktan sonra oldukça basit kullanıcı arayüzünden blog akışını (feed) girerek ve sonra da istediğiniz sosyal medyaları seçerek hemen yazılarınızı yönlendirmeye başlayabilirsiniz.

Başka yöntemler de var, ancak bu en basiti olarak gözüme çarptı. Eğer ileride daha iyisi bulursam, mutlaka onu da yazarım.

Düzenleme: Twitter hesabımda bu yönlendirme (ve ya aynı zamanda Google Feedburner'ı da denediğim için onunla) ilgili bir sorun yaşadım. Hesabımı askıya aldılar ve birkaç e-posta sonunda tekrar kullanıma açtılar. Her iki servisi de deniyor olduğum için emin değilim ama biri sorun yaratıyor. Belki de aynı anda çalıştıkları için bu oldu. Henüz emin değilim ama Twitter'da böyle bi durum var. Facebook ve LinkedIn'de sorun yok, çalışmaya devam ediyor.

Thursday, June 28, 2012

MegaBLAST - Dizilerdeki Benzerlikleri Bulma Aracı

MegaBLAST, HUSAR paketinde bulunan, BLAST (Basic Local Alignment Search Tool) paketinin bir parçası. Ayrıca BLASTN'in bir değişik türü. MegaBLAST uzun dizileri BLASTN'den daha etkili bir şekilde işliyor ve hem de çok daha hızlı işlem yapiyor ancak daha az duyarlı. Bu yüzden benzer dizileri geniş veri tabanlarında aramaya çok uygun bir araç.

Yazacağım program çoklu dizilim barındıran FASTA dosyasını alacak ve megablast komutunu çalıştıracak. Daha sonra da her okuma için bir .megablast dosyası oluşturarak bir sonraki adıma hazırlayacak.

Wednesday, June 27, 2012

Kontaminant (Kirletici) Analizi Projesi

Başlangıç olarak, araçlara, programlama diline, kısacası biyoenformatiğe alışabilmem için bana verilen bu ufak projeyi ayrıntılı olarak anlatacağım.

Biliyoruz ki, laboratuvar çalışmalarımızda ne kadar önlemeye çalışsak da kontaminant riski hep bulunuyor. Bunu ne kadar aza indirsek o kadar iyi, ki daha sonra bunun miktarını bulup, bunun üzerinden sonucumuzun bir başka değerlendirmesini de yapabiliriz. İşte bunu bulmak için bir yöntem, DNA analizi. Çalıştığınız örneğinizin DNA'sı dizileniyor ve bu DNA çeşitli programlarla analiz edilip, kirleten organizmaları DNA'larından ortaya çıkarabiliyoruz.

Literatür araştırması sonrası bu tarz programlar olduğunu gördüm ancak daha farklı açıdan yaklaşmışlar. Örneğin Tel Aviv Üniversitesi'nden bir grup, benzer bir çalışmayı insanda patojen enfeksiyonlarını bulma ve buradan tedavi geliştirme konusunda gerçekleştirmişler.1

Bu proje birkaç Perl scriptinin oluşturulmasını içeriyor. Bunların ilki, DNA diziliminin bir referans genomuyla karşılaştırılarak, eşleşebilen kısmının çıkarılması ve kalan kısmın FASTQ formatında kaydedilmesi. Daha sonra bu dosya, sonraki adımlar için FASTA formatına bir başka Perl scripti ile dönüştürülecek. Ardından Bu FASTA dosyası, MegaBLAST aracı ile taranacak ve sonuçlar kaydedilecek. Ve şimdilik son olarak da bu sonuçlar değerlendirilecek.

Proje yapım aşamasında değişebiliyor, ufak eklentiler gelebiliyor. Şu an için bu şekilde planlandı. Projeyle ilgili diğer ayrıntıları, devam ederken diğer yazılarda vereceğim. Özellikle her adımlar neleri uyguladığımı yazmak ve böylece benzer işleri yapacak arkadaşlarım için bir kaynak oluşturmayı düşünüyorum.

1. Pathogen detection using short-RNA deep sequencing subtraction and assembly, http://bioinformatics.oxfordjournals.org/content/27/15/2027.full

Monday, June 25, 2012

FASTQ Formatı - FASTQ Dosyası

Bugün programı oluştururken kullanacağım "test" dizilimini aldım. İki adet FASTQ dosyasından oluşuyor, her biri sıkıştırılmış ama buna rağmen boyutları 6 GB civarı. Ben elbette çok zaman kaybetmek istemediğim için bu dosyalardan birinin sadece bir kısmını kullanacağım.

Amacım, bu FASTQ dosyalarındaki eşleşebilen okumaları BWA aracı ile bularak, daha sonra onları çıkarmak. Ve kalan eşleşemeyen okumaları MegaBLAST aracının anlayabileceği bir dilde (FASTA formatında) kaydetmek.

Bu arada tüm projeyi bir Unix bilgisayarda hazırladığım için birçok komut öğreniyorum, daha sonra bunları ayrıca yazmaya çalışacağım. Örneğin sıkıştırılan bir dosyayı gzip -d dosya komutuyla çıkarıyorum, eğer özgün dosyayı korumak istiyorsam komuta -c ekliyorum.

Ayrıca Perl scriptlerini Emacs metin editöründe hazırlayacağım. Hiç alışkın olmadığım kısayollara ve özelliklere sahip olduğu için şimdilik zorlanıyorum. Tamamen öğrendikten sonra Emacs ile ilgili yazılar da yazmayı planlıyorum.

FASTQ Formatı


Bir FASTQ dosyasında her dizilim için (diyelim ki bir dizilim 151 bp) dört adet satır bulunuyor ve bu dört özel satır alt alta bu dosyada her dizilim için sıralanıyor. Bu satırların ilki başlık (header) ve dizilim ile ilgili özel belirteç (identifier) bulunuyor ve @ karakteri ile başlıyor. İkinci satır doğrudan baz dizilimini gösteriyor. Üçüncü satırda açıklamalar bulunuyor, boşsa sadece bir + karakteri görüyoruz. Son satır ise dizilimin kalitesini belirten kodlara sahip.

Thursday, June 21, 2012

BWA (Burrows-Wheeler Aligner) Hizalayıcı - Eşleştirici

Önceki yazımda belirttiğim gibi bir eşleştirici (aligner ya da mapper) kullanarak elimdeki verinin referans genomu ile ne derece eşlestiğini bulmaya çalışacağım. Daha sonra eşleşmeyen kısmıyla birtakım analizler yapacağım.

BWA (Burrows-Wheeler Aligner) görece kısa dizilimleri insan genomu gibi uzun referans genomlarıyla eşleştiren bir program. 200bp (bp: baz çifti) uzunluğuna kadar bwa-short algoritması, 200bp - 100kbp arası ise BWA-SW algoritması kullanılıyor.

Hizalayıcı - eşleştirici seçmede birçok faktör rol oynuyor. Birçok bu tip araç var ve farklı özelliklere sahipler. Bunlar performans ve hassasiyet (accuracy) özelliklerine göre değişiyor. Ayrıca hassasiyet, yapısal varyasyonları (structural variations - SVs) da etkiliyor.

İleride bu programın özelliklerine daha fazla değineceğim.

Wednesday, June 20, 2012

Dizileme Çalışmalarını Kirleten Organizmaları Tespit Etme

Bu yaz stajımda ilk olarak başlayacağım çalışma yavaş yavaş şekilleniyor. Bu çalışmada bir pipeline oluşturup, bunu laboratuvarlarda dizileme (sequencing) örneklerini kirleten organizmaları bulmaya çalışacağım.

Laboratuvarlarda birçok nedenden dolayı örnekler başka organizmalar ya da yabancı DNA tarafından kirlenebiliyor. Bunlar bakteri, maya olabilir ya da bir virüs DNA'sı da olabilir. Siz bir DNA'yı diziledikten sonra onun referansıyla eşleştirme çok az oranda çıkabiliyor. Bu da yabancı DNA'nın olabileceğini gösteriyor. Bir başka neden referans DNA'nın farklı olması da olabilir. Aynı zamanda algoritmada da hata olabilir. Bunlar düşük oranda eşleşme çıkaran sonuçların nedenleri olarak gösterilebilir.

Kirletici bir faktör olup olmadığını tespit etmenin çeşitli stratejileri var. Örneğin biri, genomu diziledikten ve referansıyla eşleştirdikten sonra, eşleşmeyen kısmını bu genomdan çıkarıp, onu olası tüm organizmalarla deneyerek, en fazla eşleşmeyi bulmaya çalışmak ve organizmayı belirlemek. Bunu Husar paketinde bulunan MegaBLAST ile yapabiliyoruz. Veri tabanında sahip olduğu 20.000 canlı DNA'sı üzerinden karşılaştırma yaparak, kirletici organizmaları listeleyebiliyor.

İşte yapacağım çalışmada, pipeline dizilenen genomu alacak, referansıyla eşleştirmeye çalışacak, eşleşmeyen kısmını çıkarıp bunu MegaBLAST ile diğer genomlarla karşılaştırıp bana olası organizmaların listesini döndürecek. Oldukça kolay görünen ancak verilerin çok çok büyük boyutlara sahip olması sebebiyle uzun zaman alacak bir çalışma.

Başlamadan önce Perl programlama diline ve MegaBLAST aracına alışmam gerekiyor. Elbette literatür araştırması yaparak, başka stratejiler geliştirilip geliştirilmediğini kontrol etmekte de fayda var.

Pipeline ve Pipeline Geliştirme


Bugün aldığım tanıtım derslerinin devamında, pipeline ve pipeline geliştirme ile ilgili ayrıntılı bilgiler aldım. Pipeline, aslında bildiğimiz boru hattı demek, örneğin borularla petrolün bir yerden başka bir yere taşınması için kullanılan sistem. Bunun bilgisayar terminolojisinde anlamı ise bir elementin çıktısı, diğerinin girdisi olacak şekilde oluşturulmuş işleme elementleri zinciri. Böylece çok daha komplike işlemler pipeline oluşturularak, kolay ve düzenli bir biçimde gerçekleştiriliyor. Sanırım pipeline Türkçeye ardışık düzen olarak çevriliyor, gene de ben pipeline olarak kullanacağım.

Pipeline oluşturulması birkaç adım içeriyor. Öncelikle işlem öncesi (pre-processing) yapılması gereken kontroller var, bunlar örneğin JavaScript gibi client-side bir dil ile yapılabiliyor. Amaç, kullanıcının girdilerini kontrol etmek ve henüz işlemeden herhangi bir olası uyumsuzluğu, hatayı önceden bildirmek ve kullanıcıyı uyarmak. İşlem süresince ise pipeline, gelen verinin analizini yapıyor, eğer ileride tekrar kullanmak ıstıyorsak saklıyor ve bunu çözümledikten sonra çıktı olarak veriyor. Son, çıktı adımında ise en iyi gösterim için stil dosyaları kullanılıyor ve çıktı XML dosyasına dönüştürülerek, daha sonra kolayca kullanılabilecek bir formatta gösteriliyor.

XML'e dönüştürme gerçekten çok önemli, çünkü böylece hem makina, hem de insanın okuyabileceği bir formatta sonuçlarınızı sunabiliyorsunuz. XML (Genişletilmiş İşaretleme Dili) verileri sizin belirleyebileceğiniz etiketler ve özelliklerle saklamanızı sağlayan bir standart.

Tuesday, June 19, 2012

WWW2HUSAR - HUSAR'ın Web Arayüzü

Stajımın ikinci gününde HUSAR'ın web arayüzünü konuştuk. HUSAR komut isteminden komutlarla kullanılabilen, yönetilebilen bir yazılım ancak bunu kolaylaştırmak için hazırlanmış bir web arayüzü var. WWW2HUSAR adını verdikleri bu arayüz ile listelenen araçları kolayca seçebiliyor, genetik dizinizi ekleyebiliyor ve başka birçok işlemi kolayca, birkaç tık ile yapabiliyorsunuz.

Bununla birlikte biraz daha HUSAR'ın işlevlerine göz attık. Yazılımda, yerel klasörde gen dizisi listeleri oluşturarak, bunları çoklu dizi hizalama (multiple sequence alignment) aracı ile genlerin benzerliklerini karşılastırabiliyor ve örneğin evrimsel ilişkilerini ortaya çıkarabiliyorsunuz.

Monday, June 18, 2012

DKFZ - Heidelberg Biyoenformatik Birimi'nde Staj

Erasmus programıyla yapıyor olduğum yaz stajı başladı. İlk olarak birimi yöneten bilim insanlarından birkaç saatlik tanıtım dersi aldım. Bu derste birimin kısa tarihi, birimin günümüze kadar yaptıkları projeler ve bunlarin ayrintilari konusunda bilgiler aldım.

Biyoenformatik Birimi DKFZ'nin (Deutsches Krebsforschungszentrum -- ing. German Cancer Research Center) bir çekirdek tesisi olan Genomik ve Proteomik Çekirdek Tesisi'ne bağlı bir grup. İsimleri aynı zamanda HUSAR (Heidelberg Unix Sequence Analysis Resources) ve bu isim grubun geliştirdiği dizi analizi yapma paketinin de adı olarak kullanılıyor. HUSAR hem komut satırından hem de web üzerinden yönetilebiliyor.

Adından da anlaşılacağı gibi HUSAR Unix sisteminde çalışan bir yazılım paketi. Bu yüzden her ne kadar web arayüzü olsa da üzerinde çalışmak için Unix'i öğrenmem gerekecek. Daha önce hiçbir Unix işletim sistemi kullanmadığım için ilk gün labda bulunan CentOS dağıtıcısının işletim sistemini anlamaya keşfetmeye çalıştım.

Son olarak ilk günde SRS (Sequence Retrival System) hakkında konuştuk. SRS, HUSAR'in kullandığı hızlı, kolay ve kullanıcı dostu bir data alma sistemi. Bu sistem 90'dan fazla herkese açık veri tabanına bağlanabiliyor. Böylece farkli veri tabanlarını kullanarak sorgular yapılabiliyor.

Saturday, March 17, 2012

Biyoinformatik mi? Yoksa Biyoenformatik mi?

Yazılarıma konu ararken kitaplarla birlikte interneti de karıştırıyorum. Yabancı kaynaklar elbette fazlaca var ve yeterliler, ancak Türkçe kaynaklara baktığımda ilk gözüme çarpan bu alanın isminin farklı kullanımları oldu.

Biliyorsunuz, İngilizcede bu alana bioinformatics deniyor. Gayet normal, çünkü İngilizcede informatics "ics" eki ile birlikte information sözcüğünden geliyor. Bu sözcük ise Latince kökene sahip1. Enformatik sözcüğü Türkçeye, Fransızcadan informatique sözcüğünden, "enformatik" olarak gelmiş, ayrıca bilişim olarak da Türkçesi önerilmiş2. Elbette bu Fransızca sözcük de İngilizcesi ile aynı kökene sahip.

Bu birkaç etimolojik bilgiden sonra, biyoinformatik sözcüğünün hatalı bir aktarma olduğunu görebiliriz. Doğru kullanım "biyo" ön eki ile enformatik kelimesinin birleşiminden oluşan biyoenformatik sözcüğü olacaktır.

İlginç bir not ise, Google'ın biyoenformatik aramasını "biyoinformatik mi demek istediniz" diyerek düzeltmesi tüm bunlardan şüphe duymamıza neden olabiliyor. Umarım Google da bu durumu düzeltecek bir adım atacak.

Peki, neden TDK'nin önerisi bilişim sözcüğünü kullanmıyoruz? Ne yazık ki, bunu tartışacak kadar dil bilgim yok ve alanı da bu kadar fazla tanımıyorum. Ama bu, bilişim kelimesinin zaten henüz çok fazla kullanımda olmamasından dolayı olabilir.

Kaynaklar
1. Wiktionary: information, http://en.wiktionary.org/wiki/information
2. TDK Güncel Türkçe Sözcük, http://www.tdk.gov.tr

7th International Symposium on Health Informatics and Bioinformatics

7. Sağlık Enformatiği ve Biyoenformatik üzerine Uluslararası Sempozyumu, 7th International Symposium on Health Informatics and Bioinformatics (HIBIT 2012), ilk kez 2005'te ODTÜ Enformatik Enstitüsü tarafından düzenlenmiş ve Sağlık Enformatiği, Tıbbi Enformatik, Hesaplamalı Biyoloji ve Biyoenformatik alanlarında akademisyenleri ve araştırmacıları bir araya getirmeyi ve bu alanlar hakkında yapılan çalışmaların sunulmasına ortam sağlamayı ve çalışmalar üzerine interaktif bir şekilde değerlendirmeler yapmayı amaçlamaktadır.

Bu sene, 19-22 Nisan 2012'de Ürgüp, Nevşehir Perissia Hotel'de düzenlenecek olan HIBIT 2012 organizasyonu ODTÜ, ODTÜ Enformatik Enstitüsü, ODTÜ Biyolojik Bilimler Bölümü ve ODTÜ Bilgisayar Mühendisliği Bölümü partnerliği ile gerçekleştirilmektedir.

Bu sempozyum ile ilgili tüm ayrıntıları, programı, önemli tarihleri, kayıt sayfasını, resmi web sitesinden öğrenebilirsiniz.

Friday, March 16, 2012

Hoş Geldim! Hoş Geldiniz!

Merhabalar,

Biyoloji alanında özel olarak ilgi alanım olan ve daha fazla keşfetmem, üzerine çok şey öğrenmem gereken Biyoenformatik'i, bu blog aracılığıyla (olası ziyaretçilerimle birlikte) öğreneceğim. İlk yazımı biraz önce Biyoenformatik'in çeşitli otoriteler tarafından yapılan tanımları ile tamamladım. Daha sonra, Biyoenformatik'te geçen birçok ilkelerin tanımlarından da bahsetmek istiyorum. Ayrıca, Biyoenformatik hakkında yazılım dilleri, istatiksel yöntemler de yazılarımın konularını oluşturacak. Aynı zamanda Biyoenformatik ile ilgili haberlere de yer vermek ve bu haberlerle en son gelişmeleri takip etmeyi (ettirmeyi) planlıyorum.

Umuyorum, yeterli zamanı ve gücü bulacağım ve elimden geldiğince dolu dolu, güncel bir blog yaratabileceğim. Ve umarım bu blog, Biyoenformatik ile ilgilenen diğer arkadaşlarım için de yol gösterici olacaktır.

Giriş yazım da bu olsun. :)

Tekrar görüşmek üzere!

Biyoenformatik Nedir? Biyoenformatik'in Tanımı

Birçok organizmanın ve son olarak da 2001'de insan genomunun çıkarılmasıyla, tüm 3 milyar baz çiftinin diziliminin elde edilmesiyle, karşımıza bu bilgiyi farklı şekillerde kullanacak olan alanlar çıktı.Bu genleri anlamaya çalışan, bu genlerden oluşacak proteinleri belirlemeye çalışan alanların yanında bu bilginin analizini yapma ihtiyacı da Biyoenformatik alanını doğurdu.

Biyoenformatik, biyolojik bilginin bilgisayarlar ve istatistiksel teknikler kullanılarak analiz edilmesidir; başka bir deyişle, biyoenformatik, biyolojik araştırmaları iyileştirmek ve hızlandırmak için bilgisayar veri tabanları ve algoritmaları geliştirme ve onlardan yarar sağlama bilimidir1.

NCBI (National Center for Biotechnology Information) Biyoenformatik'i şöyle tanımlamıştır: "Biyoenformatik, Biyoloji, Bilgisayar Bilimi ve Bilişim Teknolojisi bilimlerinin tek bir disiplinde birleşme alanıdır. Bu alanın nihai hedefi, yeni biyolojik keşiflere olanak sağlamak ve Biyoloji'deki birleştirici ilkelerin ayırt edilebileceği evrensel bir bakış açısı yaratmaktır. Biyoenformatik'te üç tane alt-disiplin vardır; geniş veri setleri arasındaki ilişkileri belirleyecek yeni algoritmaların ve istatiklerin geliştirilmesi, nükleotid, amino asit dizilimleri, protein domainleri ve protein yapıları gibi çeşitli türde verilerin analizi ve yorumlanması, farklı türde bilginin yönetimine ve bunlara etkili erişimine olanak sağlayacan araçların geliştirilmesi ve uygunlanması bu alt-disiplinlerdir."2.

NIH (National Institute of Health) Biomedical Information Science and Technology Initiative Consortium'a göre ise Biyoenformatik'in tanımı biyolojik, tıbbi, davranışsal ya da sağlık verilerinin, ve bu verilerin alınması, saklanması, organize edilmesi, arşivlenmesi, analiz edilmesi ve görselleştirilmesini içeren, bilgisayar araçlarının ve yaklaşımlarının araştırılması, geliştirilmesi ya da uygulanmasıdır3.

Kaynaklar
1. Definition of Bioinformatics
http://www.roseindia.net/bioinformatics/definition_of_bioinformatics.shtml
2. Bioinformatics
http://www.ncbi.nlm.nih.gov/Class/MLACourse/Modules/MolBioReview/bioinformatics.html
3. NIH WORKING DEFINITION OF BIOINFORMATICS AND
COMPUTATIONAL BIOLOGY
http://www.bisti.nih.gov/docs/CompuBioDef.pdf