基于python的dea绘图（dea模型python）

本文目录一览：

1、python的作图包有哪些
2、Python如何运用matplotlib库绘制3D图形
3、Python作图程序
4、Python气象数据处理与绘图(18)：泰勒图
5、如何使用Python绘制GWAS分析中的曼哈顿图和QQ图
6、Python气象数据处理与绘图(2)：常用数据计算方法

python的作图包有哪些

1. matplotlib

该python绘图包与matlab的绘图功能类似

2. seaborn

用来进行统计数据可视化的工具包，绘制的图像非常漂亮，该python package也是基于matplotlib的，是纯粹由python开发的

Python如何运用matplotlib库绘制3D图形

3D图形在数据分析、数据建模、图形和图像处理等领域中都有着广泛的应用，下面将给大家介绍一下如何在Python中使用 matplotlib进行3D图形的绘制，包括3D散点、3D表面、3D轮廓、3D直线（曲线）以及3D文字等的绘制。

准备工作：

python中绘制3D图形，依旧使用常用的绘图模块matplotlib，但需要安装mpl_toolkits工具包，安装方法如下：windows命令行进入到python安装目录下的Scripts文件夹下，执行： pip install --upgrade matplotlib即可；Linux环境下直接执行该命令。

安装好这个模块后，即可调用mpl_tookits下的mplot3d类进行3D图形的绘制。

下面以实例进行说明。

1、3D表面形状的绘制

这段代码是绘制一个3D的椭球表面，结果如下：

2、3D直线（曲线）的绘制

这段代码用于绘制一个螺旋状3D曲线，结果如下：

3、绘制3D轮廓

绘制结果如下：

Python作图程序

实战小程序：画出y=x^3的散点图

样例代码如下：

[python] view plain copy

#coding=utf-8

import pylab as y #引入pylab模块

x = y.np.linspace(-10, 10, 100) #设置x横坐标范围和点数

y.plot(x, x*x*x,'or') #生成图像

ax = y.gca()

ax.spines['right'].set_color('none')

ax.spines['top'].set_color('none')

ax.xaxis.set_ticks_position('bottom')

ax.spines['bottom'].set_position(('data', 0))

ax.yaxis.set_ticks_position('left')

ax.spines['left'].set_position(('data', 0))

ax.set_yticks([-1000, -500, 500, 1000])

y.xlim(x.min() , x.max() ) #将横坐标设置为x的最大值和最小值

y.show() #显示图像

[python] view plain copy

import pylab as y

程序中引入的pylab属于matplotlib的一个模块，将其名字用y代替，其中包括了许多NumPy和pyplot模块中常用的函数，方便用户快速进行计算和绘图，十分适合在IPython交互式环境中使用。

[python] view plain copy

y.np.linspace(-10, 10, 100)

此为numpy中的一个函数，返回的是等间距的值，numpy.linspace(a,b,c):a指的是开始位置，b表示的是结束位置，c表示产生点的个数（默认为50）

举例：

[python] view plain copy

np.linspace(2.0, 3.0, num=5)

array([ 2. , 2.25, 2.5 , 2.75, 3. ])

[python] view plain copy

y.plot(x, x*x*x,'or') #生成图像

后面加上‘o'表示为散点图

'r'可设置颜色为红色，基本上和matlab的操作很像。

[python] view plain copy

y.xlim(x.min(), x.max())

这条语句使用了xlim函数，将横坐标设置为x的大小

Python气象数据处理与绘图(18)：泰勒图

泰勒图绘制的核心思想是设计一个只有第一象限的极坐标，并将方差，相关系数进行捆绑，通过转化为极坐标系坐标进行绘制。为了实现泰勒图的绘制，我设计了两个函数：

set_tayloraxes(fig, location=111) 和plot_taylor(axes, refsample, sample, args, *kwargs)

set_tayloraxes()函数用于建立一个泰勒图的坐标系，这个自定义函数一般情况下不建议修改，每一个参数都是经过多次调试得到的，很可能牵一发动全身。因此，将绘图部分的独立成为了plot_taylor函数()，这部分函数较为简单，目的就是将需要绘图的数据，转换为极坐标系坐标，通过plot函数将散点打在泰勒图上，这个函数模块较为简单，可以根据自己的输入数据情况进行调整。

下面介绍下函数的具体用法：

输入：

fig：需要绘图的figure

