srcmini本文概述
生物信息学是一个跨学科领域, 致力于开发用于分析和理解生物数据的方法和软件工具。
简而言之, 你可以简单地将其视为生物学的数据科学。
在许多类型的生物学数据中, 基因组数据是分析最广泛的数据之一。尤其是随着下一代DNA测序(NGS)技术的飞速发展, 基因组学数据的数量呈指数增长。根据Zachary D等人的Stephens的说法, 基因组数据的采集量为每年EB级。
SciPy收集了许多用于科学计算的Python模块, 非常适合许多生物信息学需求。
鸣叫
在本文中, 我演示了一个使用SciPy Stack分析人类基因组的GFF3文件的示例。通用特征格式版本3(GFF3)是用于存储基因组特征的当前标准文本文件格式。特别是, 在本文中, 你将学习如何使用SciPy堆栈回答有关人类基因组的以下问题：
- 有多少基因组不完整？
- 基因组中有多少个基因？
- 一个典型的基因多长时间？
- 染色体之间的基因分布是什么样的？
可以从此处下载有关人类基因组的最新GFF3文件。此目录中的README文件提供了此数据格式的简要说明, 并在此处找到了更详尽的规范。
我们将使用SciPy堆栈的主要组件Pandas提供快速, 灵活和富于表现力的数据结构, 以操纵和理解GFF3文件。
设定
首先, 让我们在安装了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搜索。让我们激活虚拟环境并启动IPython。
source activate venv/
ipython
IPython是默认Python解释器接口的强大得多的替代品, 因此你在IPython中使用默认Python解释器所做的任何事情也可以完成。我强烈建议所有尚未使用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文件的开头, 你会看到许多以##或＃！开头的元数据/编译指示/指令行。
根据自述文件, ##表示元数据是稳定的, 而＃！表示元数据是稳定的。表示这是实验性的。
稍后, 你还将看到###, 这是基于规范的另一个具有更微妙含义的指令。
人们可读的注释应该在单个＃之后。为简单起见, 我们将以＃开头的所有行都视为注释, 并在分析过程中将其忽略。
##gff-version 3
##sequence-region 1 1 248956422
##sequence-region 10 1 133797422
##sequence-region 11 1 135086622
##sequence-region 12 1 133275309
...
##sequence-region MT 1 16569
##sequence-region X 1 156040895
##sequence-region Y 2781480 56887902
#!genome-build GRCh38.p7
#!genome-version GRCh38
#!genome-date 2013-12
#!genome-build-accession NCBI:GCA_000001405.22
#!genebuild-last-updated 2016-06
第一行表示此文件中使用的GFF格式的版本为3。
接下来是所有序列区域的摘要。正如我们稍后将看到的, 此类信息也可以在文件的正文部分中找到。
以＃开头的行！显示有关此注释文件适用的基因组特定构建GRCh38.p7的信息。
Genome Reference Consortium(GCR)是一个国际性的联盟, 负责监督几个参考基因组组件的更新和改进, 包括用于人类, 小鼠, 斑马鱼和鸡的那些。
扫描此文件, 这里是前几行注释行。
1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2, chr1, NC_000001.11
###
1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg
1 . biological_region 10650 10657 0.999 + . logic_name=eponine
1 . biological_region 10655 10657 0.999 - . logic_name=eponine
1 . biological_region 10678 10687 0.999 + . logic_name=eponine
1 . biological_region 10681 10688 0.999 - . logic_name=eponine
...
列是顺序, 来源, 类型, 开始, 结束, 得分, 链, 阶段, 属性。其中一些很容易理解。以第一行为例：
1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2, chr1, NC_000001.11
这是第一个染色体的注释, 后继序列为1, 从第一个碱基开始到第24, 895, 622个碱基。
换句话说, 第一条染色体长约2500万个碱基。
我们的分析将不需要来自三列的值为的信息。 (即得分, 分流和阶段), 因此我们暂时可以将其忽略。
最后一个属性列表示染色体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数据。
熊猫还具有分别用于1D和3D数据的系列和面板数据结构。请参阅文档以获取有关熊猫数据结构的介绍。
让我们看一下.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设置为不省略长字符串。此外, 熊猫还有许多可以定制的选项。
接下来, 让我们使用.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列。
对于每一列, 它还显示其数据类型。
开头和结尾是int64类型, 整数表示基因组中的位置。
其他列均为对象类型, 这可能意味着它们的值由整数, 浮点数和字符串组成。
存储在内存中的所有信息的大小约为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个独特的序列, 其中包括1至22号染色体, X, Y和线粒体(MT)DNA, 以及其他169个序列。
以KI和GL开头的片段是基因组中尚未成功组装成染色体的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方法使用示例, 对于快速计数分类变量非常有用。
从结果中, 我们可以看到此列有七个可能的值, 并且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 []只会返回其对应值为True的那些条目。
df []中有194个键, 其中df.source ==’GRCh38’。
如我们先前所见, seqid列中还有194个唯一值, 这意味着gdf中的每个条目都对应于一个特定的seqid。
然后, 我们使用样本方法随机选择10个条目以进行仔细研究。
你会看到未组装的序列属于supercontig类型, 而其他序列属于染色体。要计算不完整的基因组比例, 我们首先需要知道整个基因组的长度, 即所有序列的长度之和。
In [90]: gdf = gdf.copy()
In [91]: gdf['length'] = gdf.end - gdf.start + 1
In [93]: gdf.head()
Out[93]:
seqid source type start end score strand phase attributes length
0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2, chr1, NC_00000... 248956421
235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2, chr10, NC_000... 133797421
328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2, chr11, NC_000... 135086621
483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2, chr12, NC_000... 133275308
634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2, chr13, NC_000... 114364327
In [97]: gdf.length.sum()
Out[97]: 3096629532
In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT']
In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum()
Out[101]: 0.0037021917421198327
首先在上面的代码段中, 我们使用.copy()复制了gdf。否则, 原始gdf只是df的一部分, 直接对其进行修改将导致SettingWithCopyWarning(更多信息, 请参见此处)。
然后, 我们计算每个条目的长度, 并将其添加回gdf, 作为名为” length”的新列。结果总长度约为31亿, 未组装序列的比例约为0.37％。
这是最后两个命令中切片的工作方式。
首先, 我们创建一个字符串列表, 其中包含装配良好的序列的所有序列, 它们都是染色体和线粒体。然后, 我们使用isin方法过滤scrid在chrs列表中的所有条目。
在索引的开头添加了一个减号(-)以反转选择, 因为我们实际上希望列表中没有的所有内容(即我们希望以KI和GL开头的未汇编的内容)…
注意：由于组装后的序列和未组装的序列由type列区分, 因此可以替代地如下重写最后一行以获得相同的结果。
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表示基因。
对于没有任何说明的条目, 我们将返回一个空字符串。提取所有内容后, 我们将不再使用attribute列, 因此, 请使用.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...
在上面的调用中, 属性指示我们要删除的特定列。
axis = 1表示我们要删除列而不是行(默认情况下轴= 0)。
inplace = True表示放置操作是在DataFrame本身上进行的, 而不是返回放置了指定列的新副本。
快速浏览.head可以看到attribute列确实消失了, 并且添加了三个新列：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个碱基
-
最小和最大基因长度分别约为八和230万个碱基。
由于均值比中位数大得多, 这意味着长度分布向右偏斜。为了更具体的外观, 让我们绘制分布图。
import matplotlib as plt
ndf.length.plot(kind='hist', bins=50, logy=True)
plt.show()
熊猫为matplotlib提供了一个简单的界面, 使使用DataFrames或series进行绘图非常方便。
在这种情况下, 它表示我们想要一个具有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相关的蛋白样2的简称。根据其Wikipedia页面,
该基因几乎涵盖了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返回的第一个DataFrame 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...
此示例还显示了如何在过滤时使用＆;组合两个条件。 “或”的逻辑运算符为|。
请注意, 每个条件都需要用括号括起来, Pandas中语法的这一部分与Python不同, 后者应为文字和和或。
接下来, 让我们借鉴上一节中的gdf DataFrame作为每个染色体长度的来源：
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
为了清楚起见, 删除了与分析无关的列。
是的, .drop还可以采用一列列表并将其完全删除。
注意, 带有MT序列的行仍然存在；我们稍后会再讲。我们将执行的下一个操作是基于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
由于chr_gene_counts仍然是Series对象, 不支持合并操作, 因此我们需要首先使用.to_frame将其转换为DataFrame对象。
.reset_index()将原始索引(即seqid)转换为新列, 并将当前索引重置为基于0的增量数字。
cdf.head(2)的输出显示了外观。接下来, 我们使用.merge方法在seqid列(on =’seqid’)上合并两个DataFrame。
合并gdf和cdf之后, MT条目仍然丢失。这是因为, 默认情况下, .merge操作内部联接, 而通过调整how参数可以使用左联接, 右联接或外联接。
请参阅文档以获取更多详细信息。
稍后, 你可能会发现还有一个相关的.join方法。 .merge和.join相似, 但具有不同的API。
据官方文件说
相关的DataFrame.join方法在内部使用index-on和index-on-column联接合并, 但是默认情况下在索引上联接, 而不是尝试在公共列上联接(合并的默认行为)。如果你要加入索引, 则可能希望使用DataFrame.join节省一些输入。
基本上, .merge更通用, 由.join使用。
最后, 我们准备计算染色体长度与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
默认情况下, .corr计算数据帧中所有列对之间的Pearson相关性。
但是在这种情况下, 我们只有一对列, 相关性为正, 为0.73。
换句话说, 染色体越大, 拥有更多基因的可能性就越大。在按长度对值对进行排序之后, 我们还要绘制两列：
ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count', style='o-')
# 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))
如上图所示, 即使总体上呈正相关, 但并非所有染色体都适用。特别是, 对于17号, 16号, 15号, 14号, 13号染色体, 相关性实际上是负的, 这意味着染色体上的基因数量随着染色体大小的增加而减少。
研究结果和未来研究
到此, 我们结束了使用SciPy堆栈操作GFF3格式的人类基因组注释文件的教程。我们主要使用的工具包括IPython, Pandas和matplotlib。在本教程中, 我们不仅学习了熊猫中一些最普通和有用的操作, 还回答了一些有关基因组的非常有趣的问题。综上所述：
- 即使15年前提出了第一稿, 仍约有0.37％的人类基因组不完整。
- 根据我们使用的特定GFF3文件, 人类基因组中约有42, 000个基因。
- 一个基因的长度范围可以从几十个到超过两百万个碱基。
- 基因在染色体之间分布不均。总体而言, 染色体越大, 其宿主的基因越多, 但是对于染色体的一个子集, 相关性可能为负。
GFF3文件中的注释信息非常丰富, 而我们才刚刚开始。如果你对进一步的探索感兴趣, 可以考虑以下几个问题：
- 一个基因通常有多少个转录本？超过1个转录本的基因百分比是多少？
- 成绩单通常有几个同工型？
- 笔录通常有多少个外显子, CDS和UTR？它们是什么尺寸？
- 是否可以根据描述栏中描述的功能对基因进行分类？
评论前必须登录！
注册