การจัดกลุ่มข้อมูลสาย RNA โดยไม่ต้องทำ Alignment และ Folding

หมายเหตุ: ผู้เขียนไม่ใช่ผู้เชี่ยวชาญทางด้าน Bioinformatics และไม่มีความรู้พื้นฐานด้านชีววิทยา

พื้นฐานของวิธีในที่นี้มาจากผลงานที่ชื่อว่า NoFold: RNA structure clustering without folding or alignment (Kim and Middleton, 2014) ซึ่งได้อธิบายวิธีการจัดกลุ่มข้อมูล RNA โดยไม่ต้องทำ alignment และ folding ในบนความนี้ผู้เขียนได้เลือกเอาเฉพาะส่วนที่เป็นพื้นฐานของวิธีนี้ โดยไม่ได้ให้ความสำคัญกับรายละเอียดยิบย่อยมากนัก

เรารู้จักสาย DNA ถ้ามองในเชิงคอมพิวเตอร์หน่อยก็คือสาย string ของข้อมูล A, C, G, T ซึ่งต่อกันยาว ๆ ไปเรื่อย ๆ ใช้บ่งบองคุณลักษณะต่าง ๆ ของสิ่งมีชีวิต ในธรรมชาติแล้วสาย DNA จะ transcribe ไปเป็น RNA ในขั้นตอนต่อไป ซึ่งก็คือสาย string ของข้อมูล A, U, G, C ผู้เขียนไม่มีความรู้มากนักที่จะพูดในเรื่องนี้เราจะข้ามมันไป ฮ่ะ ๆ

ปัญหาการจัดกลุ่ม RNA เป็นปัญหาในเชิง Bioinformatics ซึ่ง ประโยชน์ของการจัดกลุ่ม RNA เหล่านี้ก็คือ หากเรานิยามกลุ่มของ RNA คือ RNA ที่มีคุณลักษณะหรือคุณสมบัติในทางต่าง ๆ ใกล้เคียงกัน ดังนั้นหากรู้คุณสมบัติของ RNA เส้นหนึ่งในกลุ่มแล้ว เราก็พอจะเห็นภาพของ RNA เส้นอื่น ๆ ในกลุ่มเดียวกันด้วย ซึ่งจะช่วยให้เราสามารถทำงานกับ RNA ได้ง่ายขึ้น ประหยัด และรวดเร็วมากขึ้น เช่น RNA เส้นที่ก่อโรค A ก็ควรจะอยู่ในกลุ่มเดียวกับ RNA ที่ก่อโรค A ด้วย (ซึ่งอาจจะไม่จริง) แต่อย่างไรก็ดีหากเราศึกษาเรื่องโรค A แล้ว เราก็จะสามารถจำกัดขอบเขตในการศึกษา RNA ได้แคบลงอย่างมาก

ปัญหาหลักของการจัดกลุ่ม RNA ก็คือ สาย RNA ไม่ได้อยู่บน euclidean space (เพราะว่าเป็น sequence ของอักขระ) ทำให้หลาย ๆ วิธีในการทำ clustering เช่น K-Means ไม่สามารถใช้ได้ ยังไม่ต้องกล่าวถึงว่าจะหาค่า K มาได้อย่างไร ถึงตรงนี้ทำให้การจัดกลุ่ม RNA เหมือนกับว่าใช้ได้เพียงกับ clustering ที่อาศัยเพียงระยะห่างระหว่างคู่จุดก็พอ เช่น hierarchical clustering

มีวิธีการมากมายที่ได้นำเสนอในการจัดกลุ่มของเส้น RNA เข้าด้วยกัน ซึ่งเกือบทุกวิธีจะอาศัยเทคนิกไม่ว่าการทำ alignment หรือ folding หรือทั้งคู่ (ซึ่งวิธีนี้ไม่ใช้) เข้ามาร่วมด้วยทั้งสิ้น ในการวัดความเหมือนต่างระหว่างคู่สาย RNA แล้วก็อาศัยการจัดกลุ่มโดย UPGMA ซึ่งก็คือการทำ hierarchical clustering โดยใช้ linkage แบบ average นั่นเองโดยเราจะอธิบายถึงความหมายของสิ่งดังกล่าวซักครู่ต่อไปนี้

