บทนำที่ครอบคลุมเกี่ยวกับจีโนมของคุณด้วย SciPy Stack

เผยแพร่แล้ว: 2022-03-11

ชีวสารสนเทศเป็นสาขาสหวิทยาการที่พัฒนาวิธีการและเครื่องมือซอฟต์แวร์สำหรับการวิเคราะห์และทำความเข้าใจข้อมูลทางชีววิทยา

กล่าวอย่างง่าย ๆ ว่าคุณสามารถคิดว่ามันเป็นวิทยาศาสตร์ข้อมูลสำหรับชีววิทยา

ในบรรดาข้อมูลทางชีววิทยาหลายประเภท ข้อมูลจีโนมเป็นหนึ่งในข้อมูลที่มีการวิเคราะห์อย่างกว้างขวางที่สุด โดยเฉพาะอย่างยิ่งกับความก้าวหน้าอย่างรวดเร็วของเทคโนโลยีการจัดลำดับดีเอ็นเอ (NGS) ยุคหน้า ปริมาณข้อมูลจีโนมได้เพิ่มขึ้นอย่างทวีคูณ จากข้อมูลของ Stephens Zachary D et al. การได้มาซึ่งข้อมูลจีโนมอยู่ในระดับเอกซาไบต์ต่อปี

บทนำที่ครอบคลุมเกี่ยวกับจีโนมของคุณด้วย SciPy

SciPy รวบรวมโมดูล Python จำนวนมากสำหรับการคำนวณทางวิทยาศาสตร์ ซึ่งเหมาะสำหรับความต้องการด้านชีวสารสนเทศจำนวนมาก

ทวีต

ในโพสต์นี้ ฉันสาธิตตัวอย่างการวิเคราะห์ไฟล์ GFF3 สำหรับจีโนมมนุษย์ด้วย SciPy Stack รูปแบบคุณลักษณะทั่วไป เวอร์ชัน 3 (GFF3) คือรูปแบบไฟล์ข้อความมาตรฐานปัจจุบันสำหรับการจัดเก็บคุณลักษณะจีโนม โดยเฉพาะอย่างยิ่ง ในโพสต์นี้ คุณจะได้เรียนรู้วิธีใช้ SciPy stack เพื่อตอบคำถามต่อไปนี้เกี่ยวกับจีโนมมนุษย์:

จีโนมไม่สมบูรณ์เท่าไร?
ในจีโนมมีกี่ยีน?
ยีนทั่วไปนานแค่ไหน?
การกระจายยีนระหว่างโครโมโซมมีลักษณะอย่างไร

ไฟล์ GFF3 ล่าสุดสำหรับจีโนมมนุษย์สามารถดาวน์โหลดได้จากที่นี่ ไฟล์ README ที่มาในไดเร็กทอรีนี้ให้คำอธิบายสั้นๆ เกี่ยวกับรูปแบบข้อมูลนี้ และพบข้อมูลจำเพาะที่ละเอียดยิ่งขึ้นที่นี่

เราจะใช้ Pandas ซึ่งเป็นองค์ประกอบหลักของ SciPy stack ที่ให้โครงสร้างข้อมูลที่รวดเร็ว ยืดหยุ่น และแสดงออก เพื่อจัดการและทำความเข้าใจไฟล์ GFF3

ติดตั้ง

อันดับแรก มาตั้งค่าสภาพแวดล้อมเสมือนด้วยการติดตั้ง SciPy stack กระบวนการนี้อาจใช้เวลานานหากสร้างจากแหล่งที่มาด้วยตนเอง เนื่องจากสแต็กเกี่ยวข้องกับแพ็คเกจจำนวนมาก ซึ่งบางส่วนขึ้นอยู่กับรหัส FORTRAN หรือ C ภายนอก ที่นี่ ฉันแนะนำให้ใช้ Miniconda ซึ่งทำให้การตั้งค่าง่ายมาก

 wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b

แฟล็ก -b บน bash line บอกให้ดำเนินการในโหมดแบตช์ หลังจากที่ใช้คำสั่งด้านบนเพื่อติดตั้ง Miniconda สำเร็จแล้ว ให้เริ่มสภาพแวดล้อมเสมือนใหม่สำหรับจีโนม แล้วติดตั้งกอง SciPy

 mkdir -p genomics cd genomics conda create -p venv ipython matplotlib pandas

โปรดทราบว่าเราได้ระบุเฉพาะ 3 แพ็คเกจที่เราจะใช้ในโพสต์นี้ หากคุณต้องการแพ็คเกจทั้งหมดที่อยู่ในสแต็ก SciPy ให้ผนวกแพ็คเกจเหล่านั้นต่อท้ายคำสั่ง conda create

หากคุณไม่แน่ใจในชื่อแพ็กเกจที่แน่นอน ให้ลอง conda search มาเปิดใช้งานสภาพแวดล้อมเสมือนและเริ่ม IPython

 source activate venv/ ipython

IPython เป็นการแทนที่ที่มีประสิทธิภาพมากขึ้นอย่างมากสำหรับอินเทอร์เฟซล่าม Python เริ่มต้น ดังนั้นสิ่งที่คุณเคยทำในล่าม python เริ่มต้นก็สามารถทำได้ใน IPython ฉันขอแนะนำโปรแกรมเมอร์ Python ทุกคนที่ยังไม่เคยใช้ IPython ให้ลองใช้เลย

ดาวน์โหลดไฟล์คำอธิบายประกอบ

เมื่อการตั้งค่าของเราเสร็จสิ้นแล้ว เรามาดาวน์โหลดไฟล์คำอธิบายประกอบจีโนมมนุษย์ในรูปแบบ GFF3

มีขนาดประมาณ 37 MB ซึ่งเป็นไฟล์ขนาดเล็กมากเมื่อเทียบกับเนื้อหาข้อมูลของจีโนมมนุษย์ ซึ่งมีขนาดประมาณ 3 GB ในข้อความธรรมดา นั่นเป็นเพราะไฟล์ GFF3 มีเพียงคำอธิบายประกอบของลำดับ ในขณะที่ข้อมูลลำดับมักจะจัดเก็บในรูปแบบไฟล์อื่นที่เรียกว่า FASTA หากคุณสนใจ คุณสามารถดาวน์โหลด FASTA ได้ที่นี่ แต่เราจะไม่ใช้ข้อมูลลำดับในบทช่วยสอนนี้

 !wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz

คำนำหน้า ! บอก IPython ว่านี่เป็นคำสั่งเชลล์แทนที่จะเป็นคำสั่ง Python อย่างไรก็ตาม IPython ยังสามารถประมวลผลคำสั่งเชลล์ที่ใช้บ่อยได้เช่น ls , pwd , rm , mkdir , rmdir แม้จะไม่มีคำนำหน้า ! .

