مقدمة شاملة عن الجينوم الخاص بك باستخدام SciPy Stack

المعلوماتية الحيوية هي مجال متعدد التخصصات يطور أساليب وأدوات برمجية لتحليل وفهم البيانات البيولوجية.

وبعبارة أكثر بساطة ، يمكنك ببساطة التفكير في الأمر على أنه علم بيانات للبيولوجيا.

من بين الأنواع العديدة من البيانات البيولوجية ، تعد بيانات الجينوم واحدة من أكثر البيانات التي تم تحليلها على نطاق واسع. خاصة مع التقدم السريع لتقنيات تسلسل الحمض النووي من الجيل التالي (NGS) ، كان حجم بيانات الجينوم ينمو بشكل كبير. وفقًا لستيفنس وزاكاري دي وآخرون ، يتم الحصول على بيانات الجينوم على مقياس إكسابايت سنويًا.

في هذا المنشور ، أعرض مثالاً على تحليل ملف GFF3 للجينوم البشري باستخدام SciPy Stack. الإصدار 3 من تنسيق الميزة العامة (GFF3) هو تنسيق الملف النصي القياسي الحالي لتخزين الميزات الجينية. على وجه الخصوص ، ستتعلم في هذا المنشور كيفية استخدام مكدس SciPy للإجابة على الأسئلة التالية حول الجينوم البشري:

1. كم من الجينوم غير مكتمل؟
2. كم عدد الجينات الموجودة في الجينوم؟
3. ما هي مدة الجين النموذجي؟
4. كيف يبدو توزيع الجينات بين الكروموسومات؟

يمكن تنزيل أحدث ملف GFF3 للجينوم البشري من هنا. يوفر ملف README الذي يأتي في هذا الدليل وصفًا موجزًا ​​لتنسيق البيانات هذا ، ويتم العثور هنا على مواصفات أكثر شمولاً.

سنستخدم Pandas ، وهو مكون رئيسي في مكدس SciPy يوفر هياكل بيانات سريعة ومرنة ومعبرة ، لمعالجة ملف 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 بالتنفيذ في وضع الدُفعات. بعد استخدام الأوامر أعلاه لتثبيت Miniconda بنجاح ، ابدأ بيئة افتراضية جديدة لعلم الجينوم ، ثم قم بتثبيت حزمة SciPy.

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

لاحظ أننا حددنا فقط الحزم الثلاثة التي سنستخدمها في هذا المنشور. إذا كنت تريد جميع الحزم المدرجة في مكدس SciPy ، فما عليك سوى إلحاقها بنهاية الأمر conda create .

إذا لم تكن متأكدًا من الاسم الدقيق للحزمة ، فجرّب conda search . لنقم بتنشيط البيئة الافتراضية ونبدأ IPython.

 source activate venv/ ipython

يعد IPython بديلاً أقوى بشكل ملحوظ لواجهة مترجم Python الافتراضية ، لذلك كل ما كنت تفعله في مترجم Python الافتراضي يمكن أيضًا القيام به في IPython. أوصي بشدة كل مبرمج بايثون ، لم يستخدم IPython حتى الآن ، لتجربته.

قم بتنزيل ملف التعليقات التوضيحية

بعد اكتمال الإعداد الآن ، لنقم بتنزيل ملف شرح الجينوم البشري بتنسيق GFF3.

يبلغ حجمه حوالي 37 ميجابايت ، وهو ملف صغير جدًا مقارنة بمحتوى المعلومات الخاص بجينوم بشري ، والذي يبلغ حجمه حوالي 3 جيجابايت في نص عادي. ذلك لأن ملف GFF3 يحتوي فقط على تعليق توضيحي للتسلسلات ، بينما يتم تخزين بيانات التسلسل عادةً في تنسيق ملف آخر يسمى FASTA. إذا كنت مهتمًا ، يمكنك تنزيل FASTA هنا ، لكننا لن نستخدم بيانات التسلسل في هذا البرنامج التعليمي.

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

