生物信息学 是一个跨学科领域，开发分析和理解生物数据的方法和软件工具.
更简单地说，你可以简单地把它想象成 data science for biology.
在众多类型的生物数据中，基因组学数据是分析最广泛的数据之一. 特别是随着下一代DNA测序(NGS)技术的快速发展, 基因组数据量呈指数级增长. 根据 Stephens, Zachary D等人.在美国，基因组学数据的获取以每年eb的规模进行.
在这篇文章中，我演示了一个用SciPy Stack分析人类基因组的GFF3文件的例子. 通用功能格式版本3 (GFF3)是目前存储基因组特征的标准文本文件格式. 特别是, 在这篇文章中，你将学习如何使用SciPy堆栈来回答以下关于人类基因组的问题:
最新的人类基因组GFF3文件可从 here. The README 该目录中的文件提供了该数据格式的简要描述, 并且找到了更全面的规范 here.
We will use Pandas, SciPy堆栈的一个主要组件，提供快速, flexible, 表达性的数据结构, 操纵和理解GFF3文件.
首先，让我们安装一个安装了SciPy堆栈的虚拟环境. 如果从源代码手动构建，这个过程可能会很耗时, 因为堆栈涉及许多包-其中一些依赖于外部FORTRAN或C代码. 在这里，我建议使用 Miniconda，这使得设置非常简单.
wget http://repo.continuum.io / miniconda / Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh -b
The -b bash行上的标志告诉它以批处理模式执行. 使用以上命令成功安装Miniconda, 为基因组学开创一个新的虚拟环境, 然后安装SciPy堆栈.
Mkdir -p基因组学
cd genomics
请创建-p venv python matplotlib pandas
请注意，我们只指定了这篇文章中要使用的3个包.
如果希望在SciPy堆栈中列出所有包，只需将它们附加到 conda create command.
如果您不确定包的确切名称，请尝试 conda search. 让我们激活虚拟环境并启动ippython.
源激活venv/
ipython
IPython 是否有比默认值更强大的替代品 Python解释器接口因此，您过去在默认python解释器中所做的一切也可以在IPython中完成. 我强烈建议所有还没有使用过IPython的Python程序员都去尝试一下.
现在我们的设置已经完成，让我们下载GFF3格式的人类基因组注释文件.
大约37mb, 与人类基因组的信息量相比，这是一个非常小的文件, 纯文本大约是3gb. 这是因为GFF3文件只包含序列的注释, 而序列数据通常以另一种文件格式存储，称为 FASTA. 如果您感兴趣，可以下载FASTA here，但我们不会在本教程中使用序列数据.
!wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz
The prefixed ! 告诉IPython这是一个shell命令而不是Python命令. 但是，IPython还可以处理一些常用的shell命令，例如 ls, pwd, rm, mkdir, rmdir 即使没有前缀 !.
看一下GFF文件的头部, 您将看到许多元数据/pragmas/指令行以 ## or #!.
根据 README, ## 表示元数据是稳定的，而 #! 意味着它是实验性的.
稍后你也会看到 ###，这是另一个指令，具有更微妙的含义 规范.
人类可读的注释应该在单个 #. 为简单起见，我们将处理所有以。开头的行 # 作为评论，并在我们的分析过程中简单地忽略它们.
# # gff-version 3
##sequence-region 1 1 248956422
##sequence-region 10 1 133797422
##sequence-region 11 1 135086622
##sequence-region 12 1 133275309
...
##序列区域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)是一个国际联盟, 负责监督几个参考基因组组装的更新和改进, 包括人类, mouse, zebrafish, and chicken.
浏览这个文件，下面是前几行注释.
1 GRCh38染色体1 248956422 . . . 染色体:ID = = CM000663 1;别名.2 chr1 NC_000001.11
###
1 . 生物区10469 11240.3e+03 . . external_name=oe %3D 0.79; logic_name = cpg
1 . 生物区域10650 10657 0.999 + . 爱潘妮logic_name =
1 . 生物区10655 10657 0.999 - . 爱潘妮logic_name =
1 . 生物区10678 10687 0.999 + . 爱潘妮logic_name =
1 . 生物区10681 10688 0.999 - . 爱潘妮logic_name =
...
列是 seqid, source, type, start, end, score, strand, phase, attributes. 其中一些很容易理解. 以第一行为例:
1 GRCh38染色体1 248956422 . . . 染色体:ID = = CM000663 1;别名.2 chr1 NC_000001.11
这是第一条染色体的注释，它的种子数为1, 哪个是从一垒到24垒,895,622nd base.
换句话说，第一条染色体大约有2500万个碱基长.
我们的分析不需要来自值为的三列的信息 . (i.e. Score, strand和phase)，所以我们现在可以简单地忽略它们.
最后一个属性列显示1号染色体也有三个别名，即CM000663.2、chr1和NC_000001.11. 这基本上就是GFF3文件的样子, 但我们不会一行一行地检查它们, 现在是时候将整个文件加载到Pandas中了.
Pandas非常适合处理GFF3格式，因为它是一个以制表符分隔的文件, Pandas对读取类似csv的文件有很好的支持.
注意制表符分隔格式的一个例外是当GFF3包含 ##FASTA .
根据 规范, ##FASTA 指示注释部分的结束, 然后是一个或多个FASTA(非制表符分隔)格式的序列. 但对于我们将要分析的GFF3文件，情况并非如此.
在[1]中:将pandas导入为pd
In [2]: pd.__version__
Out[2]: '0.18.1'
In [3]: col_names = ['seqid'], 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', “属性”)
Out[3]: df = pd.read_csv(“Homo_sapiens.GRCh38.85.gff3.广州”,压缩= ' gzip ',
sep='\t'， comment='#'， low_memory=False，
头=没有名字= col_names)
上面的最后一行将加载整个GFF3文件 pandas.read_csv method.
由于它不是标准的CSV文件，因此我们需要对调用进行一些定制.
首先，我们通知Pandas在GFF3中头信息不可用 header=None，然后使用。指定每个列的确切名称 名称= col_names.
If the names 参数时，Pandas将使用增量数字作为每个列的名称.
sep='\t' 告诉Pandas列是tab分隔的，而不是逗号分隔的. 的值 sep= 实际上可以是正则表达式(regex). 如果手头的文件对每个列使用不同的分隔符(哦，是的, 出现这种情况). comment='#' 表示以 # 被考虑的评论和将被忽略.
压缩= ' gzip ' 告诉Pandas输入文件是gzip压缩的.
In addition, pandas.read_csv 有一组丰富的参数，允许读取不同类型的类似csv的文件格式.
返回值的类型是a DataFrame，这是Pandas中最重要的数据结构，用于表示2D数据.
熊猫也有一个 Series and Panel 分别为一维和三维数据的数据结构. 请参阅 文档 以了解Pandas的数据结构.
让我们看一下前几个条目 .head method.
In [18]: df.head()
Out[18]:
Seqid源类型开始结束得分链相位属性
1 GRCh38染色体1 248956422 . . . 染色体:ID = = CM000663 1;别名.2 chr1 NC_00000...
1 1 . 生物区10469 11240.3e+03 . . external_name=oe %3D 0.79; logic_name = cpg
2 1 . 生物区域10650 10657 0.999 + . 爱潘妮logic_name =
3 1 . 生物区10655 10657 0.999 - . 爱潘妮logic_name =
4 1 . 生物区10678 10687 0.999 + . 爱潘妮logic_name =
输出很好地格式化为表格格式，属性列中的较长字符串部分替换为 ....
可以将Pandas设置为不省略长字符串 pd.set_option(“显示.max_colwidth ', 1). 此外，熊猫还有很多 options 可以定制.
方法获取关于该数据框的一些基本信息 .info method.
In [20]: df.info()
RangeIndex: 2601849个条目，0 ~ 2601848
数据列(共9列):
seqid对象
源对象
类型 对象
开始int64
结束 int64
评分对象
链对象
阶段对象
属性对象
Dtypes: int64(2)， object(7)
内存使用率:178.7+ MB
这表明GFF3有2,601,848行带注释的行，每行有9列.
对于每一列，它还显示了它们的数据类型.
起始和结束是int64型，代表基因组位置的整数.
其他列都是type object，这可能意味着它们的值由整数、浮点数和字符串混合组成.
所有信息的大小大约是178.7+ MB存储在内存中. 事实证明，这比未压缩文件更紧凑，未压缩文件约为402 MB. 下面显示了一个快速验证.
- hom_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && Du -s /tmp/tmp.gff3
402 / tmp / tmp.gff3
从高层次的角度来看, 我们已经将整个GFF3文件加载到Python中的DataFrame对象中, 我们接下来的所有分析都将基于这个物体.
现在，让我们看看第一列seqid是关于什么的.
In [29]: df.seqid.unique()
Out[29]:
阵列([' 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 =对象)
In [30]: df.seqid.unique().shape
出[30]:(194)
df.seqid 是否有一种从数据框架访问列数据的方法. 另一种方法是 df['seqid']，这是更通用的语法，因为如果列名是Python保留关键字(例如.g. class)或包含 . 或者空格字符，第一种方式(df.seqid)行不通.
输出结果显示有194个唯一的seqids, 包括1号到22号染色体, X, Y, 线粒体(MT) DNA以及其他169个序列.
以KI和GL开头的序列是基因组中尚未成功组装到染色体中的DNA序列或支架.
对于那些不熟悉基因组学的人来说，这很重要.
尽管第一个人类基因组草图在15年前就已经出来了, 目前的人类基因组仍然不完整. 组装这些序列的困难主要是由于 基因组中复杂的重复区域.
接下来，让我们看一下源列.
The README 表示源是一个自由文本限定符，用于描述生成此功能的算法或操作过程.
In [66]: df.source.value_counts ()
Out[66]:
哈瓦那 1441093
ensembl_havana 745065
运用 228212年
. 182510
mirbase 4701年
GRCh38 194年
insdc 74年
这是使用的一个例子 value_counts 方法，该方法对于快速计数分类变量非常有用.
从结果来看, 我们可以看到这一列有7个可能的值, GFF3文件中的大部分条目来自哈瓦那, ensemble和ensembl_havana.
您可以了解更多关于这些来源的含义以及它们之间的关系 this post.
为了让事情简单化, 我们将重点关注GRCh38来源的条目, havana, ensembl, 和ensembl_havan.a.
每条完整染色体的信息都在GRCh38源的条目中, 我们先把剩下的过滤掉, 并将过滤后的结果赋值给一个新变量 gdf.
在[70]中:gdf = df[df].source == 'GRCh38']
In [87]: gdf.shape
[87]: [194,9]
In [84]: gdf.sample(10)
Out[84]:
Seqid源类型开始结束得分链相位属性
2511585 KI270708.1 GRCh38 supercong 1 127682 . . . ID = supercontig: KI270708.= chr1_KI270708v 1;别名...
2510840 GL000208.1 GRCh38超结构1 92689 . . . ID = supercontig: GL000208.= chr5_GL000208v 1;别名...
990810 17 GRCh38染色体1 83257441 . . . = CM000679 ID =染色体:17;别名.2 chr17 NC_000...
2511481 KI270373.1 GRCh38超结构1 1451 . . . ID = supercontig: KI270373.= chrUn_KI270373 1;别名...
2511490 KI270384.1 GRCh38超结构1 1658 . . . ID = supercontig: KI270384.= chrUn_KI270384 1;别名...
2080148 6 GRCh38染色体1 170805979 . . . = CM000668 ID =染色体:6;别名.2 chr6 NC_00000...
2511504 KI270412.1 GRCh38超结构1 1179 . . . ID = supercontig: KI270412.= chrUn_KI270412 1;别名...
1201561 19 GRCh38染色体1 58617616 . . . = CM000681 ID =染色体:19;别名.2 chr19 NC_000...
2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID = supercontig: KI270340.= chrUn_KI270340 1;别名...
2594560 Y GRCh38染色体2781480 56887902 . . . = CM000686 ID =染色体:Y;别名.2 chrY NC_00002...
在Pandas中过滤很容易.
如果检查从表达式求出的值 df.source == 'GRCh38'，它是一系列 True and False 索引相同的每个条目的值 df. 将其传递给 df[] 将只返回那些其对应值为True的条目.
一共有194把钥匙 df[] for which df.source == 'GRCh38'.
中也有194个唯一的值 seqid 列，表示每一项 gdf 对应于一个特定的种子.
然后我们随机选择10个条目 sample 方法要仔细看看.
你可以看到，未组装的序列是超contig型，而其他的是染色体型. 来计算基因组中不完整的部分, 我们首先需要知道整个基因组的长度, 哪个是所有序列长度的和.
[90]中:gdf = gdf.copy()
In [91]: gdf['length'] = gdf.end - gdf.start + 1
In [93]: gdf.head()
Out[93]:
Seqid源类型开始结束得分链相位属性长度
1 GRCh38染色体1 248956422 . . . 染色体:ID = = CM000663 1;别名.2 chr1 NC_00000... 248956421
235068 10 GRCh38染色体1 133797422 . . . 染色体:ID = = CM000672 10;别名.2 chr10 NC_000... 133797421
328938 11 GRCh38染色体1 135086622 . . . = CM000673 ID =染色体:11;别名.2 chr11 NC_000... 135086621
483370 12 GRCh38染色体1 133275309 . . . = CM000674 ID =染色体:12;别名.2 chr12 NC_000... 133275308
13 GRCh38染色体1 114364328 . . . 染色体:ID = = CM000675 13;别名.2 chr13 NC_000... 114364327
In [97]: gdf.length.sum()
出[97]:3096629532
在[99]:装备= (str(_) _的范围(23)]+[‘X’,‘Y’,‘太’)
In [101]: gdf[-gdf].seqid.isin(chrs)].length.sum() / gdf.length.sum()
Out[101]: 0.0037021917421198327
在上面的代码片段中，我们首先复制了 gdf with .copy(). 否则，原始的 gdf 只是一片 df，直接修改会导致 SettingWithCopyWarning (see here 欲知详情).
然后计算每个条目的长度，并将其加回 gdf 作为一个名为length的新列. 总长度大约是3.10亿，未组装序列的比例约为0.37%.
下面是最后两个命令中切片的工作原理.
First, 我们创建了一个字符串列表，它涵盖了所有良好组合序列的序列, 哪些都是染色体和线粒体. 然后我们使用 isin 方法过滤序号在 chrs list.
负号(-)被添加到索引的开头以反转选择, 因为我们想要所有的东西 not 在列表中(i ..e. 我们想要从KI和GL开始的未组装的…
注:由于已组装和未组装序列是通过类型列来区分的, 最后一行也可以重写如下，以获得相同的结果.
Gdf [(Gdf ['type'] == 'supercontig')].length.sum() / gdf.length.sum()
在这里，我们关注来自源集成的条目, 和ensembl_havana，因为它们是大多数注释条目所属的地方.
在[109]中:edf = df[df].source.Isin (['ensembl'， 'havana'， 'ensembl_havana'])]
[111];编辑.sample(10)
Out[111]:
Seqid源类型开始结束得分链相位属性
915996 16哈瓦那CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0...
2531429 X哈瓦那外显子41196251 41196359 . + . 父母=成绩单:ENST00000462850; Name = ENSE000...
1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0...
243070 10哈瓦那外显子13116267 13116340 . + . 父母=成绩单:ENST00000378764; Name = ENSE000...
2413583 8 ensembl_havana外显子144359184 144359423 . + . 父母=成绩单:ENST00000530047; Name = ENSE000...
2160496哈瓦那外显子111322569 111322678 . - . 父母=成绩单:ENST00000434009; Name = ENSE000...
8399952 15哈瓦那外显子76227713 76227897 . - . 父母=成绩单:ENST00000565910; Name = ENSE000...
957782 16 ensembl_havana外显子67541653 67541782 . + . 父母=成绩单:ENST00000379312; Name = ENSE000...
1632979 21 ensembl_havana外显子37840658 37840709 . - . 父母=成绩单:ENST00000609713; Name = ENSE000...
1953399 4哈瓦那外显子165464390 165464586 . + . 父母=成绩单:ENST00000511992; Name = ENSE000...
[123];译.type.value_counts ()
Out[123]:
外显子 1180596
cd 704604年
five_prime_UTR 142387年
three_prime_UTR 133938年
成绩单 96375年
基因 42470年
processed_transcript 28228年
...
名称:type, dtype: int64
The isin 方法再次用于过滤.
Then, 快速值计数显示，大多数条目是外显子, 编码序列, 和未翻译区(UTR).
这些都是亚基因元素，但我们主要是在寻找基因计数. 如图所示，有42,470个，但我们想知道更多.
具体来说，他们的名字是什么，他们做什么? 要回答这些问题，我们需要仔细查看属性列中的信息.
In [127]: ndf = edf[edf ..类型== '基因']
In [173]: ndf = ndf.copy()
In [133]: ndf.sample(10).attributes.values
Out[133]:
array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=肝细胞核因子4 γ伪基因1[来源: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=核糖体蛋白S6激酶A3[来源:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12',
“ID =基因:ENSG00000231748; Name = rp11 - 227 h15.5、生物型=反义gene_id = ENSG00000231748; havana_gene = OTTHUMG00000018373; havana_version = 1; logic_name =哈瓦那;version = 1 ',
'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1受体33伪基因[来源: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[来源: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=锌指DHHC-type containing 22[来源: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=小Cajal体特异性RNA22[来源: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[来源:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16',
“ID =基因:ENSG00000229224; Name = AC105398.3;生物型=反义;gene_id = ENSG00000229224; havana_gene = OTTHUMG00000152025; havana_version = 1; logic_name =哈瓦那;version = 1 ',
'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=淋巴细胞抗原6复合物%2C位点G6E(假基因)[来源:HGNC符号%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype =对象)
它们被格式化为以分号分隔的标记值对列表. 我们最感兴趣的信息是基因名称, 基因识别与描述, 我们将用正则表达式(regex)提取它们.
import re
RE_GENE_NAME = re.编译(r 'Name = (?P.+?);')
def extract_gene_name (attributes_str):
res = RE_GENE_NAME.搜索(attributes_str)
return res.集团(“gene_name”)
Ndf ['gene_name'] = Ndf.attributes.应用(extract_gene_name)
首先，我们提取基因名称.
在正则表达式中 Name=(?P , +? 是用来代替 + because we want it to be non-greedy and let the search stop at the first semicolon; otherwise, 结果将匹配到最后一个分号.
此外，regex首先使用 re.compile 而不是直接用在 re.search 为了更好的性能，因为我们将把它应用于数千个属性字符串.
extract_gene_name 中使用的辅助函数 pd.apply, 当需要将函数应用于数据框架或系列的每个条目时，使用哪种方法.
在这个特殊的例子中，我们想要提取每个条目的基因名称 ndf.attributes，并将这些名称添加回 ndf 在一个叫做 gene_name.
基因id和描述以类似的方式提取.
RE_GENE_ID = re.编译(r 'gene_id = (?PENSG.+?);')
def extract_gene_id (attributes_str):
res = RE_GENE_ID.搜索(attributes_str)
return res.集团(“gene_id”)
Ndf ['gene_id'] = Ndf.attributes.应用(extract_gene_id)
RE_DESC = re.编译(描述= (?P.+?);')
def extract_description (attributes_str):
res = RE_DESC.搜索(attributes_str)
如果res为None:
return ''
else:
return res.集团(desc)
Ndf ['desc'] = Ndf.attributes.应用(extract_description)
正则表达式 RE_GENE_ID 更具体一点，因为我们知道每个 gene_id 必须以 ENSG, where ENS means ensembl and G means gene.
对于没有任何描述的条目，我们将返回一个空字符串. 在一切都被提取出来之后, 我们将不再使用attributes列, 我们把它去掉，让这个方法更简洁 .drop:
文献[224]:ndf.drop('attributes'， axis=1, 原地= True)
文献[225]:ndf.head()
Out[225]:
Seqid源类型start end score strand phase gene_id gene_name desc
16哈瓦那基因11869 14409 . + . eng00000223972 DDX11L1 DEAD/H-box解旋酶11 like 1[来源:HGNC Sym...
28 1哈瓦那基因14404 29570 . - . eng00000227232 WASH7P WAS蛋白家族同源7假基因...
71哈瓦那基因52473 53312 . + . eng00000268020 OR4G4P嗅觉受体家族4亚家族G成员...
74 .哈瓦那基因62948 63887 . + . eng00000240361 OR4G11P嗅觉受体家族4亚家族G成员...
77 1 ensembl_havana基因69091 70008 . + . eng00000186092 OR4F5嗅觉受体家族4亚家族F成员...
在上述通话中， attributes 指示要删除的特定列.
axis=1 意味着我们要删除一列而不是一行(axis=0 by default).
原地= True 意味着删除是在DataFrame本身上操作的，而不是返回一个删除指定列的新副本.
A quick .head 查看显示属性列确实消失了，并且有三个新列: gene_name, gene_id, and desc 已添加.
出于好奇，我们来看看 gene_id and gene_name are unique:
参见[232]:ndf.shape
[232] [c]; [c]; [c];
In [233]: ndf.gene_id.unique().shape
[233]: (42470)
[234]: ndf.gene_name.unique().shape
[234]: (42387)
Surprisingly, 基因名称的数量小于基因id的数量, 表示某个gene_name必须对应多个基因id. 让我们来看看它们是什么.
In [243]: count_df = ndf.groupby(“gene_name”).count().ix[:, 0].sort_values ().ix[::-1]
在[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
名称:seqid，类型:int64
In [262]: count_df[count_df > 1].shape
[262]: (63)
在[263]中:count_df.shape
[263]: (42387)
In [264]: count_df[count_df > 1].形状[0]/ count_df.shape[0]
Out[264]: 0.0014863047632528842
我们根据的值对所有的项进行分组 gene_name，然后用。来计算每组的项目数 .count().
的输出 ndf.groupby(“gene_name”).count()，对每个组的所有列进行计数，但大多数列具有相同的值.
注意，计数时不会考虑NA值，因此只计算第一列的计数。 seqid ( we use .ix[:, 0] 以确保没有NA值).
然后对计数值进行排序 .sort_values 然后反过来 .ix[::-1].
因此，一个基因名称可以与多达七个基因id共享.
[255]中:ndf[ndf ..gene_name == 'SCARNA20']
Out[255]:
Seqid源类型start end score strand phase gene_id gene_name desc
179399 1 ensembl基因171768070 171768175 . + . SCARNA20 Small Cajal body specific RNA20[来源:RFAM%3BAcc:RF00601]
201037 1 ensembl基因204727991 204728106 . + . SCARNA20 Small Cajal body specific RNA20[来源:RFAM%3BAcc:RF00601]
349203 11 ensembl基因8555016 8555146 . + . SCARNA20 Small Cajal body specific RNA20[来源:RFAM%3BAcc:RF00601]
14 ensembl基因63479272 63479413 . + . eng00000252800 SCARNA20 Small Cajal body specific RNA20[来源:RFAM%3BAcc:RF00601]
837233 15 ensembl基因75121536 75121666 . - . eng00000252722 SCARNA20 Small Cajal body specific RNA20[来源:RFAM%3BAcc:RF00601]
1039874 17 ensembl基因28018770 28018907 . + . eng00000251818 SCARNA20 Small Cajal body specific RNA20[来源:RFAM%3BAcc:RF00601]
17 ensembl基因60231516 60231646 . - . eng00000252577 SCARNA20小Cajal体特异性RNA20[来源:HGNC符号%3BAcc:HGNC:32578]
仔细观察所有的SCARNA20基因就会发现它们确实都是不同的.
虽然它们有相同的名字，但它们位于基因组的不同位置.
然而，他们的描述似乎对区分它们没有太大帮助.
这里的重点是要知道，基因名称并不是所有基因id唯一的，大约为0.15%的名字是由多个基因共享的.
就像我们试图理解基因组的不完整性时所做的一样, 我们可以很容易地加上a length column to ndf:
In [277]: ndf['length'] = ndf.end - ndf.start + 1
[278].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
名称:length, dtype: float64
.describe() 根据长度值计算一些简单的统计数据:
一个基因的平均长度约为36000个碱基
基因的中位长度约为5200个碱基
最小和最大基因长度分别为8和2.分别有300万个碱基长.
因为均值比中位数大很多, 这意味着长度分布向右倾斜. 为了更具体地观察，让我们绘制分布图.
导入matplotlib为PLT
ndf.length.plot(kind='hist'， bins=50, logy=True)
plt.show()
Pandas为matplotlib提供了一个简单的接口，可以非常方便地使用dataframe或series进行绘图.
在这种情况下，它说我们想要一个直方图(kind='hist')，设y轴为对数刻度(logy=True).
从直方图中，我们可以看到大多数基因都在第一个箱子里. 然而，一些基因的长度可以超过200万个碱基. 让我们来看看它们是什么:
在[39]中:ndf[ndf ..length > 2e6].sort_values(长度).ix[::-1]
Out[39]:
Seqid源类型start end score strand phase gene_name gene_id desc length
2309345 7 ensembl_havana基因146116002 148420998 . + . CNTNAP2 eng00000174469接触蛋白相关蛋白样2[来源:HG... 2304997
2422510 9 ensembl_havana基因8314246 10612723 . - . PTPRD eng00000153707蛋白酪氨酸磷酸酶%2C受体类型 ... 2298478
2527169 X ensembl_havana基因31097677 33339441 . - . DMD ENSG00000198947肌营养不良蛋白[来源:HGNC符号%3BAcc:HGNC:2928] 2241765
440886 11 ensembl_havana基因83455012 85627922 . - . DLG2 eng00000150672圆盘大MAGUK支架蛋白2[来源:H... 2172911
2323457 8 ensembl_havana基因2935353 4994972 . - . CSMD1 ENSG00000183117 CUB和Sushi多域1[来源:HGNC ... 2059620
1569914 20 ensembl_havana基因13995369 16053197 . + . MACROD2 eng00000172264宏域包含2[来源:HGNC Symbol%... 2057829
如你所见, 最长的基因被命名为CNTNAP2, 这是接触蛋白相关蛋白2的简称. 根据它 维基百科页面,
这个基因包含近1.占7号染色体的6%，是人类基因组中最大的基因之一.
Indeed! 我们刚刚验证过了. 相比之下，最小的基因呢? 结果是它们可以短到8个碱基.
在[40]:ndf中.sort_values(长度).head()
Out[40]:
Seqid源类型start end score strand phase gene_name gene_id desc length
682278 14哈瓦那基因22438547 22438554 . + . TRDD1 eng00000223997 T细胞受体δ多样性1[来源:HGNC... 8
682282 14哈瓦那基因22439007 22439015 . + . TRDD2 ENSG00000237235 T细胞受体δ多样性2[来源:HGNC... 9
2306836 7哈瓦那基因142786213 142786224 . + . TRBD1 ENSG00000282431 T细胞受体β多样性1[来源:HGNC ... 12
682286 14哈瓦那基因22449113 22449125 . + . TRDD3 ENSG00000228985 T细胞受体δ多样性3[来源:HGNC... 13
1879625 4哈瓦那基因10238213 10238235 . - . AC006499.9 eng00000271544
这两个极端情况的长度相差5个数量级.300万vs. 8)，这是一个巨大的，这可以表明生命的多样性水平.
一个基因可以通过一种称为选择性剪接的过程翻译成许多不同的蛋白质, 一些我们还没有探索过的东西. 这些信息也在GFF3文件中，但超出了本文的范围.
我想讨论的最后一件事是染色体中的基因分布, 这也可以作为介绍 .merge 组合两个数据框架的方法. 直觉上，较长的染色体可能承载更多的基因. 我们来看看这是不是真的.
在[53]:ndf = ndf[ndf.seqid.isin(chrs)]
在[54]中:chr_gene_counts = ndf.groupby(“seqid”).count().ix[:, 0].sort_values ().ix[::-1]
出[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
名称:source, dtype: int64
我们借了 chrs 变量，并使用它过滤掉未组装的序列. 根据输出，最大的1号染色体确实拥有最多的基因. 虽然Y染色体的基因数量最少，但它并不是最小的染色体.
请注意，线粒体(MT)中似乎没有基因，这是不正确的.
对第一个DataFrame进行更多的过滤 df returned by pd.read_csv 结果表明，所有MT基因均来自源insc(在产生前被过滤掉) edf 我们只考虑了哈瓦那，ensembl，或ensembl_havana的来源).
In [134]: df[(df ..类型== '基因') & (df.seqid == 'MT')]
Out[134]:
Seqid源类型开始结束得分链相位属性
2514003 MT insdc基因648 1601 . + . ID =基因:ENSG00000211459; Name = MT-RNR1;生物型= M...
2514009 MT insdc基因1671 3229 . + . ID =基因:ENSG00000210082; Name = MT-RNR2;生物型= M...
MT insdc基因3307 4262 . + . ID =基因:ENSG00000198888; Name = MT-ND1;生物型=公关...
2514029 MT insdc基因4470 5511 . + . ID =基因:ENSG00000198763; Name = MT-ND2;生物型=公关...
MT insdc基因5904 7445 . + . ID =基因:ENSG00000198804; Name = MT-CO1;生物型=公关...
2514058 MT insdc基因7586 8269 . + . ID =基因:ENSG00000198712; Name = MT-CO2;生物型=公关...
2514065 MT insdc基因8366 8572 . + . ID =基因:ENSG00000228253; Name = MT-ATP8;生物型= p...
2514069 MT insdc基因85279207 . + . ID =基因:ENSG00000198899; Name = MT-ATP6;生物型= p...
2514073 MT insdc基因9207 9990 . + . ID =基因:ENSG00000198938; Name = MT-CO3;生物型=公关...
2514080 MT insdc基因10059 10404 . + . ID =基因:ENSG00000198840; Name = MT-ND3;生物型=公关...
2514087 MT insdc基因10470 10766 . + . ID =基因:ENSG00000212907; Name = MT-ND4L;生物型= p...
2514091 MT insdc基因10760 12137 . + . ID =基因:ENSG00000198886; Name = MT-ND4;生物型=公关...
2514104 MT insdc基因12337 14148 . + . ID =基因:ENSG00000198786; Name = MT-ND5;生物型=公关...
MT insdc基因14149 14673 . - . ID =基因:ENSG00000198695; Name = MT-ND6;生物型=公关...
MT insdc基因14747 15887 . + . ID =基因:ENSG00000198727; Name = MT-CYB;生物型=公关...
此示例还展示了在筛选过程中如何组合两个条件 &; the logical operator for “or” would be |.
注意，每个条件周围的括号都是必需的, Pandas中的这部分语法与Python不同, 哪个是字面意思 and and or.
接下来，我们借 gdf 来自上一节的DataFrame作为每条染色体长度的来源:
在[61]中:gdf = gdf[gdf].seqid.isin(chrs)]
文献[62]:gdf.drop(['start'， 'end'， 'score'， 'strand'， 'phase'，“属性”)， axis=1, 原地= True)
文献[63]:gdf.sort_values(长度).ix[::-1]
Out[63]:
Seqid源类型长度
0 1 GRCh38染色体248956422
2 GRCh38染色体242193529
3 GRCh38染色体198295559
4 GRCh38染色体190214555
5 GRCh38染色体181538259
2080148 6 GRCh38染色体170805979
7 GRCh38染色体159345973
X GRCh38染色体156040895
8 GRCh38染色体145138636
2416560 9 GRCh38染色体138394717
328938 11 GRCh38染色体135086622
235068 10 GRCh38染色体133797422
483370 12 GRCh38染色体133275309
13 GRCh38染色体114364328
674767 14 GRCh38染色体107043718
767312 15 GRCh38染色体101991189
865053 16 GRCh38染色体90338345
990810 17 GRCh38染色体83257441
1155977 18 GRCh38染色体80373285
1559144 20 GRCh38染色体64444167
1201561 19 GRCh38染色体58617616
Y GRCh38染色体54106423
1647482 22 GRCh38染色体50818468
21 GRCh38染色体46709983
2513999 MT GRCh38染色体16569
为了清晰起见，删除了与分析不相关的列.
Yes, .drop 还可以在一个操作中获取列列表并将它们一起删除吗.
Note that the row with a seqid of MT is still there; we will get back to it later. 我们将执行的下一个操作是基于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
[78]中:merged = gdf.合并(cdf, =“seqid”)
[79]:合并
Out[79]:
Seqid源类型长度gene_count
1 GRCh38染色体248956422 3902
1 10 GRCh38染色体133797422 1600
2 11 GRCh38染色体135086622 2561
3 12 GRCh38染色体133275309 2140
4 13 GRCh38染色体114364328 872
5 14 GRCh38染色体107043718 1449
6 15 GRCh38染色体101991189 1476
7 16 GRCh38染色体90338345 1881
8 17 GRCh38染色体83257441 2280
9 18 GRCh38染色体80373285 766
10 19 GRCh38染色体58617616 2412
11 2 GRCh38染色体242193529 2806
12 20 GRCh38染色体64444167 965
13 21 GRCh38染色体46709983 541
14 22 GRCh38染色体50818468 996
15 3 GRCh38染色体198295559 2204
164 GRCh38染色体190214555 1751
5 GRCh38染色体181538259 2002
18 6 GRCh38染色体170805979 2154
19 7 GRCh38染色体159345973 2106
8 GRCh38染色体145138636 1628
21 9 GRCh38染色体138394717 1659
22 X GRCh38染色体156040895 1852
23 Y GRCh38染色体54106423 436
Since chr_gene_counts 仍然是一个系列对象吗, 哪个不支持合并操作, 我们需要将其转换为DataFrame对象 .to_frame.
.reset_index () 转换原始索引(i.e. seqid)放入新列，并将当前索引重置为基于0的增量数字.
的输出 cdf.head(2) 显示了它的样子. 接下来，我们使用 .merge 方法组合id列上的两个DataFrame (on='seqid').
合并后 gdf and cdf, the MT 条目仍然缺失. 这是因为，默认情况下， .merge 操作内连接，而左连接、右连接或外连接可通过调优 how parameter.
请参阅 文档 欲知详情.
后来，你可能会发现，还有一个相关的 .join method. .merge and .join 相似但有不同的api.
根据官方说法 文档 says
相关的数据框.join method, 在内部对索引对索引和索引对列的连接使用merge, 但是默认情况下是在索引上连接，而不是在公共列上连接(合并的默认行为)。. 如果您要在索引上联接，您可能希望使用DataFrame.加入可以节省输入的时间.
Basically, .merge 是更通用和使用的 .join.
最后，我们准备计算染色体之间的相关性 length and gene_count.
In [81]: merged[['length'， 'gene_count']].corr()
Out[81]:
长度gene_count
长度是1.000000 0.728221
gene_count 0.728221 1.000000
By default .corr 计算 皮尔森相关 在数据框中的所有列对之间.
但是在这种情况下我们只有一对列, 相关系数是正- 0.73.
换句话说，染色体越大，就越有可能拥有更多的基因. 在按值对排序之后，我们再绘制这两列 length:
Ax = merged[['length'， 'gene_count']].sort_values(长度).Plot (x='length'， y='gene_count'， style='o-')
在x轴的两端添加一些边距
xlim = ax.get_xlim()
Margin = xlim[0] * 0.1
ax.Set_xlim ([xlim[0] - margin, xlim[1] + margin])
标记图形上的每个点
For (s, x, y) in merged[['seqid'， 'length'， 'gene_count']].sort_values(长度).values:
ax.Text (x, y - 100, str(s))
如上图所示, 尽管总体上是正相关的, 它并不适用于所有染色体. 特别是, 对于第17号染色体, 16, 15, 14, 13, 相关性实际上是负的, 这意味着染色体上的基因数量随着染色体大小的增加而减少.
这就是我们关于使用SciPy堆栈操纵GFF3格式的人类基因组注释文件的教程. 我们主要使用的工具包括IPython、Pandas和matplotlib. 在教程期间, 我们不仅学习了Pandas中一些最常见和最有用的操作, 我们还回答了一些关于我们基因组的有趣问题. In summary:
GFF3文件中包含非常丰富的注释信息，我们只是触及了表面. 如果你对进一步的探索感兴趣，这里有几个问题你可以试试:
Zhuyi是一名熟练的Python开发人员，拥有超过七年的经验. 他还精通JavaScript和Scala.
作者都是各自领域经过审查的专家，并撰写他们有经验的主题. 我们所有的内容都经过同行评审，并由同一领域的Toptal专家验证.
世界级的文章，每周发一次.
世界级的文章，每周发一次.
