Saturday, October 30, 2010

Python Programing

Reverse complement สาย DNA ด้วย Python (ภาคจบ)


ประพัฒน์ สุริยผล


ภาพยนตร์ส่วนใหญ่มักจะจบกันที่ไตรภาค แต่ reverse complement ของเราผ่านไปแล้ว 3 ภาคยังจบไม่ลงสักที เราก็จะมาจบกันที่ครั้งนี้ครับ หวังว่าภาคจบนี้ เราคงได้เรียนรู้อะไรกันอีกพอสมควร


Reverse sequence revisited


ผมมีจังหวะแวะเข้าไปดู blog ของ THAI Bioinformatics e-magazine แล้วพบว่ามี comment มาจากคุณ Tanawut แนะนำให้ผมใช้ trick หนึ่งของ Python ซึ่งผมก็ลืมนึกถึงไปเลย จึงขออนุญาตเอามาลงในครั้งนี้ด้วยครับ


code สุดท้ายที่เราเขียนไว้สำหรับ reverse_sequence เป็นดังด้านล่าง

def reverse_sequence(dna):

____return ''.join(reversed(list(dna)))

แต่ code ที่คุณ Tanawut แนะนำมาจะเป็น


def reverse_sequence(dna):

____return dna[::-1]

สั้นลงไปอีก แต่อาจจะสับสนเล็กน้อย สำหรับคนที่ไม่คุ้นเคยกับ Python ตัวผมเองคิดว่า code นี้สั้นและเป็นสไตล์ Python แท้ๆ เลยครับ โดยใช้คุณสมบัติของ string ใน Python ที่สามารถมองเป็น list ได้ ต้องขอบคุณคุณ Tanawut มากครับ


Slicing


เรามาทำความเข้าใจเรื่อง slicing กันเล็กน้อย จากตัวอย่างครับ


>>> a = ['a', 'b', 'c', 'd']

>>> a

['a', 'b', 'c', 'd']

>>> a[0:2]

['a', 'b']

>>> a[:2]

['a', 'b']

>>> a[2:4]

['c', 'd']

>>>a[2:]

['c', 'd']


การ slicing คือการดึงข้อมูลของ list ในส่วนที่เราต้องการ จากตัวอย่าง เราสร้างตัวแปร a ให้เป็น list ของตัวอักษร เมื่อเราต้องการใช้ slicing ให้เรากำหนดตำแหน่งที่เราต้องการภายในวงเล็บแบบสี่เหลี่ยม [] โดยกำหนดตำแหน่งเริ่มต้น และตำแหน่งสุดท้าย


อย่างเช่น a[0:2] แปลว่า เราต้องการ slicing ที่ตำแหน่ง 0 (ตำแหน่งแรกสุด) ไปจนถึงก่อนตำแหน่งที่ 2 (ก็คือตัว 'c' เพราะเรานับจาก 0) สิ่งเราได้มาก็จะเป็น list ของ ['a', 'b']


ถ้าหากเราละไว้ ไม่ใส่ตัวเลขที่ตำแหน่งเริ่มต้น ค่าปกติก็จะเป็น 0 เพราะฉะนั้น คำสั่ง a[:2] จึงให้คำตอบเหมือนกัน a[0:2] อ่านว่า เราต้องการ slicing list ตั้งแต่ตำแหน่งเริ่มต้นถึงก่อนตำแหน่งที่ 2


ในทำนองเดียวกัน เราสามารถละตำแหน่งสุดท้ายได้ ซึ่งจะเป็นการบอกว่า เราต้องการ slicing ไปจนถึงตัวสุดท้ายของ list เพราะฉะนั้น ค่าของ a[2:4] และ a[2:] จึงเหมือนกัน


