B i o i n f o r m a t i c s - เราเก็บข้อมูลทาง Genomics กันอย่างไร ?

เคยสงสัยกันมั้ยว่า เราเก็บข้อมูลทาง Genomic พวก DNA ของเราเป็นแบบ Digital ได้อย่างไร วันนี้เราจะพามาทำความรู้จักกับ Format ที่เราใช้งานอย่าง FASTA และ FASTQ ที่เป็นตัวหลัก ๆ ที่เราใช้งานกัน และ ถ้าเราอยากจะเขียนโปรแกรมมาอ่านไฟล์พวกนี้ เราจะเขียนออกมาอย่างไรดี

ปล. พวกรายละเอียดบางอย่าง ถ้าไม่รู้จักไม่เป็นไรเด้อ คิดซะว่าเหมือนดู Star Wars แหละ ตอนแรกว่าจะเขียนแค่เรื่องนี้ขำ ๆ แต่ไหน ๆ ก็ไหน ๆ แล้ว เขียนเป็น Series หลาย ๆ ตอนละกัน มีเรื่องเล่าเยอะ ก็ค่อยมาย้อนอ่านตอนนี้อีกทีก็ได้

ปล 2. ความรู้พวกนี้เป็นเรื่องจริง ที่เขียนโดยคนปลอม ๆ ที่ ใ จ เ ก เ ร มาก

FASTA Format

เริ่มจาก Format ที่น่าจะเข้าใจง่ายที่สุดก่อนดีกว่า คือ FASTA Format เป็น Text-Based Format ทั่ว ๆ ไป นั่นแปลว่า เราสามารถเปิดเข้ากับ Text Editor แล้วอ่านตรง ๆ ได้เลย ซึ่งจะเป็น Format ที่สามารถเก็บได้ทั้ง Nucleotide Sequence และ Peptide Sequence เอาตามใจสั่งเลยจ๊ะแม่ !

>Sample 1
ACTTGGGAGACCATTTAAGC

ใน FASTA ต่อ 1 Record จะประกอบด้วย 2 อย่าง บรรทัดแรก เราจะเริ่มต้นด้วยเครื่องหมายมากกว่าเสมอ เป็นการบอกว่า เราจะเริ่ม Record ใหม่แล้วนะ แล้วตามด้วย Header ที่ประกอบด้วย Infomation ของ Sequence นั้น ๆ เช่น Severe acute respiratory syndrome coronavirus 2 อะไรก็ว่ากันไป แล้วแต่เราจะตั้งชื่อเลย แล้วมันก็จะเป็น Pattern แบบนี้ลงไปเรื่อย ๆ แล้วแต่ที่เราต้องการจะเก็บ เช่น 10 Records หรือจะล้าน Records ก็ไม่ว่ากัน (ถ้าคอมพี่ไหว ไปก่อนเลย)

ดังนั้นในไฟล์ FASTA Format เราจะเก็บข้อมูลได้เพียงเท่านี้เลย ถามว่า แล้วถ้าเราจะใส่ Metadata อื่น ๆ เราจะทำได้มั้ยคำตอบคือ เราทำได้ แต่เราจะต้องใส่ไว้ใน Header แทน เช่น NC12345 | xxxx อะไรก็ว่ากันไป ก็ได้เหมือนกัน

input_file = open('data.fasta', 'r')

fasta_recs = {'headers' : [], 'sequences' : []}
n_lines = 0
line = ""

current_header = ""

while (line := input_file.readline()) != '' :
    n_lines += 1
    
    if n_lines % 2 == 0 :
    	current_header = line.rstrip()[1:]
    else : 
    	fasta_recs['headers'].append(current_header)
        fasta_recs['sequences'].append(line.rstrip())

input_file.close()

พอมาเขียนโปรแกรมเพื่ออ่านเข้ามาด้วยมือ มันก็ไม่ได้ยากขนาดนั้น เขียนเป็น Python Script ง่าย ๆ เรารู้อยู่แล้วว่า Record นึง บรรทัดแรกจะเป็น Header เราก็อ่านแล้วเก็บไว้ก่อน ซึ่ง ถ้า Header ของเรามันมี Format เราก็อาจจะเขียน Function อีกตัวมา Parse เป็น Format ที่เราต้องการได้ จากนั้น พอเราอ่านบรรทัดต่อไปเราก็รู้ละว่า เป็น Sequence ข้อมูลครบละ เราก็ยัดใส่ Record ได้เลย

