Python Programming

นับจำนวน nucleotide ด้วย Python

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

สวัสดีครับ ครั้งนี้ในส่วนของ Python programming จะมีการเปลี่ยนแปลงไปพอสมควร เพราะผมคิดว่าเราเรียนรู้พื้นฐาน Python กันมาพอสมควรแล้ว ความจริงก็คือผมคิดว่าบทเรียนที่ผมเขียนค่อนข้างหยาบสำหรับให้น้องๆ ในหน่วยเรียนรู้ Python เริ่มจะหมด stock แล้ว และพอหวนกลับไปดูในส่วนที่เหลือ ผมคิดว่าอาจจะไม่ค่อยมีประโยชน์เท่าไรนัก เลยตัดสินใจเขียนขึ้นมาใหม่ สำหรับท่านที่สนใจบทความที่เหลือใน stock แวะเข้าไปดูได้ที่ http://bdm.si.mahidol.ac.th/document/python/python-tutorial ครับ ถ้ามีคำแนะนำ ติชมประการใด จะส่งเมล์หา editor คนเก่งของเรา หรือส่งเมล์ถึงผมโดยตรงก็ได้ครับ

โจทย์ - นับจำนวน bp ใน FASTA file

โจทย์ครั้งนี้ไม่ยากครับ แต่น่าจะช่วยให้เข้าใจการเขียนโปรแกรมได้ดีขึ้น และเราจะได้เรียนคำสั่ง Python เพิ่มเติมจากที่เคยกล่าวถึงหลายคำสั่งทีเดียว

มาเริ่มที่โจทย์นะครับ มีคนนำไฟล์ FASTA format ที่มี sequence ของ nucleotide อยู่มาให้ครับ แล้วบอกว่าอยากทราบว่า nucleotide ในแต่ละไฟล์นั้นมีความยาวกี่ bp ไม่ยากใช่ไหมครับ ทำเองก็ได้ แต่ถ้ามีสัก 10 ไฟล์ คงไม่ค่อยอยากทำเองสักเท่าไหร่ เรามาเขียนโปรแกรมช่วยกันดีกว่า

ศึกษาโจทย์

ก่อนที่เราจะเขียนโปรแกรม เราต้องเข้าใจโจทย์ให้ชัดเจนก่อน ไฟล์ FASTA format จะประกอบด้วย 2 ส่วนคือ ส่วนที่เป็น header และส่วนที่เป็นเนื้อหา ส่วน header จะขึ้นต้นด้วยสัญลักษณ์ > แล้วตามด้วยคำบรรยาย หลังจากนั้นจะเป็นส่วนเนื้อหาทั้งหมด ที่สำคัญเนื้อหาสามารถมีมากกว่าหนึ่งบรรทัดได้

ตัวอย่าง FASTA file

>gi|296176325|emb|HC742679.1| Sequence 184 from Patent WO2010043977

TCCTATGTGCTGACTCAGCCACCCTCGGTGTCAGTGGCCCCAGGACAGACGGCCAGGATTACCTGTGGGG

GAAACAACATTGGAAGTAAAAGTGTGCACTGGTACCAGCAGAAGCCAGGCCAGGCCCCTGTGCTGGTCGT

CTATGATGATAGCGACCGGCCCTCAGGGATCCCTGAGCGACTCTCTGGCTCCAACTCTGGGAACACGGCC

ACCCTGACCATCAGCAGGGTCGAAGCCGGGGATGAGGCCGACTATTACTGTCAGGTGTGGGATAGTAGTA

GTGGTCCTTTTGTGGTTTTCGGCGGAGGGACCAAGCTGACCGTCCTAG

ถ้าเราต้องการจะนับจำนวน nucleotide เราก็จะต้องนับ nucleotide ในส่วนที่เป็นเนื้อหาทั้งหมด

เริ่มที่การวางแผน

วิธีหนึ่งที่ใช้ได้ในการวางแผนคือถามตอบตัวเอง โจทย์คือเราต้องการนับจำนวน nucleotide ในส่วนเนื้อหา

คำถาม จะรู้ได้อย่างไรว่าเป็นเนื้อหา

คำตอบ บรรทัดนั้นไม่ได้ขึ้นต้นด้วยสัญลักษณ์ >

คำถาม แล้วบรรทัดว่างที่คั่นระหว่าง header กับเนื้อหาล่ะ

คำตอบ ถ้าบรรทัดนั้นไม่มีตัวอักษรอะไร เราก็ไม่สนใจ

จากตรงนี้เราก็เริ่มได้แผนแล้วว่า การทำงานของโปรแกรมควรจะเป็นว่า อ่านข้อมูลจากไฟล์ทีละบรรทัด แล้วดูว่าบรรทัดนั้นขึ้นต้นด้วยสัญลักษณ์ > หรือไม่ ถ้าใช่ก็ข้ามไปอ่านบรรทัดถัดไป ในทำนองเดียวกัน ถ้าบรรทัดนั้นไม่มีตัวอักษรอะไร ก็ข้ามไปอ่านบรรทัดถัดไปเช่นกัน เขียนออกมาคร่าวๆ ได้เป็น

