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.