ผู้อ่านบางท่านอาจจะนึกสงสัยว่า แล้วถ้าเราละตัวเลขทั้งสองตำแหน่งล่ะ จะเป็นอย่างไร คำตอบก็คือเราจะได้ copy ของ list ตัวนั้น ซึ่งจะไม่ใช่ list ตัวเดียวกันแต่มีค่าเหมือนกัน งงใช่ไหมครับ วิธีนี้เป็น trick ที่มีประโยชน์และช่วยให้เข้าใจเรื่องของ list ได้ดีทีเดียว ลองดูตัวอย่างกันต่อนะครับ


>>> a = ['a', 'b', 'c', 'd']

>>> b = a

>>> c = a[:]

>>> a

['a', 'b', 'c', 'd']

>>> b

['a', 'b', 'c', 'd']

>>> c

['a', 'b', 'c', 'd']

จากตัวอย่าง จะเห็นว่าตัวแปร a, b, c มีค่าเหมือนกันหมดถูกไหมครับ แต่เราสร้างตัวแปร b กับ c คนละแบบกัน


ตัวแปร b เรากำหนดให้เท่ากับตัวแปร a เลย แต่ตัวแปร c เราให้เป็น copy ของตัวแปร a


ให้ลองนึกว่า list ก็คือกล่องใบนึงนะครับ เรามีกล่องที่แปะชื่อเอาไว้ว่าชื่อ a ในกล่อง a มีลูกบอลอยู่ 4 ลูก เขียนไว้ว่า a, b, c และ d


พอเรากำหนดให้ b = a ก็เหมือนกับว่า เราเอาชื่อ b ไปแปะไว้ที่กล่องเดิมนั่นแหละครับ กล่องนั้นจึงมีชื่อทั้ง a และ b โดยที่ลูกบอลในกล่องก็เป็นชุดเดียวกัน


แต่ถ้าเรากำหนดให้ c = a[:] อันนี้ เหมือนกับว่า เราหากล่องมาอีกกล่องครับ ข้างในมีลูกบอล 4 ลูก เขียนว่า a, b, c และ d เช่นกัน และกล่องนี้แปะชื่อว่า c ถ้าดูที่ลูกบอลในกล่องจะพบว่าทั้งกล่อง a, b และ c เหมือนกันหมด


ความแตกต่างจะอยู่ตรงนี้ครับ ถ้าเราเปลี่ยนแปลงของในกล่อง a ของในกล่อง b ก็จะเปลี่ยนไปด้วย เพราะมันคือกล่องใบเดียวกัน นึกออกไหมครับ ในขณะที่ของในกล่อง c ไม่ถูกกระทบไปด้วย ดูตัวอย่างนะครับ


>>> a[0] = 'e'

>>> a

['e', 'b', 'c', 'd']

>>> b

['e', 'b', 'c', 'd']

>>> c

['a', 'b', 'c', 'd']

เรื่องนี้สำคัญนะครับ โปรแกรมทำงานผิดพลาดมีสาเหตุสำคัญจากการไม่เข้าใจเรื่องของ list ในลักษณะกล่องนี่ล่ะครับ ทำให้เราเปลี่ยนค่าไปโดยที่เราไม่ตั้งใจ


ทุกครั้งที่เราเปลี่ยนค่าในกล่อง ไม่ว่าจะเป็นตัวแปร a หรือ b จะส่งผลกระกบทั้งคู่ เพราะทั้งสองชื่อหมายถึงกล่องเดียวกัน เราสามารถแยกไม่ให้มีผลกระทบได้ ก็คือใช้วิธี copy กล่องแบบที่เราทำกับตัวแปร c ครับ เราจะมาลองกับตัวแปร b กันอีกทีครับ


>>> b = a[:]

>>> a[-1] = 'f'

>>> a

['e', 'b', 'c', 'f']

>>> b

['e', 'b', 'c', 'd']

>>> c

['a', 'b', 'c', 'd']