البادئة ! يخبر IPython أن هذا أمر shell بدلاً من أمر Python. ومع ذلك ، يمكن لـ IPython أيضًا معالجة بعض أوامر shell المستخدمة بشكل متكرر مثل ls و pwd و rm و mkdir و rmdir حتى بدون بادئة ! .

بإلقاء نظرة على رأس ملف GFF ، سترى العديد من أسطر البيانات الوصفية / العملية / التوجيهات تبدأ بـ ## أو #! .

وفقًا لـ 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 ، الذي ينطبق عليه ملف التعليقات التوضيحية هذا.

اتحاد مرجع الجينوم (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 ، المصدر ، النوع ، البداية ، النهاية ، النتيجة ، الخيط ، المرحلة ، السمات . بعضها سهل الفهم. خذ السطر الأول كمثال:

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

هذا هو التعليق التوضيحي للكروموسوم الأول مع seqid من 1 ، والذي يبدأ من القاعدة الأولى إلى القاعدة 24895.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 ، ويستخدم لتمثيل البيانات ثنائية الأبعاد.

تمتلك 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 العديد من الخيارات التي يمكن تخصيصها.

بعد ذلك ، دعنا نحصل على بعض المعلومات الأساسية حول إطار البيانات هذا باستخدام طريقة .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 ميجا بايت مخزنة في الذاكرة. تبين أن هذا يكون أكثر إحكاما من الملف غير المضغوط ، والذي سيكون حوالي 402 ميغابايت. يتم عرض تحقق سريع أدناه.

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

من وجهة نظر عالية المستوى ، قمنا بتحميل ملف GFF3 بأكمله في كائن DataFrame في Python ، وستعتمد جميع تحليلاتنا التالية على هذا الكائن الفردي.

الآن ، دعنا نرى ما هو العمود الأول seqid.

 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) بالإضافة إلى 169 سلسلة أخرى.

إن seqids الذي يبدأ بـ KI و GL عبارة عن تسلسلات DNA - أو سقالات - في الجينوم الذي لم يتم تجميعه بنجاح في الكروموسومات.

بالنسبة لأولئك الذين ليسوا على دراية بعلم الجينوم ، هذا مهم.

على الرغم من ظهور أول مسودة للجينوم البشري منذ أكثر من 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 و havana و ensembl و 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...

التصفية سهلة في الباندا.

إذا قمت بفحص القيمة التي تم تقييمها من التعبير df.source == 'GRCh38' ، فهي سلسلة من قيم True و False لكل إدخال بنفس فهرس df . تمريره إلى df[] سيؤدي فقط إلى إرجاع تلك الإدخالات حيث تكون القيم المقابلة لها صحيحة.

يوجد 194 مفتاحًا في df[] من أجلها df.source == 'GRCh38' .

كما رأينا سابقًا ، هناك أيضًا 194 قيمة فريدة في عمود seqid ، مما يعني أن كل إدخال في gdf يتوافق مع seqid معين.

ثم نختار 10 إدخالات بشكل عشوائي باستخدام طريقة sample لإلقاء نظرة فاحصة.

يمكنك أن ترى أن التسلسلات غير المجمعة هي من النوع الفائق بينما الآخرون من الكروموسوم. لحساب جزء الجينوم غير المكتمل ، نحتاج أولاً إلى معرفة طول الجينوم بأكمله ، والذي يمثل مجموع أطوال جميع السيكولات.

 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 و 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 وتسلسل التشفير (CDS) والمنطقة غير المترجمة (UTR).

هذه عناصر جينية فرعية ، لكننا نبحث بشكل أساسي عن عدد الجينات. كما هو موضح ، هناك 42470 ، لكننا نريد معرفة المزيد.

على وجه التحديد ، ما هي أسمائهم ، وماذا يفعلون؟ للإجابة على هذه الأسئلة ، نحتاج إلى إلقاء نظرة فاحصة على المعلومات الموجودة في عمود السمات.

 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 ، وهي طريقة يتم استخدامها عند الحاجة إلى تطبيق دالة على كل إدخال لإطار بيانات أو سلسلة.