การทำ Alignment

การทำ alignment ของสาย RNA ก็คือการหาความเหมือนของสาย RNA เหล่านั้นในระดับ base-pair หรือ A, U, G, C นั่นแหละ ยกตัวอย่างการทำ alignment ของ RNA 2 เส้นอย่างง่ายที่สุดอาจจะเป็นการหา Hamming distance ระหว่าง RNA 2 เส้นนั้น ก็คือจำนวนสัญลักษณ์​ A, U, G, C ที่ไม่ตรงกัน

แต่วิธีการหา Hamming distance นี้ไม่ได้ทำให้เราไปได้ไกลเท่าไหร่ก็เพราะว่า “ความเหมือน” ในเชิงนี้มักจะน้อยมาก ๆ และไม่ได้บอกอะไรมากนัก เราสามารถทำให้ดีกว่านี้ด้วยการหา Minimum edit distance ซึ่งสามารถหาได้ด้วยการ dynamic programming แต่ในที่นี้นำหนักของแต่ละการแก้ไขอาจจะไม่เท่ากัน เพื่อให้สะท้อนความเป็นจริงตามธรรมชาติมากขึ้น

อย่างไรก็ตามหากเราใช้การทำ alignment มาจัดกลุ่มสาย RNA ต่าง ๆ หลายงานวิจัยก็พบว่า หาสาย RNA ทั้งคู่ไม่ได้มีความเหมือนกันในระดับสูง ๆ (70%) อยู่แล้ว ก็ยากที่จัดกลุ่มด้วยวิธีการนี้ เพราะการทำ alignment จะไม่ให้ข้อมูลที่มีความหมายเท่าไหร่สำหรับสาย RNA ที่แตกต่างกันมาก ๆ

การทำ Folding

นอกการดูการจัดเรียงอักขระ A, U, G, C แล้วเรายังสามารถดูลักษณะการ “พับ” ของสาย RNA ได้ด้วย ซึ่งการศึกษาจำนวนมากให้การสนับสนุนว่าเป็นมาตรวัดความเหมือนต่างที่ดีกว่า alignment สำหรับ RNA ที่มีความต่างของ alignment มาก ๆ หน่อย

การทำ Folding ของสาย RNA คือ การคาดเดาผลหลังจากการ “fold” ของสาย RNA นี้แล้ว เนื่องจากการไป “ส่อง” ดูจริง ๆ ว่าสาย RNA fold แล้วมีลักษณะยังไงเป็นไปได้ยากและมีค่าใช้จ่ายสูง จึงมีการพัฒนาวิธีในการคะเนผลขึ้นมาแทน

ลักษณะโดยทั่วไปของวิธีการเหล่านี้คือการวิเคราะห์โครงสร้างที่เป็นไปได้ โดยเลือกโครงสร้างที่มีความ “เสถียร” สูงสุดตามหลัก thermodynamic นั่นก็คือมีผลรวม free energy น้อยที่สุด (ผู้เขียนเองก็ไม่รู้ว่า free energy คืออะไร -.-) แต่รักษาให้มันน้อยที่สุดคือโครงสร้างที่เราคาดหวังว่าจะเห็นได้มากที่สุดในธรรมชาติ ปัญหานี้จึงสามรถแก้ไขได้ด้วย dynamic programming

ลักษณะของโครงสร้างหลังจากทำ Folding แล้วจะเป็นดังภาพต่อไปน้ี

(เครดิตภาพ Zuker and Sankoff, 1984)

จะเห็นว่ามีตัวอักขระภาษาอังกฤษ B, I, H, M และลักษณะของขั้นบันได ซึ่งขั้นบันไดพวกนี้ก็คือแต่ละคู่ A, U, G, C ที่จับกัน และบริเวณที่เป็นตัวอักษรภาษาอังกฤษคือบริเวณที่จับกันไม่ได้ โดย B คือ bulge loop, I คือ interior loop, H คือ hair pin loop และ M คือ multiloop

อ่านเพิ่มเติมได้จาก Zuker, Michael, and David Sankoff. “RNA secondary structures and their prediction.” Bulletin of mathematical biology 46.4 (1984): 591–621.