เราเริ่มต้นด้วยคำสั่ง b = a[:] เท่ากับว่าเราเอากล่องมาใหม่อีกหนึ่งใบ ข้างในมีลูกบอลเหมือนที่ a มีอยู่ในขณะนี้ และแกะชื่อ b ที่เดิมแปะอยู่กับกล่องเดียวกันกับ a มาแปะเอาไว้ที่กล่องใหม่นี้ เท่ากับว่าตอนนี้เรามี 3 กล่องที่ไม่เกี่ยวข้องกันแล้ว จากนั้น เราลองเปลี่ยนค่า item ตัวสุดท้ายของตัวแปร a โดยเรากำหนดค่า index เป็น -1 หมายถึง นับจากข้างหลังครับ -1 หมายถึงตัวสุดท้าย -2 หมายถึงตัวก่อนสุดท้าย ไล่ไปเรื่อย นึกออกไหมครับ หลายครั้งการใช้ index ในลักษณะนี้ก็จะทำให้เราสะดวกมากขึ้น เพราะเราไม่ต้องรู้ว่า list นี้ยาวเท่าไหร่ ตัวสุดท้ายอยู่ตำแหน่งอะไร


เรามาลองดูค่าสุดท้ายของตัวแปรแต่ละตัวกันครับ ตัวแปร a ก็จะเป็นตัวแปรที่มีการเปลี่ยนแปลงทั้งตัวแรกและตัวสุดท้าย ตัวแปร b มีแค่ตัวแรกที่เป็น 'e' เพราะว่าเรา copy จากตัวแปร a หลังจากที่เราเปลี่ยนตัวแรกไปแล้ว ส่วนตัวแปร c ไม่มีการเปลี่ยนแปลงใดๆ


เรามาถึง trick สุดท้ายของการทำ slicing ครับ คือนอกจากเรากำหนดว่าเราจะ slicing จากไหนถึงไหนได้แล้ว เรายังกำหนด step ในการอ่านค่า slicing ได้ด้วย งง อีกแล้วครับ ดูตัวอย่างดีกว่า


>>> c[0:2:1]

['a', 'b']

>>> c[2:0:-1]

['c', 'b']

>>> c[::-1]

['d', 'c', 'b', 'a']

>>> c[::2]

['a', 'c']

>>> c[::-2]

['d', 'b']

ผมเลือกใช้ตัวแปร c นะครับ เพราะค่าเรียงจาก a ไป d ดูตามง่ายกว่า


คำสั่งแรก c[0:2:1] จะให้ผลไม่ต่างจาก c[0:2] ที่เราเคยเรียกใช้กัน เพราะว่าค่าปกติของ step ก็คือ 1 อยู่แล้ว Python ก็จะ slicing จากตำแหน่งที่ 0 แล้วก็เพิ่มทีละ 1 ไปจนถึงก่อนตำแหน่งสุดท้ายที่กำหนด


มาดูที่คำสั่งถัดไป c[2:0:-1] ครั้งนี้ เรากำหนด step เป็น -1 (ตัวเลขสุดท้ายหลังโคลอน) โดยเริ่มต้น slicing ที่ตำแหน่งที่ 2 คือตัว 'c' ด้วย step ติดลบ ตำแหน่งถัดไปจึงเป็น 1 ก็คือตัว 'b' ตำแหน่งถัดไปคือ 0 แต่ว่า slicing จะดึงค่ามาจนถึงก่อนตำแหน่งสุดท้าย จึงไม่รวมตำแหน่งที่ 0 มา คำตอบที่ได้จึงเป็น ['c', 'b']


ถ้าหากเราละค่าตำแหน่ง เช่นเรียกใช้ c[::1] ก็จะไม่ต่างอะไรจาก c[:] เพราะค่า step 1 เป็นค่าปรกติอยู่แล้ว แต่ว่าถ้าเราใช้ c[::-1] เราก็จะได้ copy ของ list ที่ reverse เพราะจะวิ่งจากตัวสุดท้ายย้อนมาตัวหน้าสุด คุณสมบัติอันนี้ น่าสนใจมากครับ