เมื่อดูที่ส่วนหัวของไฟล์ GFF คุณจะเห็นบรรทัด metadata/pragmas/directives จำนวนมากขึ้นต้นด้วย ## หรือ #! .

ตาม README ## หมายความว่าข้อมูลเมตามีเสถียรภาพในขณะที่ #! หมายความว่ามันเป็นการทดลอง

ต่อมาคุณจะเห็น ### ซึ่งเป็นคำสั่งอื่นที่มีความหมายที่ลึกซึ้งยิ่งขึ้นตามข้อกำหนด

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

 ##gff-version 3 ##sequence-region 1 1 248956422 ##sequence-region 10 1 133797422 ##sequence-region 11 1 135086622 ##sequence-region 12 1 133275309 ... ##sequence-region MT 1 16569 ##sequence-region X 1 156040895 ##sequence-region Y 2781480 56887902 #!genome-build GRCh38.p7 #!genome-version GRCh38 #!genome-date 2013-12 #!genome-build-accession NCBI:GCA_000001405.22 #!genebuild-last-updated 2016-06

บรรทัดแรกระบุว่าเวอร์ชันของรูปแบบ GFF ที่ใช้ในไฟล์นี้คือ 3

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

บรรทัดที่ขึ้นต้นด้วย #! แสดงข้อมูลเกี่ยวกับโครงสร้างเฉพาะของจีโนม GRCh38.p7 ที่ไฟล์คำอธิบายประกอบนี้นำไปใช้

Genome Reference Consortium (GCR) เป็นสมาคมระหว่างประเทศที่ดูแลการปรับปรุงและปรับปรุงชุดจีโนมอ้างอิงหลายชุด รวมถึงกลุ่มสำหรับมนุษย์ หนูเมาส์ ปลาม้าลาย และไก่

เมื่อสแกนไฟล์นี้ ต่อไปนี้คือบรรทัดคำอธิบายประกอบสองสามบรรทัดแรก

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11 ### 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 1 . biological_region 10678 10687 0.999 + . logic_name=eponine 1 . biological_region 10681 10688 0.999 - . logic_name=eponine ...

คอลัมน์เป็น seqid , source , type , start , end , score , strand , phase , คุณลักษณะ บางส่วนของพวกเขาเข้าใจง่ายมาก ใช้บรรทัดแรกเป็นตัวอย่าง:

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11

นี่คือคำอธิบายประกอบของโครโมโซมแรกที่มีความต่อเนื่องของ 1 ซึ่งเริ่มจากฐานแรกถึงฐานที่ 24,895,622

กล่าวอีกนัยหนึ่งโครโมโซมแรกมีความยาวประมาณ 25 ล้านเบส

การวิเคราะห์ของเราจะไม่ต้องการข้อมูลจากสามคอลัมน์ . มีค่า (เช่น สกอร์ เกลียว และเฟส) ดังนั้นเราจึงสามารถเพิกเฉยได้ในตอนนี้

คอลัมน์แอตทริบิวต์สุดท้ายระบุว่าโครโมโซม 1 ยังมีชื่อนามแฝงสามชื่อ ได้แก่ CM000663.2, chr1 และ NC_000001.11 นั่นคือลักษณะของไฟล์ GFF3 โดยพื้นฐานแล้ว แต่เราจะไม่ตรวจสอบทีละบรรทัด ดังนั้นจึงถึงเวลาโหลดไฟล์ทั้งหมดลงใน Pandas

Pandas เหมาะสมดีสำหรับการจัดการกับรูปแบบ GFF3 เนื่องจากเป็นไฟล์ที่คั่นด้วยแท็บ และ Pandas รองรับการอ่านไฟล์ที่เหมือน CSV ได้เป็นอย่างดี

หมายเหตุ ข้อยกเว้นประการหนึ่งสำหรับรูปแบบที่คั่นด้วยแท็บคือเมื่อ GFF3 มี ##FASTA

ตามข้อกำหนด ##FASTA ระบุจุดสิ้นสุดของส่วนหมายเหตุประกอบ ซึ่งจะถูกตามด้วยหนึ่งลำดับขึ้นไปในรูปแบบ FASTA (ที่ไม่ใช่แท็บคั่น) แต่นี่ไม่ใช่กรณีสำหรับไฟล์ GFF3 ที่เราจะวิเคราะห์

 In [1]: import pandas as pd In [2]: pd.__version__ Out[2]: '0.18.1' In [3]: col_names = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'] Out[3]: df = pd.read_csv('Homo_sapiens.GRCh38.85.gff3.gz', compression='gzip', sep='\t', comment='#', low_memory=False, header=None, names=col_names)

บรรทัดสุดท้ายด้านบนโหลดไฟล์ GFF3 ทั้งหมดด้วยเมธอด pandas.read_csv

เนื่องจากไม่ใช่ไฟล์ CSV มาตรฐาน เราจึงต้องปรับแต่งการโทรเล็กน้อย

อันดับแรก เราแจ้ง Pandas เกี่ยวกับความไม่พร้อมใช้งานของข้อมูลส่วนหัวใน GFF3 ด้วย header=None จากนั้นเราระบุชื่อที่แน่นอนสำหรับแต่ละคอลัมน์ด้วย names=col_names

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

sep='\t' บอก Pandas ว่าคอลัมน์ต่างๆ ถูกคั่นด้วยแท็บแทนที่จะคั่นด้วยเครื่องหมายจุลภาค ค่าของ sep= สามารถเป็นนิพจน์ทั่วไปได้ (regex) สิ่งนี้จะมีประโยชน์หากไฟล์ในมือใช้ตัวคั่นที่แตกต่างกันสำหรับแต่ละคอลัมน์ (ใช่แล้ว มันเกิดขึ้น) comment='#' หมายถึงบรรทัดที่ขึ้นต้นด้วย # ถือเป็นความคิดเห็นและจะถูกละเว้น

compression='gzip' บอก Pandas ว่าไฟล์อินพุตเป็นแบบบีบอัด gzip

นอกจากนี้ pandas.read_csv ยังมีชุดพารามิเตอร์มากมายที่ช่วยให้อ่านรูปแบบไฟล์ที่เหมือน CSV ประเภทต่างๆ ได้

ประเภทของค่าที่ส่งคืนคือ DataFrame ซึ่งเป็นโครงสร้างข้อมูลที่สำคัญที่สุดใน Pandas ใช้สำหรับแสดงข้อมูล 2D

Pandas ยังมีโครงสร้างข้อมูล Series และ Panel สำหรับข้อมูล 1D และ 3D ตามลำดับ โปรดดูเอกสารประกอบสำหรับการแนะนำโครงสร้างข้อมูลของ Pandas

มาดูรายการสองสามรายการแรกที่ใช้วิธี . .head

 In [18]: df.head() Out[18]: seqid source type start end score strand phase attributes 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 1 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 2 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 3 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 4 1 . biological_region 10678 10687 0.999 + . logic_name=eponine