ส่วน Data Structure ที่เก็บ เอาจริง ๆ ขึ้นกับความต้องการเราละว่าเราจะเก็บด้านไหน ใน Code ตัวอย่างเราเลือกที่จะเก็บแยก Header และ Sequence ออกจากกัน เพราะส่วนใหญ่เวลาเราทำงาน เราก็จะใช้เฉพาะ Sequence ซะเยอะ การที่เราดึง Sequence ทั้งยวงออกมาได้เลย น่าจะทำให้ทำงานได้ง่ายกว่า เร็วกว่า แต่ถ้าใครมีวิธีการใช้งานที่ไม่เหมือนก็อาจจะปรับเปลี่ยนไปตามการใช้งานของเราได้เด้อ

อันนี้เสริมให้หน่อย บางครั้ง เราอาจจะเจอตัวแปลก ๆ ที่ไม่ได้เป็นตัวย่อใน Nucleotide A,T,C และ G หรือแก๊งค์ Peptide เช่น N,S, W และ N จริง ๆ มีตัวอื่นอีก เอาจริง ๆ ใช้ ๆ ไป มันจำได้จริง ๆ นะ พูดเป็นเล่น

FASTQ Format

จากเดิม FASTA เราเก็บข้อมูลได้ค่อนข้างจำกัด ไม่เหมาะกับงานบางอย่าง เลยทำให้มันมี Format อีกตัว ก็คือ FASTQ ขึ้นมา เก็บทุกอย่างเหมือนกับ FASTA ทุกประการ แต่สิ่งที่เพิ่มคือ มันสามารถเก็บ Quality Score ได้ ซึ่งเลือก Encode ผ่าน ASCII Table ไป

ตัว FASTQ Format จะเป็น Format ที่ได้รับความนิยมที่เราใช้ในการทำงานกับพวกเครื่อง Sequencer ทั้งหลายเช่น Illumina และ Oxford Nanopore เป็นต้น เพราะมันสามารถเก็บ Quality Score ได้นี่แหละ

ถามว่า แล้ว Quality Score คืออะไร มันเป็นคะแนน ที่เครื่อง Sequencer จะให้คะแนนว่า มันมั่นใจกับ Base หรือข้อมูลในตำแหน่งนั้น ๆ มากแค่ไหน ซึ่งเราจะใช้สิ่งที่เรียกว่า PHRED Score ซึ่งสมการของมันจะเป็น Logarithm จากสมการด้านบนเลย เช่น ถ้าเราบอกว่า คะแนนของเราอยู่ที่ 10 เราก็จะแทน Q ด้วย 10 แล้วเราก็แก้สมการหา P เราก็จะได้ออกมาเป็น 1 ใน 10 หรือแปลง่าย ๆ คือ เครื่องคิดว่า ตำแหน่งนี้ มีโอกาสผิดคือ 1 ใน 10 คิดเป็น Percentage ก็จะเป็น เครื่องคิดว่า โอกาสที่จะไม่ใช่ Base ตัวนี้ 10% อีกมุมก็คือ เครื่องมั่นใจว่า 90% น่าจะถูก

ซึ่งคะแนน ก็จะไล่ตั้งแต่ 0-42 ถ้าเราจะเอาตัวเลขใส่เข้าไปในไฟล์ตรง ๆ เลย มันก็น่าจะเป็นเรื่องที่สิ้นเปลือง และ ทำงานยากมาก ๆ เช่น 1 กับ 10 เมื่อเราเขียนในไฟล์จริง ๆ มันจะเป็น 110 ถามว่า มันคือ 1 ทั้งหมด 2 ตำแหน่ง แล้วตามด้วย 0 หรือ เป็น 11 แล้วตามด้วย 0 แล้วแยกกันหมดเลย ดูเป็นปัญหาอยู่นะ

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

เพื่อให้ง่ายกว่าเดิม เลยเลือกที่จะกลับไปหาอะไรที่มัน Basic กว่านั้นหน่อยคือ ASCII Table เราก็ใช้อักษรแทนตัวเลขไปเลย ตามตาราง ASCII ไป ซึ่งเครื่อง Sequencer สมัยนี้ส่วนใหญ่ก็จะใช้สิ่งที่เรียกว่า ASCII_Base 33 กัน ง่าย ๆ เวลาเราจะหาว่า มันคือ Score เท่าไหร่ เราเอา ตัวอักษรนั้น แทนค่าใน ASCII Table ออกมาแล้วลบด้วย 33 เราก็จะได้ Quality Score ออกมาแล้ว ทำให้เราเก็บข้อมูลได้ง่ายขึ้นเยอะมาก