في هذه الحالة بالذات ، نريد استخراج اسم الجين لكل إدخال في 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)

يعتبر التعبير المعتاد لـ RE_GENE_ID أكثر تحديدًا لأننا نعلم أن كل gene_id يجب أن يبدأ بـ ENSG ، حيث يعني ENS ensembl و G يعني الجين.

بالنسبة للإدخالات التي لا تحتوي على أي وصف ، سنعود سلسلة فارغة. بعد استخراج كل شيء ، لن نستخدم عمود السمات بعد الآن ، لذا دعنا نتخلص منه للحفاظ على الأشياء لطيفة ونظيفة باستخدام الطريقة .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,)

والمثير للدهشة أن عدد أسماء الجينات أقل من عدد معرفات الجينات ، مما يشير إلى أن بعض اسم الجين يجب أن يتوافق مع معرفات جينية متعددة. دعونا نكتشف ما هم.

 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() بعض الإحصائيات البسيطة بناءً على قيم الطول:

• متوسط ​​طول الجين حوالي 36000 قاعدة

• يبلغ متوسط ​​طول الجين حوالي 5200 قاعدة

• يبلغ الحد الأدنى والحد الأقصى لأطوال الجينات حوالي ثمانية و 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 للجمع بين اثنين من إطارات البيانات. بشكل حدسي ، من المحتمل أن تستضيف الكروموسومات الأطول المزيد من الجينات. دعونا نرى ما إذا كان هذا صحيحًا.

 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 متغير chrs من القسم السابق ، واستخدمناه لتصفية التسلسلات غير المجمعة. بناءً على المخرجات ، يحتوي أكبر كروموسوم 1 بالفعل على معظم الجينات. بينما يحتوي الكروموسوم Y على أقل عدد من الجينات ، فإنه ليس أصغر كروموسوم.

لاحظ أنه يبدو أنه لا توجد جينات في الميتوكوندريا (MT) ، وهذا ليس صحيحًا.

مزيد من التصفية على أول DataFrame df الذي تم إرجاعه بواسطة pd.read_csv يُظهر أن جميع جينات MT هي من المصدر insdc (والتي تمت تصفيتها من قبل عند إنشاء edf حيث نظرنا فقط في مصادر havana أو ensembl أو 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:

1. لا يزال حوالي 0.37٪ من الجينوم البشري غير مكتمل على الرغم من ظهور المسودة الأولى منذ أكثر من 15 عامًا.
2. يوجد حوالي 42000 جينة في الجينوم البشري بناءً على ملف GFF3 المحدد الذي استخدمناه.
3. يمكن أن يتراوح طول الجين من بضع عشرات إلى أكثر من مليوني قاعدة.
4. لا يتم توزيع الجينات بالتساوي بين الكروموسومات. بشكل عام ، كلما زاد حجم الكروموسوم ، زاد عدد الجينات التي يستضيفها ، ولكن بالنسبة لمجموعة فرعية من الكروموسومات ، يمكن أن يكون الارتباط سالبًا.

ملف GFF3 غني جدًا بمعلومات التعليقات التوضيحية ، وقد خدشنا السطح للتو. إذا كنت مهتمًا بمزيد من الاستكشاف ، فإليك بعض الأسئلة التي يمكنك اللعب بها:

1. كم عدد النسخ التي يمتلكها الجين عادة؟ ما هي النسبة المئوية للجينات التي لديها أكثر من نسخة واحدة؟
2. كم عدد الأشكال الإسوية التي يحتويها النص عادةً؟
3. كم عدد exons و CDS و UTRs التي يمتلكها النص عادةً؟ ما هي الأحجام؟
4. هل يمكن تصنيف الجينات بناءً على وظيفتها كما هو موضح في عمود الوصف؟