Güngör Budak's Blog

Bioinformatics, web programming, coding in general

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.

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.

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.

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.

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.

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);
}