@1754ed0c-e78a-4f93-a956-249e17d2c608
CAGTATGCGCCCGTTCAATCGGTGGTGTTTTATGATCATACG
+
*%%$$$$$#$$&%&$$,.10%&///7;/3;,23%+%$#$$$%

ลักษณะของ FASTQ คล้าย ๆ กับ FASTA เลย แค่จากเดิมจุดเริ่มจะเป็นเครื่องหมายมากกว่า เราก็เปลี่ยนเป็น @ จากนั้น มันก็จะเป็น Header ซึ่งก็แล้วแต่เรา หรือเครื่อง Sequencer บางตัวมันก็จะมี Format ของมันอีกว่า มันจะประกอบด้วยเลขอะไรบ้าง เช่น Flow Cell Position อะไรก็ว่ากันไป บรรทัดต่อไปก็จะเป็นตัว Sequence แล้วตามด้วยเครื่องหมาย บวก ซึ่งตรงนี้เราสามารถใส่ข้อมูลไปได้เหมือนกัน บางครั้งเราก็จะใส่ให้เหมือนกับ Header เพื่อเวลาเรา Import หรือทำ Integrety Check เราจะได้มั่นใจได้ว่า Sequence และ Quality Score มันคือของ Record เดียวกัน ไม่ได้หลงกันมา และ บรรทัดสุดท้ายก็จะเป็น Quality Score โดยที่ ความยาว มันจะต้องเท่ากับ ความยาวของ Sequence นะ

input_file = open('APNNP_105402.fastq', 'r')

fasta_recs = {'headers' : [], 'sequences' : [], 'scores': []}
n_lines = 0
line = ""

current_header = ""
current_sequence = ""

while (line := input_file.readline()) != '' :
    n_lines += 1
    
    if n_lines % 4 == 0 :
    	current_header = line.rstrip()[1:]
    elif n_lines % 4 == 1 :
    	current_sequence = line.rstrip()
    elif n_lines % 4 == 3 :
        fastq_recs['headers'].append(current_header)
        fastq_recs['sequences'].append(line.rstrip())
    	fastq_rec['scores'] = line.rstrip()


input_file.close()

การเขียน Script มาเพื่ออ่าน เอาจริง ๆ ก็ไม่ได้ยากเลย เราแค่ปรับจาก Script ที่อ่าน FASTA นิดเดียวเอง เรารู้ว่า มันมี 4 บรรทัดต่อ Record บรรทัดแรกคือ Header เราก็อ่านเข้ามาตรง ๆ แล้วเรารู้ว่าตัวแรกสุดคือ @ ซึ่งเราไม่เอา เราก็ข้ามไป แล้วตามด้วยบรรทัดต่อไปเป็น Sequence เราก็อ่านเข้ามาตรง ๆ แล้วเราจะเจอเครื่องหมายบวก ในที่นี้เราไม่ขออ่าน เราเลยข้ามไปเลย และสุดท้าย เราเจอ Quality Score เราก็แค่อ่านเข้ามาแค่นั้นเลยหมดแล้ว

อยากลองเล่น หา Dataset มาจากไหน ?

พวก Dataset ที่เป็นพวก FASTA และ FASTQ เราหาได้เยอะมาก ๆ มีแหล่งข้อมูลให้เราลองโหลดมาเล่นเพียบเลย เช่นใน NCBI เราก็สามารถเข้าไปดาวน์โหลดพวก Sequence ของสิ่งมีชีวิตหลาย ๆ อย่างมาเล่นได้ เช่น COVID-19 ก็มีคน Sequence แล้วเอามาลงเหมือนกัน

ส่วนที่เป็นพวก FASTQ อันนี้ส่วนใหญ่จะได้มาจากเครื่อง Sequencer โดยตรงเลย มันก็จะมี Project ที่เอาข้อมูลพวกนี้มาปล่อยเหมือนกัน เช่น 1000Genome Project มี Genome ของมนุษย์ด้วยนะ

ปกติ เราเก็บข้อมูลเป็น FASTA หรือ FASTQ ?