เราสามารถกำหนด step ที่ไม่ใช้ 1 และ -1 ได้ ดังตัวอย่างครับ ผมให้เป็นหน้าที่ของท่านผู้อ่านที่จะลองไล่ดูนะครับว่า ถ้า step เป็น 2 และ -2 แล้วได้คำตอบดังตัวอย่างได้อย่างไร


เกริ่นมายาวมากเลย เพื่อมาปิดท้ายที่การใช้ string ในลักษณะของ list ผู้อ่านจะเห็นว่าเราสามารถใช้ slicing กับ string เช่นกัน และผลลัพธ์ที่ได้ ก็จะถูกแปลงกลับมาให้เป็น string เหมือนเดิม สะดวกมากครับ ดังตัวอย่าง


>>> b = 'abcd'

>>> b[0:2]

'ab'

>>> b[:2]

'ab'

>>> b[2:4]

'cd'

>>> b[2:]

'cd'

เพราะฉะนั้น ถ้าเราต้องการ reverse ตัวแปร string เราทำอย่างไรครับ ก็เพียงแค่ใช้คำสั่งแบบที่คุณ Tanawut แนะนำครับ


reverse_string = b[::-1]

เรามาต่อกันที่ function complement ที่ครั้งที่แล้วเราปรับปรุงด้วยการใช้ dictionary ช่วย ซึ่งโปรแกรมครั้งที่แล้วมีจุดบอดอยู่สองสามแห่ง หนึ่งคือ key ของ dictionary ที่เราใช้เป็นเบสตัวเล็กทั้งหมด ถ้าหากเราได้ข้อมูลสาย DNA เป็นเบสตัวใหญ่ เช่น ATGC แทนที่จะเป็น atgc โปรแกรมของเราก็จะทำงานผิดพลาด อันที่สองคืออาจจะมีเบสแปลกๆ เข้ามาได้เช่น N ซึ่งโปรแกรมของเราไม่รองรับ สามคือถ้าหากเราต้องการคงสภาพตัวเล็กตัวใหญ่ไว้ เช่นถ้าเบสเข้ามาเป็น A เราก็ต้องได้ค่า complement เป็น T แต่ถ้าเข้ามาเป็น a ให้คืนค่าเป็น t เราก็ต้องมาปรับโปรแกรมกันพอสมควร จากของเดิม


def complement(dna):

____complementary_dict = {'a':'t', 't':'a', 'g':'c', 'c':'g'}

____result = ''

____for base in dna:

________result = result + complementary_dict[base]

____return result


เราต้องเพิ่ม complementary_dict เป็นดังนี้ครับ

complementary_dict = {'a':'t', 't':'a', 'g':'c',\

______________________'c':'g', 'A':'T', 'T':'A',\

______________________'G':'C', 'C':'G', 'n':'n',\

______________________'N':'N'}


ถ้าเป็นเบสตัวอื่นนอกเหนือจากนี้ โปรแกรมของเราก็จะหยุดทำงาน พร้อมกับรายงานว่าไม่มี key ตัวนั้นอยู่ ผมขอปล่อยให้เป็นแบบนั้นก็แล้วกันนะครับ เพราะปกติเราก็พบเบสเพียงแค่เท่านี้ ถ้าหากมีอะไรแปลกปลอมเข้ามา โปรแกรมก็จะโวยวายมาให้เราทราบ


ตอนนี้เรามาลองดูว่า นอกจาก dictionary สำหรับใช้ในงานนี้แล้ว Python มีอะไรให้เราใช้อีกบ้าง เนื่องจากเราต้องการเปลี่ยนค่าใน string เมื่อไปดูคู่มือ Python พบว่า Python เตรียมคำสั่ง translate เราไว้ให้เราใช้แล้ว


>>> a = 'aabbccdd'

>>> import string

>>> table = string.maketrans('ab', 'ef')