เมื่อเราได้โครงสร้างแบบนี้มาแล้ว (หลายครั้งเรียกว่า secondary structure) แทนที่จะเปรียบเทียบความเหมือนของสาย RNA ก็เปรียบเทียความเหมือนของโครงสร้างนี้แทน ซึ่งก็มีหลายวิธี วิธีนึงที่เป็นไปได้คือการเปลี่ยนโครงสร้างนี้เป็น “กราฟ” แล้ว วัดความเหมือนของกราฟแทน (ซึ่งใช้ในวิธีที่ชื่อว่า GraphClust (Heyne and Steffen, 2012)) ความเหมือนของกราฟนี้สามารถหาได้จากวิธีชื่อว่า Graph kernel ซึ่งก็เป็นวิธีที่น่าสนใจ

สามารถอ่านเพิ่มเติม graph kernel ได้ที่ Costa, Fabrizio, and Kurt De Grave. “Fast neighborhood subgraph pairwise distance kernel.” Proceedings of the 26th International Conference on Machine Learning. Omnipress, 2010.

ทั้งวิธี alignment และ folding สามารถถูกใช้งานร่วมกันได้เพื่อให้ได้คำตอบที่ดียิ่งขึ้นไปอีกได้ โดยวิธีการนึงที่ใช้หลักการนี้ชื่อ LocARNA อ่านเพิ่มเติมได้จาก Will, Sebastian, et al. “Inferring noncoding RNA families and classes by means of genome-scale structure-based clustering.” PLoS Comput Biol 3.4 (2007): e65.

แต่สำหรับบทความนี้ เราจะกล่าวถึงวิธีที่ไม่ได้ใช้ทั้ง alignment และ folding เลย แต่เราจะอาศัยความช่วยเหลือจากฐานข้อมูล Rfam แทน… ด้วยวิธีนี้เราจะสามารถเปลี่ยนสาย RNA แต่ละสายให้เป็น “จุด”​ 1 จุด บน euclidean space ได้ และหลังจากนั้นเราจะสามารถใช้ขั้นตอนวิธีที่เราชื่นชอบในการ cluster จุดข้อมูลเหล่านี้ โดยไม่ต้องสนใจความเป็น RNA ของมันมาก่อนเลย !!

ฐานข้อมูล Rfam

คือฐานข้อมูลของ RNA ซึ่งถูกจัดหมวดหมู่เป็น family ไว้เรียบร้อยแล้ว (มากกว่า 2,000 families) ซึ่งแต่ละ family นั้นจะมีสาย RNA เพียงจำนวนหนึ่งที่ถูกจำแนกโดยมนุษย์ ซึ่งมีความน่าเชื่อถือ เรียกว่า “seed” แต่อย่างไรก็ดีการจำแนกด้วยมือทั้งหมดนั้นย่อมทำได้ยาก และลงทุนสูง วิธีการที่ Rfam ใช้คือการสร้าง Covariance Model (CM) ขึ้นมาซึ่งเปรียบได้กับ Hidden Markov Model (HMM) สำหรับฐานข้อมูล Pfam สำหรับโปรตีน แต่อย่างไรก็ดี HMM กับ CM ไม่ใช่สิ่งเดียวกัน แต่เป็น model เชิงสถิติเหมือนกัน ซึ่งสามารถใช้ทำนาย “ความน่าจะเป็น” ที่สาย RNA หนึ่ง ๆ จะเป็น RNA สมาชิกใน family ของ CM หนึ่ง ๆ

หลังจากได้ CM ของแต่ละ family มาจาก seed แล้ว Rfam ใช้ CM เหล่านี้ทำนาย RNA เส้นอื่น ๆ ที่ไม่ทราบ family อีก เพื่อหาเส้นที่มีความมั่นใจว่าเป็น family นี้มารวมเข้าไปในฐานข้อมูลอีก เรียกว่า “full” ทำให้ฐานข้อมูลใหญ่ขึ้นโดยไม่ต้องใช้แรงคนมากขึ้น

สามารถอ่านเพิ่มเติมได้จาก Griffiths-Jones, Sam, et al. “Rfam: an RNA family database.” Nucleic acids research 31.1 (2003): 439–441.

Covariance Model (CM)

