ทำ reverse complement สาย DNA ด้วย Python
ประพัฒน์ สุริยผล
สวัสดีอีกครั้งครับ ครั้งที่แล้ว เราเริ่มจะใช้ Python ทำงานทางด้าน Bioinformatics กันแล้ว ถึงแม้ว่าโจทย์จะไม่ได้ยากหรือซับซ้อนอะไร แต่ก็ทำให้เราได้เรียนรู้การเขียนโปรแกรมไปใช้งานจริงได้พอสมควร ครั้งนี้ เรายังคงวุ่นวายอยู่กับสาย DNA ครับ สิ่งที่เราพบบ่อย คือเราต้องการดูเส้น reverse complement ของสาย DNA ที่เรามีอยู่ เช่น
5' ATTCGTAAGCATTAC 3'
3' TAAGCATTCGTAATG 5'
เมื่อเรามีสายด้านบนอยู่ เราต้องการแปลงให้ได้สายด้านล่าง ถ้าเส้น DNA ค่อนข้างยาว จะทำด้วยมือก็คงเสียเวลาและมีโอกาสผิดพลาดได้สูง วิธีแก้ปัญหาที่ส่วนใหญ่มักจะทำกันคือ เข้าไปใช้บริการเว็บหลายๆ แห่งที่มีคำสั่งทำ reverse complement ให้ หรือไม่ก็เรียกใช้โปรแกรมบางตัวเช่น BioEdit ขึ้นมาช่วย
ต่อไปนี้ เราจะไม่ต้องเสียเวลาเข้าเว็บหรือเปิดโปรแกรม BioEdit อีกแล้ว เพราะเราจะเขียนโปรแกรมเล็กและสั้นด้วย Python ให้เราสามารถทำ reverse complement ได้อย่างรวดเร็ว ทุกครั้งที่ต้องการ
นิยามโจทย์ให้ชัดเจน
เราต้องตัดสินใจก่อนว่า เราต้องการ input คือค่าของสาย DNA จากทางไหน จะให้โปรแกรมเปิดไฟล์ที่มีข้อมูลของ DNA แล้วพิมพ์ผล reverse complement ออกทางหน้าจอ หรือจะให้โปรแกรมอ่านข้อมูลจากทางคีย์บอร์ดเลย การแสดงผลจะให้พิมพ์ออกหน้าจอหรือให้เก็บลงไฟล์ จะเห็นว่าการเขียนโปรแกรม ถึงแม้ว่าจะเป็นโปรแกรมเล็กๆ ก็มีปัญหาเรื่องการนิยามโจทย์ให้ชัดเจน แต่หากเราเริ่มเชี่ยวชาญแล้ว เราอาจจะนิยามไว้ในใจก็ได้ ถ้าเป็นโปรแกรมขนาดเล็ก แต่ถ้าเป็นโปรแกรมขนาดใหญ่ เราจำเป็นต้องนิยามให้ชัดเจน มิฉะนั้นจะมีปัญหาในตอนเขียนโปรแกรมครับ
เนื่องจากผู้ใช้โปรแกรมนี้คือตัวเราเอง ผมก็เลยสรุปไปเลยว่า ให้อ่านจากคีย์บอร์ดน่าจะสะดวกที่สุดในตอนนี้ เพราะผมมักจะ copy สาย DNA จากเว็บหรือจากเมล์ที่คนส่งมา แล้วต้องการหาเส้น reverse complement และการแสดงผล ผมกำหนดให้เป็นการพิมพ์คำตอบออกบนหน้าจอ ถ้าหากต้องการเก็บผลลงไฟล์ เราสามารถใช้ copy & paste ไปใส่ในโปรแกรมพวก notepad หรือ Microsoft Word ได้ หรือจะใช้วิธีที่สะดวกกว่าคือการ pipe ลงไฟล์ (ผมจะกล่าวถึงเรื่อง pipe ในครั้งต่อๆ ไปครับ)
เริ่มเขียนโปรแกรม
เช่นเคยครับ ผมลองเขียน pseudocode ออกมาก่อนดังด้านล่าง
1 dna = read_input()
2 reversed_dna = reverse_sequence(dna)
3 reverse_complement_dna = complement(reversed_dna)
4 print reverse_complement_dna
ขั้นตอนการทำงานครั้งนี้ดูง่ายมากใช่ไหมครับ มีแค่ 4 ขึ้นตอน อ่านค่า reverse เส้น DNA แล้วแปลงค่า complementary base เสร็จแล้วก็พิมพ์ผลออกมา แต่สิ่งสำคัญที่เราจะเรียนครั้งนี้คือ function ครับ
ก่อนอื่น เรามาเริ่มต้นที่การรับ input จากผู้ใช้ก่อน เราจะใช้คำสั่ง raw_input ครับ เราเคยใช้กันไปครั้งหนึ่งแล้ว ตอนที่เราคำนวณพื้นรูปต่างๆ จำได้ไหมครับ ในครั้งนี้ เราเขียนโปรแกรมออกมาได้บรรทัดเดียวเป็น
dna = raw_input('Please enter DNA sequence -> ')
เขียนฟังก์ชัน (Function)
จากนั้น เรากำหนดค่าตัวแปร reversed_dna ด้วยการเรียกใช้คำสั่ง reverse_sequence ซึ่งเป็นคำสั่งที่ไม่มีอยู่ใน python เป็นหน้าที่ของเราที่จะต้องเขียนขึ้นมาเองครับ การเขียนคำสั่งใหม่หรือฟังก์ชันใหม่ใน python มีรูปแบบดังนี้ครับ
def ชื่อฟังก์ชัน(ตัวแปรต่างๆ ที่จะส่งมาให้ฟังก์ชัน):
code ของฟังก์ชัน
ในกรณีนี้ฟังก์ชันของเราต้องการสาย DNA เข้ามา แล้วจะต้องทำงานส่งคืนสาย DNA อีกสายหนี่งที่ reverse กับสายที่ส่งเข้ามา เราเขียนได้เป็น
1 def reverse_sequence(dna):
2 result = “”
3 for base in dna:
4 result = base+result
5 return result
เรามาลองไล่การทำงานของ code ในฟังก์ชันดูนะครับ
บรรทัดที่ 1 นิยามฟังก์ชัน โดยบอกว่าฟังก์ชันชื่อ reverse_sequence แล้วต้องการตัวแปรหนึ่งตัวชื่อ dna
บรรทัดที่ 2 กำหนดตัวแปรอีกหนึ่งตัวชื่อว่า result ที่ผมตั้งชื่อนี้ เพราะตัวแปรนี้จะเก็บผลลัพธ์ที่จะส่งคืน
บรรทัดที่ 3 – 4 เป็นการวนลูป โดยแต่ละครั้ง base จะเป็นตัวอักษรตัวหนึ่งใน dna (ที่เป็นข้อมูลตัวอักษรสตริง เพราะคำสั่ง raw_input จะคืนค่าที่เป็นตัวอักษรสตริงกลับคืนมา) เช่นถ้า dna มีค่าเป็น 'atgc' ในการวนลูปครั้งแรก base จะมีค่า a เมื่อทำงานมาถึงบรรทัดที่ 4 โปรแกรมจะนำเอาค่า base มาบวกกับค่าใน result ซึ่งในที่นี้ก็คือเอามาต่อกันนั่นเอง แต่เนื่องจาก result มีค่าเป็น empty string (คือไม่มีตัวอักษรอะไรอยู่เลย จากบรรทัดที่ 2) ค่าที่ได้จึงเป็น a และเก็บค่าไว้ใน result
เมื่อวนลูปครั้งที่ 2 ตัวแปร base จะมีค่าเป็นตัวอักษรถัดไปคือ t เมื่อทำงานที่บรรทัดที่ 4 ก็จะเป็นการนำค่า base ปัจจุบันคือ t ไปรวมกับค่าใน result ซึ่งตอนนี้มีค่า a อยู่ นำผลลัพธ์ไปใส่ใน result ซึ่งจะได้เป็นค่า ta
เมื่อวนลูปครั้งที่ 3 result ก็จะถูกปรับค่าเป็น gta แล้วครั้งที่ 4 ก็จะได้ค่าเป็น cgta เมื่อวนลูปครั้งที่ 4 เสร็จแล้ว ปรากฏว่าไม่มีตัวอักษรให้วนลูปได้อีก จึงออกจากลูปมาที่บรรทัดที่ 5 ซึ่งมีการใช้คำสั่ง return เพื่อคืนค่าใน result กลับไปให้กับผู้เรียกใช้ ค่าใน result ขณะที่หลุดจากลูป คือ cgta ซึ่งเป็นสาย reverse ของ atgc ที่เราใส่เข้ามา ผลลัพธ์ที่ได้ก็จะไปอยู่ในตัวแปรชื่อ reversed_dna
ฟังก์ชันถัดไปที่เราต้องการก็คือ complement เขียนได้เป็น
1 def complement(dna):
2 result = “”
3 for base in dna:
4 if base == 'a':
5 result = result + 't'
6 elif base == 't':
7 result = result + 'a'
8 elif base == 'g':
9 result = result + 'c'
10 elif base == 'c':
11 result = result + 'g'
12 else:
13 result = result + 'n'
14 return result
การทำงานของฟังก์ชัน complement จะคล้ายกับฟังก์ชัน reverse_sequence ด้วยการวนลูปทีละ base แต่โปรแกรมที่อยู่ในลูปจะต่างออกไป สิ่งที่น่าสังเกตคือชื่อตัวแปรที่ส่งมาที่ฟังก์ชันในบรรทัด
reverse_complement_dna = complement(reversed_dna)
คือ reversed_dna แต่ชื่อที่เราตั้งไว้ตอนรับค่าในฟังก์ชัน complement คือ dna จะเห็นว่าไม่จำเป็นต้องเป็นชื่อเดียวกัน เมื่อมีการเรียกใช้ฟังก์ชัน ตัวแปร dna ที่เราใช้งานในฟังก์ชันจะมีค่าเท่ากับ reversed_dna ที่ส่งมา
สมมติว่า ตัวแปร reversed_dna มีค่า cgta เมื่อวนลูปครั้งแรก ตัวแปร base จะมีค่าเป็น c เมื่อเปรียบเทียบเงื่อนไขแล้ว จะไปตรงกับบรรทัดที่ 10 เมื่อเงื่อนไขเป็นจริง จึงไปทำงานที่บรรทัดที่ 11 ทำให้ค่า result เป็น g เพราะค่า result ตอนแรกเป็น empty string
เมื่อวนลูปครั้งที่ 2 ตัวแปร base เป็น g จะไปตรงกับเงื่อนไขบรรทัดที่ 8 ทำให้ค่า result เป็น gc เพราะว่าค่า result ก่อนหน้านี้คือ g เมื่อมาทำงาน��ี่บรรทัดที่ 9 ก็จะเท่ากับว่า result = 'g' + 'c'
วนลูปครั้งที่ 3 ตัวแปร base เป็น t ตรงกับเงื่อนไขบรรทัดที่ 6 ทำงานที่บรรทัดที่ 7 ทำให้ค่า result เป็น gca เพราะค่าเดิมของ result จากลูปครั้งก่อนคือ gc เมื่อนำมาต่อด้วย a จึงได้เป็น gca
วนลูปครั้งที่ 4 เราก็จะได้ค่า result เป็น gcat
หลุดจากลูป ทำการคืนค่า result ในบรรทัดที่ 14
บรรทัดที่ 12 และ 13 เป็นการเขียนกันเอาไว้ ในกรณีที่ผู้ใช้อาจจะพิมพ์ค่าเบสอื่นๆ ที่ไม่ใช่ atgc ในกรณีนี้ผมตัดสินใจว่า ถ้าเกิดกรณีนี้ขึ้น ให้ใส่ค่า n สำหรับค่าที่เป็น complementary base เราอาจจะตัดสินใจหยุดโปรแกรมแล้วแสดงข้อความผิดพลาดเตือนผู้ใช้ก็ได้ แต่ครั้งนี้เพื่อให้ง่าย ผมจึงตัดสินใจแบบนี้ไปก่อน
การเขียนโปรแกรมจุดสำคัญอีกประการหนึ่ง โดยเฉพาะในส่วนของโปรแกรมที่เกี่ยวข้องกับการตัดสินใจ คือต้องคิดถึงกรณีต่าง ๆ ที่เป็นไปได้ให้รอบคอบ ถ้าหากเราไม่ได้คิดมากนัก เราก็อาจจะคิดว่าตัวเบสในสาย DNA มีเพียงแค่ atgc แต่เป็นไปได้ว่า ผู้ใช้อาจจะพิมพ์ผิด หรือตั้งใจพิมพ์ตัวอักษรอื่นที่ไม่ใช่ atgc ก็ได้ เราจึงต้องเขียนโปรแกรมเพื่อดักกรณีเหล่านี้ให้เรียบร้อย ทักษะในการคิดให้รอบคอบเช่นนี้มีความสำคัญและจำเป็นมากสำหรับการเขียนโปรแกรมให้ทำงานได้อย่างถูกต้อง
กลับมาที่ฟังก์ชัน เมื่อทำงานมาจนถึงบรรทัดสุดท้าย ค่าที่ประมาณได้จากฟังก์ชันก็จะถูกคืนค่ากลับมาที่ตัวแปร reverse_complement_dna ในบรรทัดที่ 3 และพิมพ์ผลลัพธ์ในบรรทัดที่ 4 ใน pseudocode ที่เราร่างไว้
3 reverse_complement_dna = complement(reversed_dna)
4 print reverse_complement_dna
สำหรับท่านผู้ที่เขียนโปรแกรมเป็นแล้วหรือพอเป็นบ้าง จะเห็นว่าโปรแกรมที่เขียนครั้งนี้ค่อนข้างยาว และไม่ “pythonic” เลย ประมาณว่าคนที่เขียน python เป็นแล้วมาเห็นโปรแกรมนี้ จะร้องยี้ว่าเขียนเหมือนกับภาษา C ช่างเขียนได้ยาวและไม่ได้ใช้ความสามารถของ python เลย ขอให้ใจเย็นๆ ครับ ครั้งหน้าเราจะมาปรับและเรียนรู้ลึกลงไปในภาษา python ด้วยกัน ดูว่า โปรแกรมของเราจะออกมารูปร่างหน้าตาเป็นอย่างไร
จุดที่น่าสนใจในการเขียนโปรแกรมครั้งนี้คือการที่เราค่อยๆ เติมเบสในลูปเข้าไปในตัวแปร result ขอให้สังเกตความแตกต่างของการปรับค่า result ของทั้งสองฟังก์ชันให้ดีครับ ถ้าหากเราเขียนสลับตำแหน่ง ผลลัพธ์ที่ได้ก็ผิดพลาดได้ทันที สำหรับครั้งนี้ เราได้โปรแกรมที่ทำงานได้ถูกต้องตามที่เราต้องการแล้ว พบกันครั้งหน้า ขอให้มีความสุขกับการเขียนโปรแกรมทุกๆ ท่านครับ
Listing program ทั้งหมด
1 def reverse_sequence(dna):
2 result = “”
3 for base in dna:
4 result = base+result
5 return result
6
7 def complement(dna):
8 result = “”
9 for base in dna:
10 if base == 'a':
11 result = result + 't'
12 elif base == 't':
13 result = result + 'a'
14 elif base == 'g':
15 result = result + 'c'
16 elif base == 'c':
17 result = result + 'g'
18 else:
19 result = result + 'n'
20 return result
21
22 dna = raw_input('Please enter DNA sequence -> ')
23 reversed_dna = reverse_sequence(dna)
24 reverse_complement_dna = complement(reversed_dna)
25 print reverse_complement_dna