1__for line in fasta_file:

2______if line.startswith('>'):

3__________continue

4______if line.empty():

5__________continue

6______count_nucleotide

โปรแกรมข้างต้นยังไม่ใช่โปรแกรม Python ที่ใช้งานได้จริง แต่เป็นการแปลงแผนของเราให้ออกมาอยู่ในรูปของโปรแกรมแบบคร่าวๆ เรียกว่า pseudocode

คำถาม แล้วเราจะอ่านไฟล์ได้อย่างไร

คำตอบ ใช้คำสั่ง open สมมติว่าไฟล์ที่เราต้องการอ่านอยู่ใน drive c ในโฟลเดอร์ data ชื่อว่า test.fasta เราจะสามารถอ่านได้โดยคำสั่ง

fasta_file = open('c:/data/test.fasta')

ให้สังเกตว่าเราใช้เครื่องหมาย slash (/) ไม่ใช่ backslash (\) ตามที่เราใช้กันปกติใน Windows

คำถาม เราเรียกใช้คำสั่ง empty ในบรรทัดที่ 4 ใน python ไม่มีคำสั่งนี้ เราจะทำอย่างไร

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

คำถาม ค่าการขึ้นบรรทัดใหม่???

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

คำถาม แล้วเราจะตรวจสอบยังไง

คำตอบ ปกติเราจะไม่ตรวจสอบค่าการขึ้นบรรทัดใหม่ เมื่อเราต้องการรู้ว่าบรรทัดนั้นว่างหรือเปล่า โดยทั่วไปแล้ว เราจะให้ Python ตัดค่าการขึ้นบรรทัดใหม่ทิ้งไป ด้วยคำสั่ง strip() แล้วตรวจสอบดูว่าไม่มีข้อมูลเหลือ (empty string) โปรแกรมในบรรทัดที่ 4 จึงเปลี่ยนเป็น

if line.strip() == '':

Note: คำสั่ง strip() เป็นคำสั่งที่ใช้ในการกำจัด white space ที่อยู่ด้านหน้าและหลังของบรรทัด คุณสามารถหาข้อมูลเพิ่มเติมได้จาก help file ของ python ในครั้งนี้ เราต้องการใช้คำสั่ง strip() ในการกำจัดค่าการขึ้นบรรทัดใหม่ซึ่งถือว่าเป็น white space ตัวหนึ่งเหมือนกัน ค่า white space คือค่าต่างๆ ที่ไม่ได้พิมพ์ออกมาให้เห็น เช่นตัว space เองที่เป็นตัวว่างๆ ค่า Tab ที่ทำให้เกิดช่องว่างหลายช่อง และค่าการขึ้นบรรทัดใหม่ ถ้ามีค่าเหล่านี้ที่ด้านหน้าหรือหลังของบรรทัด คำสั่ง strip() จะกำจัดให้เรา คำสั่ง strip() จะเป็นคำสั่งที่เราใช้งานบ่อยมากเมื่อทำงานกับไฟล์ตัวอักษร

คำถาม แล้วเราจะนับจำนวน nucleotide ยังไง

คำตอบ Python มีคำสั่ง len ที่สามารถใช้นับตัวอักษรได้อยู่แล้ว คำสั่งบรรทัดที่ 6 จึงเป็น

print len(line)

โปรแกรมเวอร์ชันแรก และการแก้ไขข้อผิดพลาด

มาถึงตรงนี้ เราก็จะได้โปรแกรมแรกดังนี้

listing 1.

1__fasta_file = open('c:/data/test.fasta')

2__for line in fasta_file:

3______if line.startswith('>'):

4__________continue

5______if line.strip() == '':

6__________continue

7______print len(line)

ลอง save file นี้ในชื่อ count_fasta.py (ไม่ต้องมีบรรทัดที่เขียนว่า listing 1. และไม่ต้องมีเลขบรรทัดนะครับ) แล้วลองเอาโปรแกรมนี้ไปใช้งานจริง และตรวจคำตอบก็จะพบข้อผิดพลาดและความไม่สะดวกหลายอย่าง ประการแรกคือจำนวน nucleotide ที่พิมพ์ออกมาของแต่ละบรรทัดเกินมา 1 ตัวอักษร และโปรแกรมไม่ได้นับจำนวน nucleotide ทั้งหมดมาให้ แต่กลับพิมพ์ความยาวของ nucleotide แต่ละบรรทัดออกมา แถมยังผิดเสียด้วย

คำถาม ทำไมนับเกินมา 1

คำตอบ เพราะว่าคำสั่ง len นับค่าการขึ้นบรรทัดใหม่ ถึงแม้ว่าเราจะมองไม่เห็นค่านี้ แต่ python มองเห็น และคำสั่ง len ก็ทำงานถูกต้อง เพียงแต่เราไม่ได้นึกถึง จึงได้คำตอบที่ผิดออกมา เราสามารถแก้ไขด้วยการใช้คำสั่ง

print len(line.strip())