เราจะไม่กล่าวถึงวิธีการสร้างเชิงทฤษฎี CM ในที่นี้ ซึ่งมีความลึกซึ้งเกินความจำเป็นในหัวข้อนี้ แต่ก็มีความจำเป็นที่จะทำให้เห็นภาพรวมของ CM ซักหน่อย

สิ่งที่ CM ทำได้ก็คล้ายกับสิ่งที่ Hidden markov model (HMM) ทำได้ สำหรับ HMM หากได้รับข้อมูล เรียกว่า​ “observation” เข้ามาแล้วสามารถตอบได้ว่าลำดับของ state ใด และลำดับของ action ใดใน state machine มีความน่าจะเป็นสูงสุดที่จะเป็นสาเหตุให้เราเห็น “observation” นั้นได้

One of the applications of HMM, is to predict the sequence of state changes, based on the sequence of observations.”

ทั้ง CM และ HMM ถูกสร้างด้วยวิธีการคล้าย ๆ กัน คือต้องออกแบบ “states” และ “transitions” ของ states เล่านั้นให้ได้ แล้วนำ “observation” ประเภทเดียวกันจำนวนมากมา “สอน”​ มัน …​แล้วมันก็จะสามารถจำรูปแบบของ observation ที่มันถูกสอนมาได้ และใช้คะเน observation อื่น ๆ ในอนาคตได้

ซึ่งในที่นี้ CM ถูกนำมาใช้ในการ ดูว่า “สาย RNA นี้เข้ากันได้ดีกับสมาชิกอื่น ๆ ในกลุ่มที่ CM นี้เป็นตัวแทนได้ดีขนาดไหน” โดยสามารถบอกคะแนนออกมาเป็น “bitscore” ได้ ซึ่งหากมีค่ามากก็แสดงว่า CM มีความมั่นใจต่อสาย RNA นี้มาก ซึ่งก็แปลว่าสาย RNA นี้มีความคล้ายคลึงกับสาย RNA อื่น ๆ ที่ใช้สร้าง CM นี้มาก

สามารถอ่านเพิ่มเติมได้จาก Eddy, Sean R., and Richard Durbin. “RNA sequence analysis using covariance models.” Nucleic acids research 22.11 (1994): 2079–2088.

วิธีการสร้างและใช้งาน CM