rect：图的位置，如111为1行1列第一个，122为1行2列第2个

输出：

polar_ax：泰勒坐标系

输入：

axes : setup_axes返回的泰勒坐标系

refsample ：参照样本

sample ：评估样本

args, *kwargs ：plt.plot()函数的相关参数，设置点的颜色，形状等等。

下面给出示例：

如何使用Python绘制GWAS分析中的曼哈顿图和QQ图

曼哈顿图和QQ图是两个在全基因组关联（GWAS）分析里面最常出现的图形，基本上已经是GWAS的标配，几乎在每篇GWAS的文章都会见到，它们的作用和所要传达出来的信息我也在上一篇关于GWAS的博文中做了些说明，在这里我们就只集中在如何用Python和geneview将其有效地展现出来。

首先，准备一些数据来作为例子。

我这里用来展现的数据是2011年丹麦人所做过的一个关于年轻人过度肥胖的GWAS研究——GOYA，数据也是从他们所发表的结果中获得，总共有5,373个样本，其中超重的个体（case）有2,633个，正常的个体（control）是2,740个，从样本量上看还算可以。为了方便使用，我对其做了相关的处理，包括从PED和MAP文件到GEN文件的生成，并重复了一次case-control的关联性分析，计算出了芯片上所研究的各个SNP位点与肥胖相关的显著性程度（即p-value），最后又将结果数据抽取出来做成数据集——放在这里供下载（15.6Mb，csv格式）。

【注】以上内容虽提及到了一些领域内术语和相关文件格式，但若不懂也请不必纠结，因为后续处理都是基于这个最终的数据集来完成的。

接着，需要将geneview软件包加入到你的Python中，有多种不同的方式，但推荐直接使用pip，以下是安装比较稳定的发布版，直接在终端命令行下(Linux or Mac)输入：

pip install geneview

第三种办法就是直接下载源码，然后自行编译，虽然不推荐这种做法（因为还有依赖包必须自行下载安装，过程会比较麻烦低效），但对于某些不能连接外网的集群也只能如此，这三种方式都是可行的。

曼哈顿图

将示例数据下载下来：

wget data/master/GOYA.csv

先简单地查看一下数据的格式:

chrID,rsID,position,pvalue

1,rs3094315,742429,0.144586

1,rs3115860,743268,0.230022

1,rs12562034,758311,0.644366

1,rs12124819,766409,0.146269

1,rs4475691,836671,0.458197

1,rs28705211,890368,0.362731

1,rs13303118,908247,0.22912

1,rs9777703,918699,0.37948

1,rs3121567,933331,0.440824

一共是4列（逗号分隔），分别为：[1]染色体编号，[2]SNP rs 编号，[3] 位点在染色体上的位置，[4]显著性差异程度（pvalue）。在本例曼哈顿图中我们只需要使用第1,3和4列；而QQ图则只需要第4列——pvalue。

下面先从绘制曼哈顿图开始。我们先将需要的数据读取到一个列表中，可以这样做：

import csv

data = []

with open("GOYA.csv") as f:

f_csv = csv.reader(f)

headers = next(f_csv)

data = [[row[0], int(row[2]), float(row[3])] for row in f_csv]

现在GOYA.csv中的数据就都存放在data列表中了，由于Python在读取文件中数据时，都是以string类型存放，因此对于第3和第4列的数据有必要事先把做点类型转换。

接下来，调用geneview中的曼哈顿图函数。

import matplotlib.pyplot as plt

from geneview.gwas import manhattanplot

ax = manhattanplot(data, xlabel="Chromosome", ylabel="-Log10(P-value)") # 这就是Manhattan plot的函数

plt.show()

只需这样的一句代码就能创建一个漂亮的曼哈顿图，有必要再次指出的是，geneview是以matplotlib为基础开发出来的，所创建的图形对象实际上仍属于matplotlib，geneview内部自定义了很多图形风格，同时封装了大量只属于基因组数据的图表类型，但图形的输出格式以及界面显示都仍和matplotlib一样，因此在这里我们使用matplotlib.pyplot的show()函数(上例中：plt.show())将所绘制出来的曼哈顿图显示出来。如果要将图形保存下来，则只需执行`plt.savefig("man.png")`，这样就会在该目录下生成一个名为『man.png』png格式的曼哈顿图，若是要存为pdf格式，则只需将所要保存的文件名后缀改成『.pdf』（plt.savefig("man.pdf")）就可以了。下面这些格式：emf,