เอาต์พุตถูกจัดรูปแบบอย่างสวยงามในรูปแบบตารางโดยมีสตริงที่ยาวกว่าในคอลัมน์แอตทริบิวต์ถูกแทนที่ด้วย ... บางส่วน

คุณสามารถตั้งค่า Pandas ไม่ให้ละเว้นสตริงที่ยาวด้วย pd.set_option('display.max_colwidth', -1) นอกจากนี้ Pandas ยังมีตัวเลือกมากมายที่สามารถปรับแต่งได้

ต่อไป มาดูข้อมูลพื้นฐานเกี่ยวกับ dataframe นี้ด้วยวิธี .info

 In [20]: df.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 2601849 entries, 0 to 2601848 Data columns (total 9 columns): seqid object source object type object start int64 end int64 score object strand object phase object attributes object dtypes: int64(2), object(7) memory usage: 178.7+ MB

นี่แสดงว่า GFF3 มี 2,601,848 บรรทัดที่มีคำอธิบายประกอบ และแต่ละบรรทัดมีเก้าคอลัมน์

สำหรับแต่ละคอลัมน์ จะแสดงประเภทข้อมูลด้วย

จุดเริ่มต้นและจุดสิ้นสุดนั้นเป็นประเภท int64 ซึ่งเป็นจำนวนเต็มแสดงถึงตำแหน่งในจีโนม

คอลัมน์อื่นๆ เป็นประเภท object ทั้งหมด ซึ่งอาจหมายความว่าค่าของคอลัมน์ประกอบด้วยจำนวนเต็ม ทศนิยม และสตริงผสมกัน

ขนาดของข้อมูลทั้งหมดอยู่ที่ประมาณ 178.7+ MB ที่เก็บไว้ในหน่วยความจำ ปรากฏว่ามีขนาดเล็กกว่าไฟล์ที่ไม่มีการบีบอัด ซึ่งจะมีขนาดประมาณ 402 MB การตรวจสอบอย่างรวดเร็วแสดงอยู่ด้านล่าง

 gunzip -c Homo_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && du -s /tmp/tmp.gff3 402M /tmp/tmp.gff3

จากมุมมองระดับสูง เราได้โหลดไฟล์ GFF3 ทั้งหมดลงในวัตถุ DataFrame ใน Python และการวิเคราะห์ต่อไปนี้ทั้งหมดของเราจะขึ้นอยู่กับวัตถุเดียวนี้

ทีนี้ มาดูกันว่าลำดับของคอลัมน์แรกเกี่ยวกับอะไร

 In [29]: df.seqid.unique() Out[29]: array(['1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6', '7', '8', '9', 'GL000008.2', 'GL000009.2', 'GL000194.1', 'GL000195.1', ... 'KI270757.1', 'MT', 'X', 'Y'], dtype=object) In [30]: df.seqid.unique().shape Out[30]: (194,)

df.seqid เป็นวิธีหนึ่งในการเข้าถึงข้อมูลคอลัมน์จากดาต้าเฟรม อีกวิธีหนึ่งคือ df['seqid'] ซึ่งเป็นรูปแบบทั่วไปมากกว่า เพราะหากชื่อคอลัมน์เป็นคำสำคัญที่สงวนไว้ของ Python (เช่น class ) หรือมี . หรืออักขระเว้นวรรค วิธีแรก ( df.seqid ) จะไม่ทำงาน

ผลลัพธ์แสดงให้เห็นว่ามีซีคิดที่ไม่ซ้ำกัน 194 ตัว ซึ่งรวมถึงโครโมโซม 1 ถึง 22, X, Y และไมโตคอนเดรีย (MT) DNA และซีคิดอื่น ๆ อีก 169 รายการ

ลำดับขั้นที่ขึ้นต้นด้วย KI และ GL เป็นลำดับดีเอ็นเอหรือโครงนั่งร้านในจีโนมที่ยังประกอบไม่สำเร็จในโครโมโซม

สำหรับผู้ที่ไม่คุ้นเคยกับจีโนม นี่เป็นสิ่งสำคัญ

แม้ว่าจีโนมมนุษย์ร่างแรกจะออกมาเมื่อกว่า 15 ปีที่แล้ว แต่จีโนมมนุษย์ในปัจจุบันก็ยังไม่สมบูรณ์ ความยากลำบากในการประกอบลำดับเหล่านี้ส่วนใหญ่เกิดจากบริเวณที่ซ้ำซ้อนที่ซับซ้อนในจีโนม

ต่อไป มาดูคอลัมน์ต้นทางกัน

README กล่าวว่าแหล่งที่มาเป็นตัวระบุข้อความอิสระที่มีจุดประสงค์เพื่ออธิบายอัลกอริทึมหรือขั้นตอนการปฏิบัติงานที่สร้างคุณลักษณะนี้

 In [66]: df.source.value_counts() Out[66]: havana 1441093 ensembl_havana 745065 ensembl 228212 . 182510 mirbase 4701 GRCh38 194 insdc 74

นี่คือตัวอย่างการใช้เมธอด value_counts ซึ่งมีประโยชน์อย่างมากสำหรับการนับตัวแปรหมวดหมู่อย่างรวดเร็ว

จากผลลัพธ์ เราจะเห็นว่ามีค่าที่เป็นไปได้เจ็ดค่าสำหรับคอลัมน์นี้ และรายการส่วนใหญ่ในไฟล์ GFF3 มาจาก havana, ensembl และ ensembl_havana

คุณสามารถเรียนรู้เพิ่มเติมเกี่ยวกับความหมายของแหล่งข้อมูลเหล่านี้และความสัมพันธ์ระหว่างแหล่งข้อมูลเหล่านี้ได้ในโพสต์นี้

เพื่อให้ง่ายขึ้น เราจะเน้นที่รายการจากแหล่งที่มา GRCh38, ฮาวานา, วงดนตรี และ ensembl_havan.a

กี่จีโนมไม่สมบูรณ์?