>>> a.translate(table)

'eeffccdd'

การใช้งาน เริ่มต้นจากการเตรียมตารางสำหรับแปลงค่า วิธีง่ายที่สุดคือเรียกใช้คำสั่ง maketrans ใน module string โดยเรา import module string เข้ามา แล้วเรียกใช้ดังตัวอย่าง พารามิเตอร์ที่เราส่งไปให้คำสั่ง maketrans จะมีสองค่า ค่าแรกจะเป็นค่าเริ่มต้น ค่าที่สองจะเป็นค่าที่ต้องการเปลี่ยน ความยาวของทั้งสองค่าต้องเท่ากัน เพราะ maketrans จะจับคู่ไล่ไป ดังเช่นในตัวอย่าง เราต้องการจะเปลี่ยนจาก a เป็น e และ b เป็น f จึงส่งค่า 'ab' และ 'ef' ไป


จากนั้น เราสามารถเรียกใช้คำสั่ง translate ได้เลย ส่วนตัวอักษรตัวอื่นที่ไม่ได้อยู่ในส่วนที่เราต้องการให้เปลี่ยนแปลงก็จะคงอยู่เหมือนเดิม


ดังนั้น ถ้าเราต้องการนำคำสั่ง translate มาใช้กับ function complement ของเราบ้าง เราก็จะได้ดังนี้


import string

complement_table = string.maketrans('atgcATGC', 'tacgTACG')

def complement(dna):

____return dna.translate(complement_table)


ผมแยกบรรทัดที่ import กับบรรทัดที่สร้างตัวแปร complement_table ออกมาข้างนอก เพราะถ้าเราใส่ไว้ใน function จะถูกเรียกใช้ทุกครั้งที่เราเรียกใช้ function complement โดยอันที่จริงแล้ว เรา import ครั้งเดียวก็พอ และตัวแปร complement_table ก็ตายตัว สร้างครั้งเดียวแล้วก็ไม่เปลี่ยนแปลง สามารถนำไปใช้ได้เรื่อยๆ ผมจึงแยกออกมา


function complement อันนี้ ต่างจากอันที่ใช้ dictionary คือมันจะไม่โวยวายถ้าได้ค่าเบสที่มันไม่รู้จัก เพราะคำสั่ง translate จะคืนค่าเดิมกลับมา เช่นในตารางเราไม่ได้ใส่เบส N เพราะว่าค่าที่คืนกลับมาเป็นค่าเดิม


สุดท้ายแล้ว โปรแกรม reverse complement ของเราจะหน้าตาเป็นแบบนี้ครับ


import string

COMPLEMENT_TABLE = string.maketrans('atgcATGC', 'tacgTACG')

def complement(dna):

____return dna.translate(COMPLEMENT_TABLE)


def reverse_sequence(dna):

____return dna[::-1]

dna = raw_input('Please enter DNA sequence -> ')

reversed_dna = reverse_sequence(dna)

reverse_complement_dna = complement(reversed_dna)

print reverse_complement_dna

ถ้าไปเปรียบเทียบกับครั้งแรกที่เราเขียน จะเห็นว่าโปรแกรมสั้นลงไปมาก และอ่านเข้าใจได้ง่ายขึ้นด้วย แต่ถ้าอยากให้สั้นกว่านี้อีก เราก็ตัดนิยาม function ทิ้ง แล้วเรียกใช้แบบ chain of function ได้ตามสะดวกครับ


import string

COMPLEMENT_TABLE = string.maketrans('atgcATGC', 'tacgTACG')

dna = raw_input('Please enter DNA sequence -> ')

print dna.translate(COMPLEMENT_TABLE)[::-1]

เหลือ 4 บรรทัดครับ และยังคงอ่านได้ไม่ยากนัก Python ก็ดีอย่างนี้นี่เอง แล้วพบกันใหม่ฉบับหน้าครับ

No comments:

Post a Comment