eps, pdf, png, jpg, ps, raw, rgba, svg,

svgz等都是支持的，至于最新的还有多少种，还请参照matplotlib文档中说明。

此外，geneview中的每个画图函数都有着足够的灵活性，我们也可以根据自己的需要做一些调整，比如：

xtick = ['1', '2','3','4','5','6','7','8','9','10','11','12','13','14','16','18', '20','22']

manhattanplot(data,

xlabel="Chromosome", # 设置x轴名字

ylabel="-Log10(P-value)", # 设置y轴名字

xtick_label_set = set(xtick), # 限定横坐标轴上的刻度显示

s=40, # 设置图中散点的大小

alpha=0.5, # 调整散点透明度

color="#f28b1e,#9a0dea,#ea0dcc,#63b8ff", # 设置新的颜色组合

)

实现新的颜色组合、限定x轴上的刻度显示和散点大小的调节。甚至还可以将散点改为线：

manhattanplot(data，

xlabel="Chromosome", # 设置x轴名字

ylabel="-Log10(P-value)", # 设置y轴名字

xtick_label_set = set(xtick), # 限定横坐标轴上的刻度显示

alpha=0.5, # 调整散点透明度

color="#f28b1e,#9a0dea,#ea0dcc,#63b8ff", # 设置新的颜色组合

kind="line"

)

其它方面的调整请查看geneview文档中的相关说明。

Q-Q图

qq图只需用到上例中的pvalue那一列：

import csv

import matplotlib.pyplot as plt

from geneview.gwas import qqplot

pvalue=[]

with open("GOYA.csv") as f:

f_csv = csv.reader(f)

headers = next(f_csv)

pvalue = [float(row[3]) for row in f_csv]

ax = qqplot(pvalue, color="#00bb33", xlabel="Expected p-value(-log10)", ylabel="Observed p-value(-log10)") # Q-Q 图

plt.show()

同样，也可以根据自己的需要对改图进行相关的调整。

以上，便是如何使用Python来制作Manhattan图和QQ图的方法，geneview的集成函数简化了这样的一个过程。

Python气象数据处理与绘图(2)：常用数据计算方法

对于气象绘图来讲，第一步是对数据的处理，通过各类公式，或者统计方法将原始数据处理为目标数据。

按照气象统计课程的内容，我给出了一些常用到的统计方法的对应函数：

在计算气候态，区域平均时均要使用到求均值函数，对应NCL中的dim_average函数，在python中通常使用np.mean()函数

numpy.mean(a, axis, dtype)

假设a为[time,lat,lon]的数据，那么

需要特别注意的是，气象数据中常有缺测，在NCL中，使用求均值函数会自动略过，而在python中，当任意一数与缺测(np.nan)计算的结果均为np.nan，比如求[1,2,3,4，np.nan]的平均值，结果为np.nan

因此，当数据存在缺测数据时，通常使用np.nanmean()函数，用法同上，此时[1,2,3,4，np.nan]的平均值为(1+2+3+4)/4 = 2.5

同样的，求某数组最大最小值时也有np.nanmax(), np.nanmin()函数来补充np.max(), np.min()的不足。

其他很多np的计算函数也可以通过在前边加‘nan’来使用。

另外，

也可以直接将a中缺失值全部填充为0。

np.std(a, axis, dtype)

用法同np.mean()

在NCL中有直接求数据标准化的函数dim_standardize()

其实也就是一行的事，根据需要指定维度即可。

皮尔逊相关系数：

相关可以说是气象科研中最常用的方法之一了，numpy函数中的np.corrcoef(x, y)就可以实现相关计算。但是在这里我推荐scipy.stats中的函数来计算相关系数：

这个函数缺点和有点都很明显，优点是可以直接返回相关系数R及其P值，这避免了我们进一步计算置信度。而缺点则是该函数只支持两个一维数组的计算，也就是说当我们需要计算一个场和一个序列的相关时，我们需要循环来实现。

其中a[time,lat,lon]，b[time]

(NCL中为regcoef()函数)

同样推荐Scipy库中的stats.linregress(x,y)函数：

slop: 回归斜率

intercept：回归截距

r_value：相关系数

p_value： P值

std_err：估计标准误差

直接可以输出P值，同样省去了做置信度检验的过程，遗憾的是仍需同相关系数一样循环计算。