使用 SciPy 堆棧全面介紹您的基因組

已發表: 2022-03-11

生物信息學是一個跨學科領域，它開髮用於分析和理解生物數據的方法和軟件工具。

更簡單地說，您可以簡單地將其視為生物學的數據科學。

在眾多類型的生物數據中，基因組數據是分析最廣泛的數據之一。尤其是隨著下一代 DNA 測序 (NGS) 技術的快速發展，基因組學數據量呈指數級增長。根據 Stephens、Zachary D 等人的說法，基因組學數據採集的規模為每年艾字節。

SciPy 收集了許多用於科學計算的 Python 模塊，非常適合許多生物信息學需求。

鳴叫

在這篇文章中，我演示了一個使用 SciPy Stack 分析人類基因組的 GFF3 文件的示例。通用特徵格式版本 3 (GFF3) 是當前用於存儲基因組特徵的標准文本文件格式。特別是，在這篇文章中，您將學習如何使用 SciPy 堆棧來回答以下有關人類基因組的問題：

有多少基因組是不完整的？
基因組中有多少個基因？
一個典型的基因有多長？
染色體之間的基因分佈是什麼樣的？

可以從這裡下載人類基因組的最新 GFF3 文件。此目錄中的 README 文件提供了此數據格式的簡要說明，更詳盡的規範可在此處找到。

我們將使用 Pandas 來操作和理解 GFF3 文件，它是 SciPy 堆棧的主要組件，提供快速、靈活和富有表現力的數據結構。

設置

首先，讓我們設置一個安裝了 SciPy 堆棧的虛擬環境。如果從源代碼手動構建此過程可能會很耗時，因為堆棧涉及許多包——其中一些依賴於外部 FORTRAN 或 C 代碼。在這裡，我推薦使用 Miniconda，它使設置變得非常容易。

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

bash 行上的-b標誌告訴它以批處理模式執行。使用上述命令成功安裝 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 中完成。我強烈建議所有尚未使用 IPython 的 Python 程序員嘗試一下。

下載註釋文件

現在我們的設置完成了，讓我們下載 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 這是一個 shell 命令而不是 Python 命令。然而，即使沒有前綴，IPython 也可以處理一些常用的 shell 命令，如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 的信息。

基因組參考聯盟 (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 、 attributes 。其中一些非常容易理解。以第一行為例：

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

這是第一個seqid為1的染色體的註釋，從第一個鹼基開始到第24,895,622個鹼基。

換句話說，第一條染色體的長度約為 2500 萬個鹼基。

我們的分析不需要來自值為的三列的信息. （即分數、鍊和相位），所以我們現在可以簡單地忽略它們。

最後一個屬性列說，Chromosome 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)

上面的最後一行使用pandas.read_csv方法加載整個 GFF3 文件。

由於它不是標準的 CSV 文件，我們需要稍微自定義調用。

首先，我們用header=None通知 Pandas GFF3 中的標題信息不可用，然後我們用names=col_names指定每一列的確切名稱。

如果未指定names參數，Pandas 將使用遞增數字作為每列的名稱。

sep='\t'告訴 Pandas 列是製表符分隔的，而不是逗號分隔的。 sep=的值實際上可以是正則表達式 (regex)。如果手頭的文件對每列使用不同的分隔符，這將變得很方便（哦，是的，這種情況發生了）。 comment='#'表示以#開頭的行被視為註釋，將被忽略。

compression='gzip'告訴 Pandas 輸入文件是 gzip 壓縮的。

此外， pandas.read_csv具有豐富的參數集，可以讀取不同種類的類似 CSV 的文件格式。

返回值的類型是DataFrame ，它是 Pandas 中最重要的數據結構，用於表示 2D 數據。

Pandas 還具有分別用於 1D 和 3D 數據的Series和Panel數據結構。有關 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

輸出以表格格式很好地格式化，屬性列中較長的字符串部分替換為... 。

您可以使用pd.set_option('display.max_colwidth', -1)將 Pandas 設置為不省略長字符串。此外，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 行註釋，每行有 9 列。

對於每一列，它還顯示了它們的數據類型。

start 和 end 是 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 文件加載到 Python 中的一個 DataFrame 對像中，我們接下來的所有分析都將基於這個單個對象。

現在，讓我們看看第一列 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 個獨特的 seqid，其中包括染色體 1 到 22、X、Y 和線粒體 (MT) DNA 以及 169 個其他 seqid。

以 KI 和 GL 開頭的 seqid 是基因組中尚未成功組裝到染色體中的 DNA 序列或支架。

對於那些不熟悉基因組學的人來說，這很重要。

儘管第一個人類基因組草圖在 15 多年前問世，但目前的人類基因組仍然不完整。組裝這些序列的困難主要是由於基因組中復雜的重複區域。

接下來，我們來看看源列。

自述文件說來源是一個自由文本限定符，旨在描述生成此功能的算法或操作過程。

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

這是使用value_counts方法的示例，該方法對於快速計算分類變量非常有用。

從結果中，我們可以看到該列有 7 個可能的值，並且 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...

在 Pandas 中過濾很容易。

如果您檢查從表達式df.source == 'GRCh38'評估的值，它是具有與df相同索引的每個條目的一系列True和False值。將其傳遞給df[]只會返回對應值為 True 的條目。

