6 Mart 2017 Pazartesi

Yeni Nesil Sekanslama Temel Veri Analizi - 3

Geçen yazıda elimizdeki ham verilerdeki adaptörleri temizleyip gerçek okumaları elde etmeyi başarmıştık, fakat ortalığı da biraz dağıtmıştık. Sadece ara aşamalarda kullandığımız ve son aşamada artık ihtiyaç duymadığımız dosyaları ya silebiliriz, ya da bu süreçli ardarda bağlayarak her bir aşamadan elde ettiğimiz çıktıyı doğrudan sonraki aşamaya gönderebiliriz. UNIX dünyasında bu yaklaşıma Pipeline adı veriliyor, Türkçe'ye boru programlama veya boru haberleşme şeklinde çevrilmiş. Bilgisayarda kayıtlı dosyaları birer depo olarak düşünelim, depolardan sıvı çekip bu sıvıları işleyen veya diğer bir deyişle dosyalardan veri çekip bu veriler üzerinden işlem yapan programları da birer işleme tesisi olarak hayal edelim. Bu benzetme üzerinden gidelim: şimdiye kadar yaptığımız şey FASTQ depolarından bilgiyi satır satır alıp işleyen programlar kullanıp, işlenmiş verileri tekrar depolara, veya dosyalara kaydetmekti. Artık aradaki depolardan kurtulup işleme tesisleri doğrudan birbirine bağlayabiliriz: Böylece hem aradaki dosyalardan kurtulur, hem de yaptığımız işlemleri daha derli toplu olarak kaydedebiliriz.

Önce basit bir örnekle başlayalım. Üzerinde çalıştığımız dosyaların ortak özelliği SRR1513652 ile başlamaları. İlk aşamada Terminal'i kullanarak bunları ls ile listeleyelim:

ls SRR1513652*

Bilgisayardan aldığımız yanıt:

SRR1513652.fastq.gz SRR1513652.paste SRR1513652.seq SRR1513652_trimmed.fastq

Yukarıda kullandığımız * karakteri, o pozisyona herhangi bir karakter veya karakter grubunun gelebileceğini ifade ediyor. Tek başına * karakterini kullanarak ls * gibi bir komut kullansaydık bulunduğumuz klasördeki tüm dosyaları listelerdik. Şimdi hem dosyaları listeleyip hem de dosya büyüklüklerini görüntüleyelim:

ls -lh SRR1513652*

Şöyle bir yanıt bekliyoruz:

-rw-r--r--  1 ahmetrasit  staff    64M Feb 22 16:13 SRR1513652.fastq.gz
-rw-r--r--@ 1 ahmetrasit  staff   280M Feb 22 16:50 SRR1513652.paste
-rw-r--r--@ 1 ahmetrasit  staff    64M Feb 22 16:51 SRR1513652.seq
-rw-r--r--@ 1 ahmetrasit  staff   280M Feb 22 16:21 SRR1513652_trimmed.fastq

Elde ettiğimiz cevapta kaç satır olduğunu nasıl öğrenirdik? Tabii ki wc kullanarak. Bu durumda, ls isimli işleme tesisinden elde edilen çıktıyı doğrudan wc isimli işleme tesisine gönderebilmenin bir yolunu bulsak aradığımız cevaba da ulaşırız. İşte tam da bu noktada konumuza geri dönüyoruz: | karakterine. Nasıl ki bir komut sonucunda elde edilen verileri bir dosyaya yazdırmak için > karakterini kullanıyorsak, bir komuttan elde ettiğimiz sonucu bir diğerine göndermek için de | karakterini kullanıyoruz. Hadi deneyelim:

ls -lh SRR1513652* | wc

Elde ettiğim yanıt:

       4      36     285

İşte bu! ls komutu sonucunda elde ettiğimiz çıktıda 4 satır, 36 kelime, ve 285 karakter yer alıyormuş. Peki bir de şunu deneyelim: ls ile elde ettiğimiz yanıtta sadece dosya büyüklüklerini ve dosya isimlerini görüntüleyelim. Burada devreye awk giriyor: ls'ten elde ettiğimiz sonucu doğrudan awk'a yönlendirip burada da 5. ve 9. sütunları ekrana yazdıracağız:

ls -lh SRR1513652* | awk '{print $5,$9}'

Ve elde ettiğimiz sonuç:

64M SRR1513652.fastq.gz
280M SRR1513652.paste
64M SRR1513652.seq
280M SRR1513652_trimmed.fastq

Görüntülediğimiz bilgilerin sırasını değiştirelim ve önce dosya adını görüntüleyelim:

ls -lh SRR1513652* | awk '{print $9,$5}'

Ve yanıt:

SRR1513652.fastq.gz 64M
SRR1513652.paste 280M
SRR1513652.seq 64M
SRR1513652_trimmed.fastq 280M

Artık bu saatten sonra bilgisayara yaptıramayacağımız şey yok :) Geçen yazıdaki adımları artık tek satırda daha derli toplu bir şekilde yazabiliriz:

cutadapt -a CTGTAGGCACCATCAAT -e 0.1 -O 5 -m 15  SRR1513652.fastq.gz | paste - - - - | awk -F\t '{print $2}' > SRR1513652.seq

Bitti :) Bu kadar. Hadi özetleyelim: cutadapt programını çalıştırıp, sonucu bir dosyaya yazmamasını ve onun yerine ekrana yazdırmasını -o ve dosya adını komuttan kaldırarak söyledik, bu durumda temizlenen okumaları ekrana yazdırmaya karar verdi, biz de hemen bu çıktıyı alıp paste komutuna yönlendirdik ve aradaki depoyu ortadan kaldırdık. paste komutu kendisine gelen verileri işleyip her dört satırı bir araya topladı ve bir satır haline getirdi, tam bunları ekrana yazdıracakken o çıktıyı da toplayıp awk komutuna yönlendirdik. awk kendisine bir dosya adı verilmediği için kendisine yönlendirilen çıktıyı aldı ve ikinci sütunu, yani sekansların yer aldığı sütunu seçip onu bize verdi. Artık bu son işlem olduğu için bu çıktıyı alıp bir dosyaya yönlendirdik ve kaydettik. Yine de her şey yolunda mı diye dosyamızın baş kısmına bir göz atalım:

head SRR1513652.seq 

Her şey yolunda:

ACTGGTCTTCATTAATCAAGAACGAAAGTCA
ACTGCACCGTGACTCTCCTCGAAGACAGTGTCAAGCGG
ACTGTGTTCAGTTTATCGTACAAAT
ACTGTTCTTAGTTGGTGGAGTGATTTGTCTGGTTT
ACTGTATGAGATAGTGACTAAAAATATAAA
ACTGTTTTATTTGGAGCATGATCAA
ACTGATTCCACTTGGGGTAAATTC
ACTGGACCATTTCGTCCAGAGTTTTTGGGT
ACTGTCAACTGGAAATATTTATCAA
ACTGCACCTTATCCGGGATCCTCATATGG

Evet. Hadi bir de oluşturduğumuz dosyada kaç okuma varmış onu da görelim:

wc SRR1513652.seq 

Yanıtımız:

2173777 2173777 67118919 SRR1513652.seq

İşte bu kadar, her şey yolunda, daha önce adım adım gerçekleştirdiğimiz işlemleri bir araya getirdiğimizde de yine sorunsuz bir şekilde çalıştırmayı başardık. Yani aslında bir önceki yazımızda prototipini geliştirdiğimiz analiz sürecini üretimde kullanılır hale getirdik. Bundan sonra tek yapmamız gereken, farklı bir dosyayla karşılaştığımızda dosya adını değiştirmek, o kadar.

Son olarak, elimizdeki okumaları hızlıca sayıp en fazla hangi sekansla karşılaştığımızı hızlı bir şekilde görüntüleyelim. Bunun için de iki yeni komut kullanacağız ve elde ettiğimiz sonucun tamamını ekrana yazdırmak yerine en sık karşılaşılan ilk on sekansı yazdıracağız:

sort SRR1513652.seq | uniq -c | sort -k1nr | head

Elde ettiğimiz yanıt:

44385 ACTGTCACCGGGTGTACATCAGCTAA
34004 ACTGTCACCGGGAGAAAAACTGGAGT
33222 ACTGTCACCGGGTGAACACTTGCAGT
32559 ACTGTCACCGGGTGAAAATTCGCATG
15257 ACTGTCACCGGGTGGAAACTAGCAGT
10642 ACTGGAGACTATCCTTTGATTGAGTTTTTTGT
9775 ACTGGACTATCCTTTGATTGAGTTTTTTGT
9282 ACTGCACCCGTACATATGTTTCCGTGCT
8789 ACTGAGACTATCCTTTGATTGAGTTTTTTGT
8175 ACTGACTATCCTTTGATTGAGTTTTTTGT

Yukarıdaki komut grubunu açıklayıp elde ettiğimiz sonuçtaki ilk 10 sekansın neler olduğunu keşfetmeyi size bırakacağım. Bu sekansların C. elegans'a ait olduğunu da hatırlatayım. sort komutu, kendisine sunulan verileri alfabetik sıraya göre sıralar. Bu ilk komutumuz olduğu için hangi dosyayı sıralaması gerektiğini sort komutundan hemen sonra yazarak belirttik. Elde ettiğimiz sıralanmış sonuçları uniq komutuna gönderdik: Bu komut kendisine gelen veride birbirinin aynısı olan satırları çıkarıp sadece birer kez sunar, -c parametresiyle kullanıldığında da her bir satırla kaç kez karşılaşıldığını belirtir. Elde ettiğimiz veriyi tekrar sıralamalıyız, bunun için tekrar sort komutuna başvuracağız. Önce uniq -c komutunun bize  normalde nasıl bir çıktı verdiğine bir göz atalım:

   1 ACTGAAAAAAAAAAAAAAAA
   2 ACTGAAAAAAAAAAAAAAAAA
   1 ACTGAAAAAAAAAAAAAAAAAA
   3 ACTGAAAAAAAAAAAAAAAAAAA
   8 ACTGAAAAAAAAAAAAAAAAAAAA
   4 ACTGAAAAAAAAAAAAAAAAAAAAA
   2 ACTGAAAAAAAAAAAAAAAAAAAAAA
   1 ACTGAAAAAAAAAAAAAAAAAAAAAAA
   1 ACTGAAAAAAAAAAAAAAAAAAAAAAAAA
   1 ACTGAAAAAAAAAAAAAAAAAAAAAAAAACTGTAGGCCCCCCCAAATTGT

Bu veriyi sort komutuyla tekrar sıralarken sort komutuna bir öncekine göre daha fazla bilgi verip tam olarak ne istediğimizi anlatmalıyız. Bunu da -k1nr parametreleri ile sağlıyoruz. sort komutu ile birlikte kullandığımız -k, veriyi boşluk kullanarak sütunlara ayırmak istediğimizi, ardından gelen 1 ile birinci sütunu kullanarak sıralama yapmak istediğimizi, ardından gelen n karakteri ile bu sütunun sayılara dönüştürülmesini (bu durum biraz garip gelebilir size, farkı görmek için n karakterini çıkarıp yukarıdaki komut grubunu bir de bu şekilde çalıştırırsanız farkı göreceksiniz), ve son olarak eklediğimiz r karakteri ile de sonuçları görüntülerken çoktan aza doğru sıralayarak görüntülemesini istediğimizi belirtiyor.

Nasıl, iyi gidiyoruz değil mi? Şimdiye kadar bahsettiğimiz şeyleri sindirmek zaman alabilir, size önerim yukarıdaki komutları ve komut gruplarını değiştirerek ve farklı dosyalarla denemeniz. 


Bu bölümde kullandığımız program: cutadapt
Bu bölümde kullandığımız UNIX komutları: ls, wc, awk, paste, awk, head, sort, uniq, | ve > .
Bu bölümde kullandığımız dosyalar: SRR1513652.fastq.gz ve SRR1513652.seq .