Güngör Budak's Blog

Bioinformatics, web programming, coding in general

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.

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.

Kaynaklar

  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.

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

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.

Kaynaklar

  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.

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.

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.

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.

Kaynaklar

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