df[]中有 194 個鍵，其中df.source == 'GRCh38' 。

正如我們之前看到的，在seqid列中也有 194 個唯一值，這意味著 gdf 中的每個條目對應於一個特定的gdf 。

然後我們用sample方法隨機選擇10個條目仔細觀察。

您可以看到未組裝的序列屬於超重疊群類型，而其他序列屬於染色體。要計算不完整基因組的比例，我們首先需要知道整個基因組的長度，即所有 seqid 長度的總和。

 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

在上面的代碼片段中，我們首先使用.copy()製作了gdf的副本。否則，原始gdf只是df的一部分，直接修改它會導致SettingWithCopyWarning （有關詳細信息，請參閱此處）。

然後我們計算每個條目的長度並將其作為名為“length”的新列添加回gdf 。結果總長度約為 31 億，未組裝序列的比例約為 0.37%。

以下是最後兩個命令中切片的工作方式。

首先，我們創建一個字符串列表，其中包含所有組裝良好的序列的 seqid，它們都是染色體和線粒體。然後我們使用isin方法過濾所有 seqid 在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方法再次用於過濾。

然後，快速值計數顯示大多數條目是外顯子、編碼序列 (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)

它們被格式化為以分號分隔的標籤-值對列表。我們最感興趣的信息是基因名稱、基因ID和描述，我們會用正則表達式（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)

首先，我們提取基因名稱。

在正則表達式中Name=(?P<gene_name>.+?); , +? 使用而不是+因為我們希望它是非貪婪的，並讓搜索在第一個分號處停止；否則，結果將匹配到最後一個分號。

此外，正則表達式首先使用re.compile進行編譯，而不是像在re.search中那樣直接使用以獲得更好的性能，因為我們會將其應用於數千個屬性字符串。

extract_gene_name用作要在pd.apply中使用的輔助函數，這是在需要將函數應用於數據幀或序列的每個條目時使用的方法。

在這種特殊情況下，我們要為ndf.attributes中的每個條目提取基因名稱，並將名稱添加回gene_name ndf新列中。

以類似的方式提取基因 ID 和描述。

 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,)

令人驚訝的是，基因名稱的數量少於基因ID的數量，這表明某些gene_name必須對應多個基因ID。讓我們找出它們是什麼。

 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]反轉順序。

結果，一個基因名稱最多可以與七個基因 ID 共享。

 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 基因表明它們確實都不同。

雖然它們具有相同的名稱，但它們位於基因組的不同位置。

然而，它們的描述似乎對區分它們沒有多大幫助。

這裡的重點是要知道基因名稱對於所有基因 ID 來說並不是唯一的，並且大約 0.15% 的名稱由多個基因共享。

一個典型的基因有多長？

與我們試圖了解基因組不完整性時所做的類似，我們可以輕鬆地在ndf中添加一個length列：

 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 個鹼基
最小和最大基因長度分別約為 8 和 230 萬個鹼基。

因為平均值遠大於中位數，這意味著長度分佈向右傾斜。為了更具體地了解，讓我們繪製分佈圖。

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

Pandas 為 matplotlib 提供了一個簡單的接口，以便使用 DataFrames 或系列進行繪圖非常方便。

在這種情況下，它表示我們想要一個具有 50 個 bin 的直方圖（ kind='hist' ），並讓 y 軸在對數刻度上（ logy=True ）。

從直方圖中，我們可以看到大部分基因都在第一個 bin 內。然而，一些基因長度可能超過 200 萬個鹼基。讓我們找出它們是什麼：

 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 associated protein-like 2 的縮寫。根據其維基百科頁面，

該基因幾乎佔 7 號染色體的 1.6%，是人類基因組中最大的基因之一。

確實！我們只是自己驗證了這一點。相反，最小的基因呢？事實證明，它們可以短至八個鹼基。

 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

兩種極端情況的長度相差五個數量級（230 萬對 8），這是巨大的，可以表明生命的多樣性水平。

單個基因可以通過稱為選擇性剪接的過程翻譯成許多不同的蛋白質，這是我們尚未探索過的。此類信息也在 GFF3 文件中，但不在本文的範圍內。

染色體間的基因分佈

我想討論的最後一件事是染色體之間的基因分佈，這也可以作為介紹合併兩個 DataFrame 的.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) 中似乎沒有基因，這是不正確的。

對 pd.read_csv 返回的第一個pd.read_csv df進行更多過濾表明所有 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. 總之：

大約 0.37% 的人類基因組仍然不完整，儘管第一稿在 15 年前問世。
根據我們使用的這個特殊的 GFF3 文件，人類基因組中有大約 42,000 個基因。
一個基因的長度可以從幾十個鹼基到超過兩百萬個鹼基不等。
基因在染色體中分佈不均。總的來說，染色體越大，它所承載的基因越多，但對於染色體的一個子集，相關性可能是負的。

GFF3 文件的註解信息非常豐富，而我們只是觸及了皮毛。如果你有興趣進一步探索，這裡有幾個問題你可以玩：

一個基因通常有多少個轉錄本？有多少百分比的基因有超過 1 個轉錄本？
轉錄本通常有多少種亞型？
轉錄本通常有多少個外顯子、CDS 和 UTR？它們的尺寸是多少？
是否可以根據描述欄中所述的功能對基因進行分類？