มีเครื่องมือที่ได้รับความนิยมในการทำงานนี้อยู่แล้วชื่อว่า Infernal (http://infernal.janelia.org) ซึ่งถูกพัฒนาโดยผู้พัฒนาเดียวกับ Rfam ดังนั้นจึงไม่แปลกใจว่า CM ที่ถูกกล่าวถึงใน Rfam ก็คือ CM ที่ถูกพัฒนามาจาก Infernal ด้วย เครื่องมือนี้จึงเป็นทางเลือกที่เหมาะสม

Infernal มาพร้อมกับโปรแกรมย่อย ๆ เช่น (และอื่น ๆ แต่ผู้เขียนไม่มีความรู้มากไปกว่านี้)

  1. cmbuild คือ โปรแกรมสำหรับการสร้าง CM จากชุดข้อมูล RNA ที่เรามี ที่เราคิดว่าเป็นกลุ่มเดียวกัน และต้องการสร้าง CM เพื่อใช้ในการสะกัดลักษณะจำเพาะของ RNA ในกลุ่มนี้ออกมา ซึ่ง CM นี้สามารถใช้ในการจำแนกสาย RNA ใหม่ ๆ เข้ามาในกลุ่มเดียวกันนี้ได้ด้วย
  2. cmscore (v 1.0.2) คือ โปรแกรมสำหรับการเปรียบเทียบสาย RNA กับ CM ที่ต้องการ เพื่อหาความคล้ายคลึงของสาย RNA นั้นกับกลุ่ม RNA ที่ CM นี้เป็นตัวแทน

วิธีการใช้งานเต็ม ๆ สามารถหาได้จาก Userguide ของโปรแกรม Infernal เอง และในคู่มือนี้เองก็ได้อธิบายวิธีการสร้าง CM ของ Infernal ไว้ด้วย ซึ่งจะช่วยให้เราเข้าใจ CM มากขึ้น

การแปลงสาย RNA ให้อยู่บน euclidean space

วิธีการแปลงนี้อาศัยหลักบอกตำแหน่งของสาย RNA บน euclidean space ด้วยการบอก “ระยะห่าง” เทียบกับ “จุดสำคัญต่าง ๆ” วิธีนี้มีหลักการทำงานคล้าย ๆ ดังต่อไปนี้ สมมติว่ามีสถานที่ที่ทุกคนรู้จักอยู่ 3 ที่ แทนที่เราจะบอกตำแหน่ง GPS ของเราเอง เราสามารถบอกระยะห่างจากจุด 3 จุดนี้แทน ซึ่งก็อาจจะเพียงพอในการระบุตำแหน่งอย่างไม่กำกวมของเราแล้ว

หากแต่ในที่นี้เราใช้จุดบอกตำแหน่งราว 2 พันจุด ซึ่งจุดเหล่านั้นก็คือ CM ที่ถูกสร้างมาสำหรับแต่ละ family ของ Rfam หากเราได้สาย RNA มาสายหนึ่ง เราสามารถหาค่านี้ได้ด้วยการหา cmscore ของทุก ๆ CM กับสาย RNA สายนี้ซึ่งเราก็จะได้เป็นค่า bitscore ราว 2 พันค่าเช่นกัน bitscore ของ CM ของ family ไหนมีค่ามากก็แสดงว่าเราใกล้ family นั้นมาก …

แต่แท้จริงแล้วเราไม่จำเป็นต้องสนใจความหมายของ bitscore เลยก็ได้ เราสนใจเพียงแค่ หากเรานำเลข bitscore 2 พันค่านี้มามัดรวมกัน เราก็เหมือนจะได้พิกัดตำแหน่งของจุดหนึ่งจุดบน euclidean space 2 พันมิติแล้ว !!

ถึงตรงนี้เราก็สามารถใช้ขั้นตอนวิธีทาง Machine Learning ต่าง ๆ ในการจัดกลุ่มข้อมูลเหล่านี้โดยไม่ต้องสนใจเลยว่าก่อนเคยเป็นสาย RNA มาก่อน

แต่เดี๋ยวก่อน! ค่า bitscore อาจจะแกว่งในช่วงที่กว้างได้อย่างมาก (และอาจจะไม่เท่ากัน) เราควรจะป้องกันหรือลดผลกระทบจาก “bias” ที่อาจจะเกิดขึ้นได้ ด้วยการ normalize ข้อมูลเหล่านี้ก่อน ซึ่งวิธีที่งานวิจัย NoFold ใช้ก็คือการทำ Z-normalize ซึ่งก็ทำได้ดังสมการต่อไปนี้ ซึ่งจะปรับให้ช่วงของค่าลดลงอย่างมาก

หาก 2 พันมิติยังเยอะเกินไป! เราสามารถใช้ principal component analysis (PCA) เพื่อลดจำนวนมิติจากนี้ลงไปได้อีก ซึ่งงานวิจัย NoFold เสนอว่าสามารถทำให้เหลือเพียง 100 มิติได้โดยไม่มีปัญหาใด ๆ

ถึงตอนนี้เราสามารถแปลงสาย RNA ให้อยู่ในรูปของ euclidean space 100 มิติได้แล้ว! และหากเราจะจัดกลุ่มต่อด้วยการใช้ hierarchical clustering แบบไหน หรือ k-means แบบใด หรือแม้แต่วิธีการเหนือชั้นอื่น ๆ ก็ย่อมทำได้ แต่นี่ก็ยังไม่ได้กล่าวถึงวิธีการเลือกจำนวนกลุ่มของข้อมูลที่ขั้นตอนวิธีเหล่านี้จำต้องกำหนด

ถึงอย่างไรก็ดีนี่ก็เป็นจุดเริ่มต้นที่ดี (และแปลกใหม่) ที่ทำให้เราสามารถมอง RNA ซึ่งเป็น string แท้ ๆ เป็นจุดบน euclidean space ได้ …

สุดท้ายนี้ …​ สำหรับคนที่สนใจลงมือทำจริง ๆ ผู้เขียนแนะนำให้ลองเริ่มจากการใช้ NoFold ซึ่งสามารถหาโค้ดได้จาก repository ดังต่อไปนี้ https://github.com/sarahmid/nofold