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.