โดยทั่ว ๆ ไปแล้ว ถ้าเราดูจาก Component ของแต่ละ Record เราจะเห็นว่า FASTQ น่าจะใหญ่กว่าแน่นอน เออ ก็ใช่แหละ อย่างที่บอกก็คือ มันได้มาจากเครื่อง Sequencer โดยตรงเลย ขึ้นกับความละเอียดในการ Sequence และ ขนาดของ Genome ที่เรา Sequence มาอีก ขนาดไฟล์ อาจจะใหญ่ไปได้ถึง 1 TB หรือมากกว่านั้นสำหรับ Genome Size แค่คนเองนะ ยิ่งถ้าใช้ในการวินิฉัยทางการแพทย์ เราอาจจะต้องการความละเอียดที่สูงมาก ๆ ทำให้ขนาดไฟล์ใหญ่เข้าไปได้อีก ขึ้นกับปัญหาที่เราต้องการเลย ไว้เราจะมาเล่าเรื่องความละเอียดให้อ่านกัน

อ่านแล้วจะรู้สึกว่า โอ้มันบ้าไปแล้ว เราจะเอาอะไรใหญ่ขนาดนั้น แต่อย่าลืมนะว่า ถ้าเรา Sequence Genome หรือ Gene ขนาดใหญ่จริง ๆ มันก็ใหญ่มากอยู่แล้ว เช่นคนเอง ก็มีถึง 3.3 Billion Base Pair ยังไม่นับว่า ถ้าเราทำแบบ Pair-end เราก็จะได้ออกมา 2 ด้าน 2 ไฟล์ เท่ากับว่า ถ้าไฟล์นึงขนาดสัก 1 TB ทำแบบ Pair-end เราจะต้องใช้พื้นที่อย่างน้อย 2 TB ในการเก็บเลยนะ

จะเห็นว่า FASTQ โคตรใหญ่ ใหญ่แบบบ้านชุบแป้งทอดเลย ทำให้ ถ้าเกิดว่า งานที่เราทำ Analysis ที่เราต้องการทำมันไม่มีความจำเป็นต้องใช้ เราก็อาจจะ Process มันให้เป็น FASTA แล้วเก็บ ส่วน FASTQ ต้นฉบับ เราก็อาจจะเก็บไว้ในพวก External HDD หรือใน Server ของ Lab (ถ้ามี) ก็ได้ พกไปมาทั้งยวงก็บ้าอยู่ ดูจะเป็นการหาทำไปหน่อยเยอะ

เราเปรียบเทียบกับไฟล์รูปดีกว่า ถ้าใครที่ถ่ายรูป ก็น่าจะคุ้นเคยกับพวก Raw File และ JPG ที่เราใช้งานกัน ในการเก็บข้อมูลพวกนี้ FASTQ ก็จะเหมือนกับ Raw File ที่มีขนาดใหญ่มาก ๆ มีข้อมูลของรูปภาพเยอะมาก ๆ ในขณะที่บางครั้งเราไม่ต้องการข้อมูลเยอะขนาดนั้น เราก็สามารถ Encode เป็น JPG ได้ เพื่อให้เราทำงานได้ง่ายขึ้นนั่นเอง

สรุป

จริง ๆ แล้วในการเก็บข้อมูลทาง Genomics เรามีหลาย Format ในการจัดเก็บ แต่ตัวที่เราใช้กันบ่อย ๆ ใช้กันเป็นหลัก ก็จะเป็นทั้ง FASTQ และ FASTA นอกจากนั้น มันก็ยังมีพวก BAM และ SAM ที่เราใช้เก็บพวก Alignment ต่าง ๆ อะยากไปมั้ยเนี่ย และพวก VCF ที่เก็บพวก Variance ที่เราหาออกมา มันก็จะมี Format สำหรับงาน และ Analysis ที่แตกต่างกันไป ในการทำงานจริง ๆ เราก็อาจจะต้องมีการ Process อะไรพวกนี้วน ๆ อยู่แถว ๆ นี้แหละ

ปล. เราจะบอกว่า Script ที่เราเขียนให้ดูสำหรับอ่าน FASTQ และ FASTA น่ะ มันมีคนเขียนมาเป็น Library แล้วนะชื่อว่า BioPython แต่ก่อนที่เราจะใช้ เราก็ควรจะเขียนเองให้ได้ก่อนใช่ม้าาา (พูดกลบเกลือนไปงั้นแหละ)