ข้อมูลเกี่ยวกับโครโมโซมทั้งหมดอยู่ในรายการจากแหล่งที่มา GRCh38 ดังนั้น ก่อนอื่นให้กรองส่วนที่เหลือออก และกำหนดผลลัพธ์ที่กรองแล้วให้กับตัวแปรใหม่ gdf

 In [70]: gdf = df[df.source == 'GRCh38'] In [87]: gdf.shape Out[87]: (194, 9) In [84]: gdf.sample(10) Out[84]: seqid source type start end score strand phase attributes 2511585 KI270708.1 GRCh38 supercontig 1 127682 . . . ID=supercontig:KI270708.1;Alias=chr1_KI270708v... 2510840 GL000208.1 GRCh38 supercontig 1 92689 . . . ID=supercontig:GL000208.1;Alias=chr5_GL000208v... 990810 17 GRCh38 chromosome 1 83257441 . . . ID=chromosome:17;Alias=CM000679.2,chr17,NC_000... 2511481 KI270373.1 GRCh38 supercontig 1 1451 . . . ID=supercontig:KI270373.1;Alias=chrUn_KI270373... 2511490 KI270384.1 GRCh38 supercontig 1 1658 . . . ID=supercontig:KI270384.1;Alias=chrUn_KI270384... 2080148 6 GRCh38 chromosome 1 170805979 . . . ID=chromosome:6;Alias=CM000668.2,chr6,NC_00000... 2511504 KI270412.1 GRCh38 supercontig 1 1179 . . . ID=supercontig:KI270412.1;Alias=chrUn_KI270412... 1201561 19 GRCh38 chromosome 1 58617616 . . . ID=chromosome:19;Alias=CM000681.2,chr19,NC_000... 2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID=supercontig:KI270340.1;Alias=chrUn_KI270340... 2594560 Y GRCh38 chromosome 2781480 56887902 . . . ID=chromosome:Y;Alias=CM000686.2,chrY,NC_00002...

การกรองเป็นเรื่องง่ายใน Pandas

หากคุณตรวจสอบค่าที่ประเมินจากนิพจน์ df.source == 'GRCh38' จะเป็นชุดค่า True และ False สำหรับแต่ละรายการที่มีดัชนีเดียวกันกับ df การส่งผ่านไปยัง df[] จะส่งคืนเฉพาะรายการเหล่านั้นโดยที่ค่าที่สอดคล้องกันเป็น True

มี 194 คีย์ใน df[] ซึ่ง df.source == 'GRCh38'

ดังที่เราได้เห็นก่อนหน้านี้ ยังมีค่าที่ไม่ซ้ำกัน 194 ค่าในคอลัมน์ seqid ซึ่งหมายความว่าแต่ละรายการใน gdf จะสอดคล้องกับ seqid ที่เจาะจง

จากนั้นเราสุ่มเลือก 10 รายการด้วยวิธี sample เพื่อดูรายละเอียด

คุณจะเห็นว่าลำดับที่ไม่ประกอบเป็นประเภท supercontig ในขณะที่ลำดับอื่นๆ เป็นโครโมโซม ในการคำนวณเศษส่วนของจีโนมที่ไม่สมบูรณ์ ก่อนอื่นเราต้องรู้ความยาวของจีโนมทั้งหมด ซึ่งก็คือผลรวมของความยาวของซีกิดทั้งหมด

 In [90]: gdf = gdf.copy() In [91]: gdf['length'] = gdf.end - gdf.start + 1 In [93]: gdf.head() Out[93]: seqid source type start end score strand phase attributes length 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 248956421 235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2,chr10,NC_000... 133797421 328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2,chr11,NC_000... 135086621 483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2,chr12,NC_000... 133275308 634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2,chr13,NC_000... 114364327 In [97]: gdf.length.sum() Out[97]: 3096629532 In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT'] In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum() Out[101]: 0.0037021917421198327

ในตัวอย่างด้านบนก่อน เราทำสำเนาของ gdf ด้วย . .copy() มิฉะนั้น gdf ดั้งเดิมเป็นเพียงส่วนหนึ่งของ df และการแก้ไขโดยตรงจะส่งผลให้ SettingWithCopyWarning (ดูรายละเอียดเพิ่มเติมที่นี่)

จากนั้นเราจะคำนวณความยาวของแต่ละรายการและเพิ่มกลับเข้าไปใน gdf เป็นคอลัมน์ใหม่ที่ชื่อ "ความยาว" ความยาวรวมจะอยู่ที่ประมาณ 3.1 พันล้าน และเศษส่วนของลำดับที่ไม่ได้ประกอบคือประมาณ 0.37%

นี่คือวิธีการแบ่งส่วนในสองคำสั่งสุดท้าย

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

เครื่องหมายลบ ( - ) ถูกเพิ่มที่จุดเริ่มต้นของดัชนีเพื่อย้อนกลับการเลือก เนื่องจากเราต้องการทุกอย่างที่ ไม่ อยู่ในรายการ (เช่น เราต้องการสิ่งที่ไม่ได้ประกอบที่ขึ้นต้นด้วย KI และ GL)...

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

 gdf[(gdf['type'] == 'supercontig')].length.sum() / gdf.length.sum()

มียีนกี่ตัว?

ที่นี่เราเน้นที่รายการจากแหล่งรวม, ฮาวานา และ ensembl_havana เนื่องจากเป็นที่ที่รายการคำอธิบายประกอบส่วนใหญ่อยู่

 In [109]: edf = df[df.source.isin(['ensembl', 'havana', 'ensembl_havana'])] In [111]: edf.sample(10) Out[111]: seqid source type start end score strand phase attributes 915996 16 havana CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0... 2531429 X havana exon 41196251 41196359 . + . Parent=transcript:ENST00000462850;Name=ENSE000... 1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0... 243070 10 havana exon 13116267 13116340 . + . Parent=transcript:ENST00000378764;Name=ENSE000... 2413583 8 ensembl_havana exon 144359184 144359423 . + . Parent=transcript:ENST00000530047;Name=ENSE000... 2160496 6 havana exon 111322569 111322678 . - . Parent=transcript:ENST00000434009;Name=ENSE000... 839952 15 havana exon 76227713 76227897 . - . Parent=transcript:ENST00000565910;Name=ENSE000... 957782 16 ensembl_havana exon 67541653 67541782 . + . Parent=transcript:ENST00000379312;Name=ENSE000... 1632979 21 ensembl_havana exon 37840658 37840709 . - . Parent=transcript:ENST00000609713;Name=ENSE000... 1953399 4 havana exon 165464390 165464586 . + . Parent=transcript:ENST00000511992;Name=ENSE000... In [123]: edf.type.value_counts() Out[123]: exon 1180596 CDS 704604 five_prime_UTR 142387 three_prime_UTR 133938 transcript 96375 gene 42470 processed_transcript 28228 ... Name: type, dtype: int64

ใช้วิธี isin อีกครั้งสำหรับการกรอง

จากนั้น การนับค่าอย่างรวดเร็วแสดงว่ารายการส่วนใหญ่เป็น exon, coding sequence (CDS) และขอบเขตที่ไม่ได้แปล (UTR)

สิ่งเหล่านี้เป็นองค์ประกอบของยีนย่อย แต่เรากำลังมองหาการนับจำนวนยีนเป็นหลัก ดังที่แสดงว่ามี 42,470 แต่เราต้องการทราบข้อมูลเพิ่มเติม

