使用 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.

有关详细信息，请参阅文档。

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？它们的尺寸是多少？
是否可以根据描述栏中所述的功能对基因进行分类？