本文概述
- 设定
- 下载注释文件
- 研究结果和未来研究
简而言之, 你可以简单地将其视为生物学的数据科学。
在许多类型的生物学数据中, 基因组数据是分析最广泛的数据之一。尤其是随着下一代DNA测序(NGS)技术的飞速发展, 基因组学数据的数量呈指数增长。根据Zachary D等人的Stephens的说法, 基因组数据的采集量为每年EB级。
文章图片
SciPy收集了许多用于科学计算的Python模块, 非常适合许多生物信息学需求。
鸣叫
在本文中, 我演示了一个使用SciPy Stack分析人类基因组的GFF3文件的示例。通用特征格式版本3(GFF3)是用于存储基因组特征的当前标准文本文件格式。特别是, 在本文中, 你将学习如何使用SciPy堆栈回答有关人类基因组的以下问题：
- 有多少基因组不完整？
- 基因组中有多少个基因？
- 一个典型的基因多长时间？
- 染色体之间的基因分布是什么样的？
我们将使用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-version3
##sequence-region1 1 248956422
##sequence-region10 1 133797422
##sequence-region11 1 135086622
##sequence-region12 1 133275309
...
##sequence-regionMT 1 16569
##sequence-regionX 1 156040895
##sequence-regionY 2781480 56887902
#!genome-buildGRCh38.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)是一个国际性的联盟, 负责监督几个参考基因组组件的更新和改进, 包括用于人类, 小鼠, 斑马鱼和鸡的那些。
扫描此文件, 这里是前几行注释行。
1GRCh38chromosome1248956422...ID=chromosome:1;
Alias=CM000663.2, chr1, NC_000001.11
###
1.biological_region10469112401.3e+03 ..external_name=oe %3D 0.79;
logic_name=cpg
1.biological_region10650106570.999+.logic_name=eponine
1.biological_region10655106570.999-.logic_name=eponine
1.biological_region10678106870.999+.logic_name=eponine
1.biological_region10681106880.999-.logic_name=eponine
...
列是顺序, 来源, 类型, 开始, 结束, 得分, 链, 阶段, 属性。其中一些很容易理解。以第一行为例：
1GRCh38chromosome1248956422...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 pdIn [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]:
seqidsourcetypestartendscore strand phaseattributes
01GRCh38chromosome1248956422...ID=chromosome:1;
Alias=CM000663.2, chr1, NC_00000...
11.biological_region10469112401.3e+03..external_name=oe %3D 0.79;
logic_name=cpg
21.biological_region10650106570.999+.logic_name=eponine
31.biological_region10655106570.999-.logic_name=eponine
41.biological_region10678106870.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):
seqidobject
sourceobject
typeobject
startint64
endint64
scoreobject
strandobject
phaseobject
attributesobject
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]:
havana1441093
ensembl_havana745065
ensembl228212
.182510
mirbase4701
GRCh38194
insdc74
这是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]:
seqidsourcetypestartend score strand phaseattributes
2511585KI270708.1GRCh38supercontig1127682...ID=supercontig:KI270708.1;
Alias=chr1_KI270708v...
2510840GL000208.1GRCh38supercontig192689...ID=supercontig:GL000208.1;
Alias=chr5_GL000208v...
99081017GRCh38chromosome183257441...ID=chromosome:17;
Alias=CM000679.2, chr17, NC_000...
2511481KI270373.1GRCh38supercontig11451...ID=supercontig:KI270373.1;
Alias=chrUn_KI270373...
2511490KI270384.1GRCh38supercontig11658...ID=supercontig:KI270384.1;
Alias=chrUn_KI270384...
20801486GRCh38chromosome1170805979...ID=chromosome:6;
Alias=CM000668.2, chr6, NC_00000...
2511504KI270412.1GRCh38supercontig11179...ID=supercontig:KI270412.1;
Alias=chrUn_KI270412...
120156119GRCh38chromosome158617616...ID=chromosome:19;
Alias=CM000681.2, chr19, NC_000...
2511474KI270340.1GRCh38supercontig11428...ID=supercontig:KI270340.1;
Alias=chrUn_KI270340...
2594560YGRCh38chromosome278148056887902...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]:
seqidsourcetypestartend score strand phaseattributeslength
01GRCh38chromosome1248956422...ID=chromosome:1;
Alias=CM000663.2, chr1, NC_00000...248956421
23506810GRCh38chromosome1133797422...ID=chromosome:10;
Alias=CM000672.2, chr10, NC_000...133797421
32893811GRCh38chromosome1135086622...ID=chromosome:11;
Alias=CM000673.2, chr11, NC_000...135086621
48337012GRCh38chromosome1133275309...ID=chromosome:12;
Alias=CM000674.2, chr12, NC_000...133275308
63448613GRCh38chromosome1114364328...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]:
seqidsourcetypestartend score strand phaseattributes
91599616havanaCDS2746354127463592.-2ID=CDS:ENSP00000457449;
Parent=transcript:ENST0...
2531429Xhavanaexon4119625141196359.+.Parent=transcript:ENST00000462850;
Name=ENSE000...
122194419ensembl_havanaCDS56417405641946.+0ID=CDS:ENSP00000467423;
Parent=transcript:ENST0...
24307010havanaexon1311626713116340.+.Parent=transcript:ENST00000378764;
Name=ENSE000...
24135838ensembl_havanaexon144359184144359423.+.Parent=transcript:ENST00000530047;
Name=ENSE000...
21604966havanaexon111322569111322678.-.Parent=transcript:ENST00000434009;
Name=ENSE000...
83995215havanaexon7622771376227897.-.Parent=transcript:ENST00000565910;
Name=ENSE000...
95778216ensembl_havanaexon6754165367541782.+.Parent=transcript:ENST00000379312;
Name=ENSE000...
163297921ensembl_havanaexon3784065837840709.-.Parent=transcript:ENST00000609713;
Name=ENSE000...
19533994havanaexon165464390165464586.+.Parent=transcript:ENST00000511992;
Name=ENSE000...
In [123]: edf.type.value_counts()
Out[123]:
exon1180596
CDS704604
five_prime_UTR142387
three_prime_UTR133938
transcript96375
gene42470
processed_transcript28228
...
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 reRE_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]:
seqidsourcetypestartend score strand phasegene_id gene_namedesc
161havanagene1186914409.+.ENSG00000223972DDX11L1DEAD/H-box helicase 11 like 1 [Source:HGNC Sym...
281havanagene1440429570.-.ENSG00000227232WASH7PWAS protein family homolog 7 pseudogene [Sourc...
711havanagene5247353312.+.ENSG00000268020OR4G4Polfactory receptor family 4 subfamily G member...
741havanagene6294863887.+.ENSG00000240361OR4G11Polfactory receptor family 4 subfamily G member...
771ensembl_havanagene6909170008.+.ENSG00000186092OR4F5olfactory 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
SCARNA207
SCARNA166
SCARNA175
SCARNA154
SCARNA214
SCARNA114
Clostridiales-13
SCARNA43
C1QTNF9B-AS12
C11orf712
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]:
seqidsourcetypestartend score strand phasegene_id gene_namedesc
1793991ensemblgene171768070171768175.+.ENSG00000253060SCARNA20Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601]
2010371ensemblgene204727991204728106.+.ENSG00000251861SCARNA20Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601]
34920311ensemblgene85550168555146.+.ENSG00000252778SCARNA20Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601]
71852014ensemblgene6347927263479413.+.ENSG00000252800SCARNA20Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601]
83723315ensemblgene7512153675121666.-.ENSG00000252722SCARNA20Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601]
103987417ensemblgene2801877028018907.+.ENSG00000251818SCARNA20Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601]
110821517ensemblgene6023151660231646.-.ENSG00000252577SCARNA20small 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]:
count4.247000e+04
mean3.583348e+04
std9.683485e+04
min8.000000e+00
25%8.840000e+02
50%5.170500e+03
75%3.055200e+04
max2.304997e+06
Name: length, dtype: float64
.describe()根据长度值计算一些简单的统计信息：
- 基因的平均长度约为36, 000个碱基
- 基因的中位长度约为5, 200个碱基
- 最小和最大基因长度分别约为八和230万个碱基。
import matplotlib as pltndf.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]:
seqidsourcetypestartend score strand phase gene_namegene_iddesclength
23093457ensembl_havanagene146116002148420998.+.CNTNAP2ENSG00000174469contactin associated protein-like 2 [Source:HG...2304997
24225109ensembl_havanagene831424610612723.-.PTPRDENSG00000153707protein tyrosine phosphatase%2C receptor type ...2298478
2527169Xensembl_havanagene3109767733339441.-.DMDENSG00000198947dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928]2241765
44088611ensembl_havanagene8345501285627922.-.DLG2ENSG00000150672discs large MAGUK scaffold protein 2 [Source:H...2172911
23234578ensembl_havanagene29353534994972.-.CSMD1ENSG00000183117CUB and Sushi multiple domains 1 [Source:HGNC ...2059620
156991420ensembl_havanagene1399536916053197.+.MACROD2ENSG00000172264MACRO domain containing 2 [Source:HGNC Symbol%...2057829
如你所见, 最长的基因被命名为CNTNAP2, 它是与contactin相关的蛋白样2的简称。根据其Wikipedia页面,
该基因几乎涵盖了7号染色体的1.6％, 是人类基因组中最大的基因之一。确实！我们只是验证了自己。相反, 最小的基因呢？事实证明, 它们可以短至八个碱基。
In [40]: ndf.sort_values('length').head()
Out[40]:
seqidsourcetypestartend score strand phasegene_namegene_iddesclength
68227814havanagene2243854722438554.+.TRDD1ENSG00000223997T cell receptor delta diversity 1 [Source:HGNC...8
68228214havanagene2243900722439015.+.TRDD2ENSG00000237235T cell receptor delta diversity 2 [Source:HGNC...9
23068367havanagene142786213142786224.+.TRBD1ENSG00000282431T cell receptor beta diversity 1 [Source:HGNC ...12
68228614havanagene2244911322449125.+.TRDD3ENSG00000228985T cell receptor delta diversity 3 [Source:HGNC...13
18796254havanagene1023821310238235.-.AC006499.9ENSG0000027154423
两种极端情况的长度相差五个数量级(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
13902
22806
112561
192412
172280
32204
62154
122140
72106
52002
161881
X1852
41751
91659
81628
101600
151476
141449
22996
20965
13872
18766
21541
Y436
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 sourcetypestartend score strand phaseattributes
2514003MTinsdcgene6481601.+.ID=gene:ENSG00000211459;
Name=MT-RNR1;
biotype=M...
2514009MTinsdcgene16713229.+.ID=gene:ENSG00000210082;
Name=MT-RNR2;
biotype=M...
2514016MTinsdcgene33074262.+.ID=gene:ENSG00000198888;
Name=MT-ND1;
biotype=pr...
2514029MTinsdcgene44705511.+.ID=gene:ENSG00000198763;
Name=MT-ND2;
biotype=pr...
2514048MTinsdcgene59047445.+.ID=gene:ENSG00000198804;
Name=MT-CO1;
biotype=pr...
2514058MTinsdcgene75868269.+.ID=gene:ENSG00000198712;
Name=MT-CO2;
biotype=pr...
2514065MTinsdcgene83668572.+.ID=gene:ENSG00000228253;
Name=MT-ATP8;
biotype=p...
2514069MTinsdcgene85279207.+.ID=gene:ENSG00000198899;
Name=MT-ATP6;
biotype=p...
2514073MTinsdcgene92079990.+.ID=gene:ENSG00000198938;
Name=MT-CO3;
biotype=pr...
2514080MTinsdcgene1005910404.+.ID=gene:ENSG00000198840;
Name=MT-ND3;
biotype=pr...
2514087MTinsdcgene1047010766.+.ID=gene:ENSG00000212907;
Name=MT-ND4L;
biotype=p...
2514091MTinsdcgene1076012137.+.ID=gene:ENSG00000198886;
Name=MT-ND4;
biotype=pr...
2514104MTinsdcgene1233714148.+.ID=gene:ENSG00000198786;
Name=MT-ND5;
biotype=pr...
2514108MTinsdcgene1414914673.-.ID=gene:ENSG00000198695;
Name=MT-ND6;
biotype=pr...
2514115MTinsdcgene1474715887.+.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]:
seqidsourcetypelength
01GRCh38chromosome248956422
13646412GRCh38chromosome242193529
17058553GRCh38chromosome198295559
18645674GRCh38chromosome190214555
19649215GRCh38chromosome181538259
20801486GRCh38chromosome170805979
21969817GRCh38chromosome159345973
2514125XGRCh38chromosome156040895
23213618GRCh38chromosome145138636
24165609GRCh38chromosome138394717
32893811GRCh38chromosome135086622
23506810GRCh38chromosome133797422
48337012GRCh38chromosome133275309
63448613GRCh38chromosome114364328
67476714GRCh38chromosome107043718
76731215GRCh38chromosome101991189
86505316GRCh38chromosome90338345
99081017GRCh38chromosome83257441
115597718GRCh38chromosome80373285
155914420GRCh38chromosome64444167
120156119GRCh38chromosome58617616
2594560YGRCh38chromosome54106423
164748222GRCh38chromosome50818468
161671021GRCh38chromosome46709983
2513999MTGRCh38chromosome16569
为了清楚起见, 删除了与分析无关的列。
是的, .drop还可以采用一列列表并将其完全删除。
注意, 带有MT序列的行仍然存在；我们稍后会再讲。我们将执行的下一个操作是基于seqid的值合并两个数据集。
In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index()
In [75]: cdf.head(2)
Out[75]:
seqidgene_count
013902
122806
In [78]: merged = gdf.merge(cdf, on='seqid')
In [79]: merged
Out[79]:
seqidsourcetypelengthgene_count
01GRCh38chromosome2489564223902
110GRCh38chromosome1337974221600
211GRCh38chromosome1350866222561
312GRCh38chromosome1332753092140
413GRCh38chromosome114364328872
514GRCh38chromosome1070437181449
615GRCh38chromosome1019911891476
716GRCh38chromosome903383451881
817GRCh38chromosome832574412280
918GRCh38chromosome80373285766
1019GRCh38chromosome586176162412
112GRCh38chromosome2421935292806
1220GRCh38chromosome64444167965
1321GRCh38chromosome46709983541
1422GRCh38chromosome50818468996
153GRCh38chromosome1982955592204
164GRCh38chromosome1902145551751
175GRCh38chromosome1815382592002
186GRCh38chromosome1708059792154
197GRCh38chromosome1593459732106
208GRCh38chromosome1451386361628
219GRCh38chromosome1383947171659
22XGRCh38chromosome1560408951852
23YGRCh38chromosome54106423436
由于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]:
lengthgene_count
length1.0000000.728221
gene_count0.7282211.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个基因。
- 一个基因的长度范围可以从几十个到超过两百万个碱基。
- 基因在染色体之间分布不均。总体而言, 染色体越大, 其宿主的基因越多, 但是对于染色体的一个子集, 相关性可能为负。
- 一个基因通常有多少个转录本？超过1个转录本的基因百分比是多少？
- 成绩单通常有几个同工型？
- 笔录通常有多少个外显子, CDS和UTR？它们是什么尺寸？
- 是否可以根据描述栏中描述的功能对基因进行分类？
推荐阅读
- 使用R提升数据处理能力
- JSON中的双向关系支持
- 如何修复Chrome和Edge上的RESULT_CODE_HUNG错误（解决办法）
- 如何修复DX11 Feature Level 10.0需要运行引擎错误（解决办法）
- 如何修复Firefox右键单击不起作用（解决办法教程）
- 如何修复Windows 10亮度不工作（解决办法分步教程）
- 如何修复Windows 10 0xC00D36D5未连接摄像头（解决办法）
- 如何修复Windows 10 Firefox没有声音的问题（解决办法）
- App Technical Support For "Just diary"