โดยเฉพาะพวกเขาชื่ออะไรและพวกเขาทำอะไร? ในการตอบคำถามเหล่านี้ เราต้องดูข้อมูลในคอลัมน์แอตทริบิวต์อย่างละเอียด

 In [127]: ndf = edf[edf.type == 'gene'] In [173]: ndf = ndf.copy() In [133]: ndf.sample(10).attributes.values Out[133]: array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=hepatocyte nuclear factor 4 gamma pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:35417];gene_id=ENSG00000228611;havana_gene=OTTHUMG00000016986;havana_version=2;logic_name=havana;version=2', 'ID=gene:ENSG00000177189;Name=RPS6KA3;biotype=protein_coding;description=ribosomal protein S6 kinase A3 [Source:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12', 'ID=gene:ENSG00000231748;Name=RP11-227H15.5;biotype=antisense;gene_id=ENSG00000231748;havana_gene=OTTHUMG00000018373;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1 receptor 33 pseudogene [Source:HGNC Symbol%3BAcc:HGNC:37353];gene_id=ENSG00000227426;havana_gene=OTTHUMG00000154474;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000087250;Name=MT3;biotype=protein_coding;description=metallothionein 3 [Source:HGNC Symbol%3BAcc:HGNC:7408];gene_id=ENSG00000087250;havana_gene=OTTHUMG00000133282;havana_version=3;logic_name=ensembl_havana_gene;version=8', 'ID=gene:ENSG00000177108;Name=ZDHHC22;biotype=protein_coding;description=zinc finger DHHC-type containing 22 [Source:HGNC Symbol%3BAcc:HGNC:20106];gene_id=ENSG00000177108;havana_gene=OTTHUMG00000171575;havana_version=3;logic_name=ensembl_havana_gene;version=5', 'ID=gene:ENSG00000249784;Name=SCARNA22;biotype=scaRNA;description=small Cajal body-specific RNA 22 [Source:HGNC Symbol%3BAcc:HGNC:32580];gene_id=ENSG00000249784;logic_name=ncrna;version=1', 'ID=gene:ENSG00000079101;Name=CLUL1;biotype=protein_coding;description=clusterin like 1 [Source:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16', 'ID=gene:ENSG00000229224;Name=AC105398.3;biotype=antisense;gene_id=ENSG00000229224;havana_gene=OTTHUMG00000152025;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=lymphocyte antigen 6 complex%2C locus G6E (pseudogene) [Source:HGNC Symbol%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype=object)

จัดรูปแบบเป็นรายการคู่ค่าแท็กที่คั่นด้วยเครื่องหมายอัฒภาค ข้อมูลที่เราสนใจมากที่สุดคือชื่อยีน รหัสยีน และคำอธิบาย และเราจะดึงข้อมูลเหล่านี้ออกมาด้วยนิพจน์ทั่วไป (regex)

 import re RE_GENE_NAME = re.compile(r'Name=(?P<gene_name>.+?);') def extract_gene_name(attributes_str): res = RE_GENE_NAME.search(attributes_str) return res.group('gene_name') ndf['gene_name'] = ndf.attributes.apply(extract_gene_name)

ขั้นแรก เราแยกชื่อยีน

ใน regex Name=(?P<gene_name>.+?); , +? ใช้แทน + เพราะเราต้องการให้ไม่โลภและให้การค้นหาหยุดที่เครื่องหมายอัฒภาคแรก มิฉะนั้น ผลลัพธ์จะจับคู่กับอัฒภาคสุดท้าย

นอกจากนี้ regex ยังถูกคอมไพล์ด้วย re.compile แทนที่จะใช้โดยตรงใน re.search เพื่อประสิทธิภาพที่ดีขึ้น เพราะเราจะนำไปใช้กับสตริงแอตทริบิวต์นับพัน

extract_gene_name ทำหน้าที่เป็นฟังก์ชันตัวช่วยที่ใช้ใน pd.apply ซึ่งเป็นวิธีที่จะใช้เมื่อจำเป็นต้องใช้ฟังก์ชันกับทุกรายการของ dataframe หรือ series

ในกรณีนี้ เราต้องการแยกชื่อยีนของทุกรายการใน ndf.attributes และเพิ่มชื่อกลับเข้าไปใน ndf ในคอลัมน์ใหม่ชื่อ gene_name

รหัสยีนและคำอธิบายจะถูกดึงออกมาในลักษณะเดียวกัน

 RE_GENE_ID = re.compile(r'gene_id=(?P<gene_id>ENSG.+?);') def extract_gene_id(attributes_str): res = RE_GENE_ID.search(attributes_str) return res.group('gene_id') ndf['gene_id'] = ndf.attributes.apply(extract_gene_id) RE_DESC = re.compile('description=(?P<desc>.+?);') def extract_description(attributes_str): res = RE_DESC.search(attributes_str) if res is None: return '' else: return res.group('desc') ndf['desc'] = ndf.attributes.apply(extract_description)

regex สำหรับ RE_GENE_ID มีความเฉพาะเจาะจงมากกว่าเล็กน้อย เนื่องจากเรารู้ว่าทุก gene_id ต้องขึ้นต้นด้วย ENSG โดยที่ ENS หมายถึง ensembl และ G หมายถึง gene

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

 In [224]: ndf.drop('attributes', axis=1, inplace=True) In [225]: ndf.head() Out[225]: seqid source type start end score strand phase gene_id gene_name desc 16 1 havana gene 11869 14409 . + . ENSG00000223972 DDX11L1 DEAD/H-box helicase 11 like 1 [Source:HGNC Sym... 28 1 havana gene 14404 29570 . - . ENSG00000227232 WASH7P WAS protein family homolog 7 pseudogene [Sourc... 71 1 havana gene 52473 53312 . + . ENSG00000268020 OR4G4P olfactory receptor family 4 subfamily G member... 74 1 havana gene 62948 63887 . + . ENSG00000240361 OR4G11P olfactory receptor family 4 subfamily G member... 77 1 ensembl_havana gene 69091 70008 . + . ENSG00000186092 OR4F5 olfactory receptor family 4 subfamily F member...

ในการเรียกข้างต้น attributes ระบุคอลัมน์เฉพาะที่เราต้องการวาง

axis=1 หมายถึง เรากำลังวางคอลัมน์แทนที่จะเป็นแถว ( axis=0 โดยค่าเริ่มต้น)

inplace=True หมายความว่าการดรอปดำเนินการบน DataFrame เอง แทนที่จะส่งคืนสำเนาใหม่โดยมีการดร็อปคอลัมน์ที่ระบุ

การดู .head อย่างรวดเร็วแสดงให้เห็นว่าคอลัมน์แอตทริบิวต์หายไปแล้ว และมีการเพิ่มคอลัมน์ใหม่สามคอลัมน์: gene_name , gene_id และ desc

ด้วยความอยากรู้ มาดูกันว่า gene_id และ gene_name ทั้งหมดไม่ซ้ำกันหรือไม่:

 In [232]: ndf.shape Out[232]: (42470, 11) In [233]: ndf.gene_id.unique().shape Out[233]: (42470,) In [234]: ndf.gene_name.unique().shape Out[234]: (42387,)

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

 In [243]: count_df = ndf.groupby('gene_name').count().ix[:, 0].sort_values().ix[::-1] In [244]: count_df.head(10) Out[244]: gene_name SCARNA20 7 SCARNA16 6 SCARNA17 5 SCARNA15 4 SCARNA21 4 SCARNA11 4 Clostridiales-1 3 SCARNA4 3 C1QTNF9B-AS1 2 C11orf71 2 Name: seqid, dtype: int64 In [262]: count_df[count_df > 1].shape Out[262]: (63,) In [263]: count_df.shape Out[263]: (42387,) In [264]: count_df[count_df > 1].shape[0] / count_df.shape[0] Out[264]: 0.0014863047632528842

เราจัดกลุ่มรายการทั้งหมดตามค่าของ gene_name จากนั้นนับจำนวนรายการในแต่ละกลุ่มด้วย .count()

หากคุณตรวจสอบผลลัพธ์จาก ndf.groupby('gene_name').count() คอลัมน์ทั้งหมดจะถูกนับสำหรับแต่ละกลุ่ม แต่ส่วนใหญ่มีค่าเหมือนกัน

โปรดทราบว่าจะไม่มีการพิจารณาค่า NA ในการนับ ดังนั้นให้นับเฉพาะคอลัมน์แรก ตาม seqid (เราใช้ .ix[:, 0] เพื่อให้แน่ใจว่าไม่มีค่า NA)

จากนั้นเรียงลำดับค่าการนับด้วย .sort_values และกลับลำดับด้วย .ix[::-1]

ผลลัพธ์ที่ได้ สามารถแชร์ชื่อยีนกับรหัสยีนได้ถึงเจ็ดรหัส

 In [255]: ndf[ndf.gene_name == 'SCARNA20'] Out[255]: seqid source type start end score strand phase gene_id gene_name desc 179399 1 ensembl gene 171768070 171768175 . + . ENSG00000253060 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 201037 1 ensembl gene 204727991 204728106 . + . ENSG00000251861 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 349203 11 ensembl gene 8555016 8555146 . + . ENSG00000252778 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 718520 14 ensembl gene 63479272 63479413 . + . ENSG00000252800 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 837233 15 ensembl gene 75121536 75121666 . - . ENSG00000252722 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1039874 17 ensembl gene 28018770 28018907 . + . ENSG00000251818 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1108215 17 ensembl gene 60231516 60231646 . - . ENSG00000252577 SCARNA20 small Cajal body-specific RNA 20 [Source:HGNC Symbol%3BAcc:HGNC:32578]

การพิจารณายีน SCARNA20 อย่างละเอียดถี่ถ้วนแสดงให้เห็นว่าพวกมันต่างกันโดยสิ้นเชิง

แม้ว่าพวกมันจะใช้ชื่อเดียวกัน แต่พวกมันจะอยู่ที่ตำแหน่งต่างๆ ของจีโนม

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

ประเด็นคือต้องรู้ว่าชื่อยีนนั้นไม่ซ้ำกันสำหรับรหัสยีนทั้งหมด และประมาณ 0.15% ของชื่อที่แชร์โดยยีนหลายตัว

ยีนทั่วไปนานแค่ไหน?

คล้ายกับที่เราทำเมื่อเราพยายามทำความเข้าใจความไม่สมบูรณ์ของจีโนม เราสามารถเพิ่มคอลัมน์ length ลงใน ndf ได้อย่างง่ายดาย:

 In [277]: ndf['length'] = ndf.end - ndf.start + 1 In [278]: ndf.length.describe() Out[278]: count 4.247000e+04 mean 3.583348e+04 std 9.683485e+04 min 8.000000e+00 25% 8.840000e+02 50% 5.170500e+03 75% 3.055200e+04 max 2.304997e+06 Name: length, dtype: float64

.describe() คำนวณสถิติง่ายๆ ตามค่าความยาว:

ความยาวเฉลี่ยของยีนประมาณ 36,000 เบส
ความยาวมัธยฐานของยีนยาวประมาณ 5,200 เบส
ความยาวยีนต่ำสุดและสูงสุดคือประมาณแปดและ 2.3 ล้านเบสตามลำดับ

เนื่องจากค่าเฉลี่ยมากกว่าค่ามัธยฐานมาก มันจึงหมายความว่าการกระจายความยาวเบ้ไปทางขวา เพื่อให้มีลักษณะที่เป็นรูปธรรมมากขึ้น ให้พลอตการแจกแจงกัน

 import matplotlib as plt ndf.length.plot(kind='hist', bins=50, logy=True) plt.show()

Pandas มีอินเทอร์เฟซที่เรียบง่ายสำหรับ matplotlib เพื่อให้การวางแผนสะดวกยิ่งขึ้นด้วย DataFrames หรือซีรีส์

ในกรณีนี้ มันบอกว่าเราต้องการพล็อตฮิสโตแกรม ( kind='hist' ) ที่มีถังขยะ 50 ช่อง และให้แกน y อยู่ในมาตราส่วนบันทึก ( logy=True )

จากฮิสโตแกรม เราจะเห็นว่ายีนส่วนใหญ่อยู่ภายในถังแรก อย่างไรก็ตาม ความยาวของยีนบางตัวอาจมีมากกว่าสองล้านเบส มาดูกันว่าพวกเขาคืออะไร:

 In [39]: ndf[ndf.length > 2e6].sort_values('length').ix[::-1] Out[39]: seqid source type start end score strand phase gene_name gene_id desc length 2309345 7 ensembl_havana gene 146116002 148420998 . + . CNTNAP2 ENSG00000174469 contactin associated protein-like 2 [Source:HG... 2304997 2422510 9 ensembl_havana gene 8314246 10612723 . - . PTPRD ENSG00000153707 protein tyrosine phosphatase%2C receptor type ... 2298478 2527169 X ensembl_havana gene 31097677 33339441 . - . DMD ENSG00000198947 dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928] 2241765 440886 11 ensembl_havana gene 83455012 85627922 . - . DLG2 ENSG00000150672 discs large MAGUK scaffold protein 2 [Source:H... 2172911 2323457 8 ensembl_havana gene 2935353 4994972 . - . CSMD1 ENSG00000183117 CUB and Sushi multiple domains 1 [Source:HGNC ... 2059620 1569914 20 ensembl_havana gene 13995369 16053197 . + . MACROD2 ENSG00000172264 MACRO domain containing 2 [Source:HGNC Symbol%... 2057829

อย่างที่คุณเห็น ยีนที่ยาวที่สุดชื่อ CNTNAP2 ซึ่งย่อมาจาก contactin ที่สัมพันธ์กับโปรตีนเหมือน 2 ตามหน้าวิกิพีเดีย

ยีนนี้ครอบคลุมเกือบ 1.6% ของโครโมโซม 7 และเป็นหนึ่งในยีนที่ใหญ่ที่สุดในจีโนมมนุษย์

อย่างแท้จริง! เราเพิ่งตรวจสอบว่าตัวเอง ในทางตรงกันข้าม ยีนที่เล็กที่สุดล่ะ? ปรากฎว่าพวกมันสั้นถึงแปดฐาน

 In [40]: ndf.sort_values('length').head() Out[40]: seqid source type start end score strand phase gene_name gene_id desc length 682278 14 havana gene 22438547 22438554 . + . TRDD1 ENSG00000223997 T cell receptor delta diversity 1 [Source:HGNC... 8 682282 14 havana gene 22439007 22439015 . + . TRDD2 ENSG00000237235 T cell receptor delta diversity 2 [Source:HGNC... 9 2306836 7 havana gene 142786213 142786224 . + . TRBD1 ENSG00000282431 T cell receptor beta diversity 1 [Source:HGNC ... 12 682286 14 havana gene 22449113 22449125 . + . TRDD3 ENSG00000228985 T cell receptor delta diversity 3 [Source:HGNC... 13 1879625 4 havana gene 10238213 10238235 . - . AC006499.9 ENSG00000271544 23

ความยาวของกรณีสุดโต่งทั้งสองกรณีนี้ห่างกันห้าลำดับความสำคัญ (2.3 ล้านเทียบกับ 8) ซึ่งถือว่าใหญ่มากและสามารถบ่งบอกถึงระดับความหลากหลายของชีวิตได้

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

การกระจายยีนระหว่างโครโมโซม

สิ่งสุดท้ายที่ฉันต้องการจะพูดถึงคือการกระจายยีนระหว่างโครโมโซม ซึ่งยังทำหน้าที่เป็นตัวอย่างสำหรับการแนะนำวิธีการ .merge เพื่อรวม .merge สองชุดเข้าด้วยกัน ตามสัญชาตญาณ โครโมโซมที่ยาวขึ้นมีแนวโน้มที่จะเป็นโฮสต์ของยีนมากกว่า มาดูกันว่าจริงหรือไม่

 In [53]: ndf = ndf[ndf.seqid.isin(chrs)] In [54]: chr_gene_counts = ndf.groupby('seqid').count().ix[:, 0].sort_values().ix[::-1] Out[54]: chr_gene_counts seqid 1 3902 2 2806 11 2561 19 2412 17 2280 3 2204 6 2154 12 2140 7 2106 5 2002 16 1881 X 1852 4 1751 9 1659 8 1628 10 1600 15 1476 14 1449 22 996 20 965 13 872 18 766 21 541 Y 436 Name: source, dtype: int64

เรายืมตัวแปร chrs จากส่วนก่อนหน้า และใช้เพื่อกรองลำดับที่ยังไม่ได้ประกอบ จากผลลัพธ์ โครโมโซม 1 ที่ใหญ่ที่สุดมียีนมากที่สุดแน่นอน แม้ว่าโครโมโซม Y จะมีจำนวนยีนที่น้อยที่สุด แต่ก็ไม่ใช่โครโมโซมที่เล็กที่สุด

โปรดทราบว่าดูเหมือนว่าจะไม่มียีนในไมโตคอนเดรีย (MT) ซึ่งไม่เป็นความจริง

การกรองเพิ่มเติมอีกเล็กน้อยใน DataFrame df แรกที่ส่งคืนโดย pd.read_csv แสดงว่ายีน MT ทั้งหมดมาจาก insdc แหล่งที่มา (ซึ่งถูกกรองออกก่อนหน้านี้เมื่อสร้าง edf โดยที่เราพิจารณาเฉพาะแหล่งที่มาของฮาวานา วงดนตรี หรือ ensembl_havana)

 In [134]: df[(df.type == 'gene') & (df.seqid == 'MT')] Out[134]: seqid source type start end score strand phase attributes 2514003 MT insdc gene 648 1601 . + . ID=gene:ENSG00000211459;Name=MT-RNR1;biotype=M... 2514009 MT insdc gene 1671 3229 . + . ID=gene:ENSG00000210082;Name=MT-RNR2;biotype=M... 2514016 MT insdc gene 3307 4262 . + . ID=gene:ENSG00000198888;Name=MT-ND1;biotype=pr... 2514029 MT insdc gene 4470 5511 . + . ID=gene:ENSG00000198763;Name=MT-ND2;biotype=pr... 2514048 MT insdc gene 5904 7445 . + . ID=gene:ENSG00000198804;Name=MT-CO1;biotype=pr... 2514058 MT insdc gene 7586 8269 . + . ID=gene:ENSG00000198712;Name=MT-CO2;biotype=pr... 2514065 MT insdc gene 8366 8572 . + . ID=gene:ENSG00000228253;Name=MT-ATP8;biotype=p... 2514069 MT insdc gene 8527 9207 . + . ID=gene:ENSG00000198899;Name=MT-ATP6;biotype=p... 2514073 MT insdc gene 9207 9990 . + . ID=gene:ENSG00000198938;Name=MT-CO3;biotype=pr... 2514080 MT insdc gene 10059 10404 . + . ID=gene:ENSG00000198840;Name=MT-ND3;biotype=pr... 2514087 MT insdc gene 10470 10766 . + . ID=gene:ENSG00000212907;Name=MT-ND4L;biotype=p... 2514091 MT insdc gene 10760 12137 . + . ID=gene:ENSG00000198886;Name=MT-ND4;biotype=pr... 2514104 MT insdc gene 12337 14148 . + . ID=gene:ENSG00000198786;Name=MT-ND5;biotype=pr... 2514108 MT insdc gene 14149 14673 . - . ID=gene:ENSG00000198695;Name=MT-ND6;biotype=pr... 2514115 MT insdc gene 14747 15887 . + . ID=gene:ENSG00000198727;Name=MT-CYB;biotype=pr...

This example also shows how to combine two conditions during filtering with & ; the logical operator for “or” would be | .

Note that the parentheses around each condition are required, and this part of the syntax in Pandas is different from Python, which would have been literal and and or .

Next, let's borrow the gdf DataFrame from the previous section as a source for the length of each chromosome:

 In [61]: gdf = gdf[gdf.seqid.isin(chrs)] In [62]: gdf.drop(['start', 'end', 'score', 'strand', 'phase' ,'attributes'], axis=1, inplace=True) In [63]: gdf.sort_values('length').ix[::-1] Out[63]: seqid source type length 0 1 GRCh38 chromosome 248956422 1364641 2 GRCh38 chromosome 242193529 1705855 3 GRCh38 chromosome 198295559 1864567 4 GRCh38 chromosome 190214555 1964921 5 GRCh38 chromosome 181538259 2080148 6 GRCh38 chromosome 170805979 2196981 7 GRCh38 chromosome 159345973 2514125 X GRCh38 chromosome 156040895 2321361 8 GRCh38 chromosome 145138636 2416560 9 GRCh38 chromosome 138394717 328938 11 GRCh38 chromosome 135086622 235068 10 GRCh38 chromosome 133797422 483370 12 GRCh38 chromosome 133275309 634486 13 GRCh38 chromosome 114364328 674767 14 GRCh38 chromosome 107043718 767312 15 GRCh38 chromosome 101991189 865053 16 GRCh38 chromosome 90338345 990810 17 GRCh38 chromosome 83257441 1155977 18 GRCh38 chromosome 80373285 1559144 20 GRCh38 chromosome 64444167 1201561 19 GRCh38 chromosome 58617616 2594560 Y GRCh38 chromosome 54106423 1647482 22 GRCh38 chromosome 50818468 1616710 21 GRCh38 chromosome 46709983 2513999 MT GRCh38 chromosome 16569

The columns that are not relevant to the analysis are dropped for clarity.

Yes, .drop can also take a list of columns and drop them altogether in one operation.

Note that the row with a seqid of MT is still there; we will get back to it later. The next operation we will perform is merge the two datasets based on the values of seqid.

 In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index() In [75]: cdf.head(2) Out[75]: seqid gene_count 0 1 3902 1 2 2806 In [78]: merged = gdf.merge(cdf, on='seqid') In [79]: merged Out[79]: seqid source type length gene_count 0 1 GRCh38 chromosome 248956422 3902 1 10 GRCh38 chromosome 133797422 1600 2 11 GRCh38 chromosome 135086622 2561 3 12 GRCh38 chromosome 133275309 2140 4 13 GRCh38 chromosome 114364328 872 5 14 GRCh38 chromosome 107043718 1449 6 15 GRCh38 chromosome 101991189 1476 7 16 GRCh38 chromosome 90338345 1881 8 17 GRCh38 chromosome 83257441 2280 9 18 GRCh38 chromosome 80373285 766 10 19 GRCh38 chromosome 58617616 2412 11 2 GRCh38 chromosome 242193529 2806 12 20 GRCh38 chromosome 64444167 965 13 21 GRCh38 chromosome 46709983 541 14 22 GRCh38 chromosome 50818468 996 15 3 GRCh38 chromosome 198295559 2204 16 4 GRCh38 chromosome 190214555 1751 17 5 GRCh38 chromosome 181538259 2002 18 6 GRCh38 chromosome 170805979 2154 19 7 GRCh38 chromosome 159345973 2106 20 8 GRCh38 chromosome 145138636 1628 21 9 GRCh38 chromosome 138394717 1659 22 X GRCh38 chromosome 156040895 1852 23 Y GRCh38 chromosome 54106423 436

Since chr_gene_counts is still a Series object, which doesn't support a merge operation, we need to convert it to a DataFrame object first with .to_frame .

.reset_index() converts the original index (ie seqid ) into a new column and resets current index as 0-based incremental numbers.

The output from cdf.head(2) shows what it looks like. Next, we used the .merge method to combine the two DataFrame on the seqid column ( on='seqid' ).

After merging gdf and cdf , the MT entry is still missing. This is because, by default, .merge operates an inner join, while left join, right join, or outer join are available by tuning the how parameter.

Please refer to the documentation for more details.

Later, you may find that there is also a related .join method. .merge and .join are similar yet have different APIs.

According to the official documentation says

The related DataFrame.join method, uses merge internally for the index-on-index and index-on-column(s) joins, but joins on indexes by default rather than trying to join on common columns (the default behavior for merge). If you are joining on index, you may wish to use DataFrame.join to save yourself some typing.

Basically, .merge is more general-purpose and used by .join .

Finally, we are ready to calculate the correlation between chromosome length and gene_count .

 In [81]: merged[['length', 'gene_count']].corr() Out[81]: length gene_count length 1.000000 0.728221 gene_count 0.728221 1.000000

By default .corr calculates the Pearson correlation between all pairs of columns in a dataframe.

But we have only a single pair of columns in this case, and the correlation turns out to be positive – 0.73.

In other words, the larger the chromosome, the more likely it is to have more genes. Let's also plot the two columns after sorting the value pairs by length :

 ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count',) # add some margin to both ends of x axis xlim = ax.get_xlim() margin = xlim[0] * 0.1 ax.set_xlim([xlim[0] - margin, xlim[1] + margin]) # Label each point on the graph for (s, x, y) in merged[['seqid', 'length', 'gene_count']].sort_values('length').values: ax.text(x, y - 100, str(s))

As seen in image above, even though it is a positive correlation overall, it does not hold for all chromosomes. In particular, for Chromosome 17, 16, 15, 14, 13, the correlation is actually negative, meaning the number of genes on the chromosome decreases as the chromosome size increases.

Findings and Future Research

That ends our tutorial on the manipulation of an annotation file for human genome in GFF3 format with the SciPy stack. The tools we've mainly used include IPython, Pandas, and matplotlib. During the tutorial, not only have we learned some of the most common and useful operations in Pandas, we also answered some very interesting questions about our genome. In summary:

จีโนมมนุษย์ประมาณ 0.37% ยังไม่สมบูรณ์แม้ว่าร่างแรกจะออกมาเมื่อ 15 ปีที่แล้วก็ตาม
มียีนประมาณ 42,000 ยีนในจีโนมมนุษย์โดยอิงจากไฟล์ GFF3 เฉพาะที่เราใช้
ความยาวของยีนอาจมีตั้งแต่ไม่กี่โหลไปจนถึงมากกว่าสองล้านเบส
ยีนไม่กระจายไปตามโครโมโซม โดยรวมแล้ว ยิ่งโครโมโซมมีขนาดใหญ่เท่าใด ยีนก็ยิ่งโฮสต์มากขึ้นเท่านั้น แต่สำหรับชุดย่อยของโครโมโซม ความสัมพันธ์อาจเป็นลบได้

ไฟล์ GFF3 มีข้อมูลคำอธิบายประกอบมากมาย และเราเพิ่งขีดข่วนพื้นผิว หากคุณสนใจที่จะสำรวจเพิ่มเติม ต่อไปนี้คือคำถามสองสามข้อที่คุณสามารถเล่นได้:

ยีนปกติมีกี่สำเนา? เปอร์เซ็นต์ของยีนที่มีมากกว่า 1 การถอดเสียง?
การถอดเสียงโดยทั่วไปมีกี่ไอโซฟอร์ม?
โดยทั่วไปแล้วการถอดเสียงจะมี exons, CDS และ UTR จำนวนเท่าใด มีขนาดใดบ้าง?
เป็นไปได้ไหมที่จะจัดหมวดหมู่ยีนตามหน้าที่ของพวกมันตามที่อธิบายไว้ในคอลัมน์คำอธิบาย?