แต่คำสั่งด้านบน ยังไม่ตอบโจทย์เรื่องการนับจำนวน nucleotide ทั้งหมด เราต้องแก้ไข ด้วยการยังไม่พิมพ์ค่าที่นับได้ของแต่ละบรรทัด แต่ให้สะสมยอดไว้ก่อน เพื่อพิมพ์ผลในภายหลัง เราจึงแก้ไขโปรแกรมด้วยการเพิ่มตัวแปรตัวหนึ่งสำหรับนับจำนวน nucleotide

listing 2.

1__fasta_file = open('c:/data/test.fasta')

2__count = 0

3__for line in fasta_file:

4______if line.startswith('>'):

5__________continue

6______if line.strip() == '':

7__________continue

8______count = count + len(line.strip())

9__print count

ให้สังเกตว่า บรรทัดที่ 9 ไม่ได้ย่อหน้าเข้าไป เพราะเราต้องการให้ออกจากลูป for ก่อน จึงค่อยพิมพ์ค่าใน count การทำงานของโปรแกรมคือให้ค่า count = 0 เสียก่อนในบรรทัดที่ 2 จากนั้นทุกครั้งที่การทำงานมาถึงบรรทัดที่ 8 Python ก็จะนำเอาค่า count ที่มีอยู่ไปบวกกับความยาวของ nucleotide ของบรรทัดนั้น แล้วใส่กลับไปที่ตัวแปร count อีกครั้งหนึ่ง (ให้อ่าน code จากทางด้านขวาของเครื่องหมาย = ก่อน เพื่อดูว่าค่าใหม่ที่จะใส่ให้ตัวแปรที่อยู่ทางซ้ายคืออะไร) เมื่อสิ้นสุดการอ่านไฟล์ ก็จะออกมาจากลูป for Python ก็จะพิมพ์ค่าที่อยู่ใน count ออกมา

ปรับแต่งโปรแกรม ให้ดูเป็นมืออาชีพ

อันที่จริงโปรแกรมใน listing 2. ก็ทำงานได้ถูกต้องแล้ว แต่สามารถปรับปรุงให้โปรแกรมดูดีขึ้นอ่านง่ายขึ้นได้

จุดแรกคือเราจะเห็นว่าเราเรียกใช้ line.strip() 2 แห่ง เราสามารถนำค่า line.strip() ใส่ไว้ในตัวแปรตัวหนึ่งแล้วนำไปใช้ได้ จะทำให้โปรแกรมกระชับขึ้น จุดที่ 2 คือคำสั่ง if ที่บรรทัดที่ 4 และ 6 นั้นถ้าเป็นจริง จะทำงานเหมือนกันคือ continue เราควรจะรวบเงื่อนไขบรรทัดที่ 4 และ 6 ได้ สุดท้ายจุดที่ 3 คือเมื่อตัวแปร count ถูกใช้ในลักษณะนี้ เราสามารถเขียนแบบย่อได้ ซึ่งเมื่อเราคุ้นเคยกับ python แล้วจะอ่านได้เข้าใจง่ายกว่า เพราะจะทราบได้ทันทีว่าเราต้องการ update ตัวแปร count เราจึงได้โปรแกรมสุดท้ายออกมาเป็น

listing 3.

1__fasta_file = open('c:/data/ test.fasta')

2__count = 0

3__for line in fasta_file:

4______data = line.strip()

5______if line.startswith('>') or data == '':

6__________continue

7______count += len(data)

8__print count

เพิ่มความสะดวก ด้วยการอ่านชื่อไฟล์จาก command line

การที่เราจะต้องแก้ชื่อไฟล์ทุกครั้งที่จะนับจำนวน nucleotide ช่างไม่สะดวกเอาเสียเลย ทางที่ดี เราน่าจะสามารถเรียกใช้งานได้ในลักษณะนี้ คือ เปิด command box ขึ้นมา (รันโปรแกรม cmd.exe) แล้วเรียกใช้โปรแกรมได้ โดยใส่ชื่อไฟล์ fasta ต่อท้าย ดังตัวอย่าง

python count_fasta.py c:\data\test.fasta

listing 4.

1__import sys

2__fasta_file = open(sys.argv[1])

3__count = 0

4__for line in fasta_file:

5______data = line.strip()

6______if line.startswith('>') or data == '':

7__________continue

8______count += len(data)

9__print count

ให้สังเกตบรรทัดที่ 1 และ 2 ครับ ผมขออธิบายคร่าวๆ ในครั้งนี้ก่อนว่า บรรทัดที่ 1 เป็นการบอกให้ Python นำ module ที่ชื่อว่า sys เข้ามา และบรรทัดที่สอง เราเปิดไฟล์ที่มีชื่ออยู่ใน sys.argv[1] แทนที่เราจะกำหนดไว้ตายตัว ท่านที่สนใจรายละเอียดสามารถค้นข้อมูลเพิ่มเติมได้ครับ เราจะกล่าวถึงรายละเอียดของ module sys และการอ่านค่าจาก command line ในภายหลัง

แต่ครั้งนี้ คุณได้โปรแกรมนับจำนวน nucleotide ที่ใช้งานได้จริงเก็บไว้ใน stock เรียบร้อยแล้ว