Yoğun geçen birkaç haftadan sonra tekrar merhaba. Son yazımızda elimizdeki veriden adaptörleri çıkarmış ve her bir okumadan kaçar tane olduğunu hesaplamıştık. Bu yazıda birkaç ufak temizliğin ardından elimizdeki okumaları nihayet kullanılabilir bir FASTA formatına dönüştüreceğiz, böylece bir sonraki yazımızda bu verileri C. elegans genomuna hizalayıp üzerinden biyolojik çıkarımlar yapabileceğimiz bir veriye sahip olabileceğiz.
Geçen yazıda bir şeyin dikkatinizi çektiğini tahmin ediyorum: ACTG istisnasız her okumanın başında yer alıyor. Yeni Nesil Sekanslama'nın en tatsız taraflarından biri, teknolojinin göreceli olarak yeniliğinden ötürü zaman içerisinde standartların birçok kez değişmesi. Bu nedenle birkaç yıl geriye gittiğinizde bile farklı uygulamalarla karşılaşabiliyorsunuz. Üzerinde çalıştığımız bu dosya, aslında ana projeye ait olan ve çok daha büyük bir dosya grubunun sadece bir parçası: Burada bahsettiğim uygulamaları kolaylıkla ve kısa sürede gerçekleştirebilmek adına ufak bir dosyayla çalışmayı tercih ettim. Bu çalışmada birden fazla örnek aynı anda sekanslandığı için her bir örnek bir barkod ile işaretlenmiş, yani her bir okumanın başına bu barkod deneysel olarak eklenerek sekanslama öyle gerçekleştirilmiş. Bu dosyada kullanılan barkod da karşımıza ACTG olarak çıkıyor. Daha güncel veya farklı konfigürasyonlarda genelde barkod FASTQ dosyasında @ ile başlayan tanımlama satırlarda yer alabiliyor, birden fazla örnek içeren çalışmalarda bu tarz farklılıklara dikkat etmek önem arzediyor.
Gelelim nasıl devam edeceğimize. İlk olarak yapacağımız şey her bir okumadaki ilk dört bazdan kurtulmak. Bunun için bir awk fonksiyonu olan substr'yi kullanacağız. Bu fonksiyon üç parametre ile kullanılabilir: üzerinde işlem yapmak istediğimiz veri, bu verinin kaçıncı karakterinden itibaren almak istediğimiz, ve bunu kaçıncı karaktere kadar yapmak istediğimiz. Yani substr($0, 1, 4) bize her bir satırın ilk dört karakterini verir. Eğer durmak istediğimiz karakter sayısını, yani son parametreyi belirtmezsek, bu durumda o satırın sonuna kadar olan kısmı elde ederiz. Bizim yapmak istediğimiz şey tam da bu: ilk dört karakteri bırakıp kalanlar üzerinden çalışmak, bu durumda aşağıdaki gibi bir komuta ihtiyacımız var:
awk '{print substr($0,5)}' SRR1513652.seq | more
Komutun son kısmına eklediğim | more ifadesi olmasaydı dosyada yeniden düzenlenen her bir satır ekrana yazılacaktı, bunu istemediğim için her seferinde sadece bir ekran dolusu veriyle devam etmek için more komutunu kullandım. Boşluk tuşuna [space bar] bastığınızda bir sonraki sayfayı görüntüleyebilir, q harfine basarak da görüntülemeyi durdurabilirsiniz.
Dosyayı bu haliyle yeniden kaydedebiliriz, ancak işleri pratikleştirmek adına komutları ardarda ekleyerek devam etmek istiyorum. Yapacağımız şey bu verileri yeniden saymak, ardından da FASTA formatına dönüştürmek olacak. Saymak için bir önceki yazıdaki komutları yukarıdaki işlemin sonuna ekleyeceğiz (ilk sort komutunun yerine dikkat, ve tabii ki more komutunu da çıkarıp sona ekledim):
awk '{print substr($0,5)}' SRR1513652.seq | sort | uniq -c | sort -k1nr | more
FASTQ verilerini FASTA formatına dönüştürürken her bir okuma sayısını da okumanın yanına koymayı tercih ediyorum, böylece daha sonra okuma sayısına ulaşmak istediğimde fazladan bir adımla uğraşmak zorunda kalmıyorum. FASTA formatındaki okumaların şu şekilde olmasını istiyorum:
>TCACCGGGTGTACATCAGCTAA:44384
TCACCGGGTGTACATCAGCTAA
>GACTATCCTTTGATTGAGTTTTTTGT:9775
GACTATCCTTTGATTGAGTTTTTTGT
Her bir okumaya farklı bir isim vermek yerine, doğrudan okumanın içeriğini ve kaç kez okunduğunu yazmak çok daha pratik ve bizi birçok farklı aşamadan kurtarıyor. Okumayı ve okuma sayısını birbirinden ayırmak için : karakterini kullanıyorum: Daha sonra hizalanan her bir okumaya ilişkin bilgi sahibi olmak istediğimde ilgili kısmı pratik olarak alabileceğim.
Bunu yapmak için, aşağıdaki gibi elde ettiğim satırları iki sütuna bölüp, önce > karakterini, ardından her bir satırın içinci sütununda yer alan okumayı, ardından : karakterini, ve her bir satırın ilk sütununda yer alan okuma sayısını yan yana getireceğim. Bunların yanısıra, istediğim veriyi oluşturmak amacıyla bir alt satıra geçmek için \n karakterine ve sekansın kendisini tekrar koymaya ihtiyacım var. awk kullanırken yan yana getirmek istediğim karakter grupları veya değişkenleri aralarına herhangi bir karakter koymadan yazabiliriz. Yani mesela şöyle bir ifadeyi yazdırmaya ihtiyacımız olacak: ">"$2":"$1"\n"$2 . Artık bir önceki komutun sonuna da bunu ekleyerek herşeyi toparlayabiliriz:
awk '{print substr($0,5)}' SRR1513652.seq | sort | uniq -c | sort -k1nr | awk '{print ">"$2":"$1"\n"$2}' > SRR1513652.fa
Yukarıdaki satırda italik olarak gösterdiğim kısım yeni eklediğimiz kısım. Burada awk komutunu tekrar kullandık, böylece bir önceki komuttan elde edilen çıktıyı tekrar başka bir awk komut dizisine göndererek istediğimiz çıktıyı elde etmiş olduk. İşimiz bittiği için de elde ettiğimiz çıktıyı bir dosyaya yazdırdık. Haydi dosyanın içeriğine bir göz atalım ve ilk 20 satırı görüntüleyelim:
head -n 20 SRR1513652.fa
Elde ettiğimiz çıktı:
>TCACCGGGTGTACATCAGCTAA:44384
TCACCGGGTGTACATCAGCTAA
>TCACCGGGAGAAAAACTGGAGT:34004
TCACCGGGAGAAAAACTGGAGT
>TCACCGGGTGAACACTTGCAGT:33222
TCACCGGGTGAACACTTGCAGT
>TCACCGGGTGAAAATTCGCATG:32559
TCACCGGGTGAAAATTCGCATG
>TCACCGGGTGGAAACTAGCAGT:15257
TCACCGGGTGGAAACTAGCAGT
>GAGACTATCCTTTGATTGAGTTTTTTGT:10642
GAGACTATCCTTTGATTGAGTTTTTTGT
>GACTATCCTTTGATTGAGTTTTTTGT:9775
GACTATCCTTTGATTGAGTTTTTTGT
>CACCCGTACATATGTTTCCGTGCT:9282
CACCCGTACATATGTTTCCGTGCT
>AGACTATCCTTTGATTGAGTTTTTTGT:8789
AGACTATCCTTTGATTGAGTTTTTTGT
>ACTATCCTTTGATTGAGTTTTTTGT:8175
ACTATCCTTTGATTGAGTTTTTTGT
İşte bu! Artık verilerimiz C. elegans genomuna hizalanmaya hazır. Bir sonraki yazımızda bunu nasıl yapacağımızı anlatacağım, ancak öncesinde elimizdeki sekansların uzunluğuna bir göz atmak istiyorum. Bu kısmı keşfetmeyi size bırakacağım ve ne yaptığımızı benzer bir örnekle ilerideki bir yazıda anlatacağım:
awk '{print substr($0,5)}' SRR1513652.seq | awk '{count[length($0)]++}END{for(seq in count)print seq"\t"count[seq]}' | sort -k2nr
Elde ettiğim çıktı, elimizdeki verilerde en fazla 22 baz uzunluğunda okumaları gördüğümüzü ve bunu 21 ve 23 bazlık okumaların takip ettiğini söylüyor. Dikkatli okuyucuların hemen farkettiği üzere kullandığımız veri kısa RNA'ların [small RNAs] tespiti amacıyla üretilmişti, bu nedenle böyle bir tabloyla karşılaşmayı zaten bekliyorduk.
Bu bölümde kullandığımız UNIX komutları: awk, more, sort, uniq, ve head.
Bu bölümde kullandığımız dosyalar: SRR1513652.seq ve SRR1513652.fa