Python科研统计作图Plotnine+Seaborn+matplotlib替代R ggplot2系列（2）

系列文章的第二讲主要关注更加plotnine作图更加细节的地方。首先作为科研工作者，我们很大一部分工作，无论是某个算法相比于之前的算法更加节省时间，稳定性更高，还是生命科学里面一个基因敲除之后蛋白表达水平降低，这些数据结果的呈现方式都需要按照不同的类别分组，然后再进一步呈现，如果用生命科学中的术语来说，我们需要和control进行比较，从而得出我们的结论。这就免不了需要分组做图，这样更加直观，说服力更强。The gramma of graphics这门书，也就是ggplot2的灵感来源，实际上就是告诉我们通过什么样的数据格式去展现不同情况下的数据。

首先，传统的参数统计以及方差分析等作图的方式相对简单，数据量较小，而且使用柱状图和折线图即可，简单的分组和事后比较。然而，非参数统计，或者多元统计方面，则需要一些复杂的方法，比如降维，排序等。下图以python非度量多维尺度为例，详细展示了如何使用sklearn中的MDS函数进行nmds分析并且将数据转化为dataframe，用于plotnine作图。在这个案例中，不同的样品点之间，或者不同的处理之间，可以通通过不同的形状和颜色区别开来，对于嵌套的实验设计同样可以这样做。而且同一个处理之间通过置信椭圆区别开来。这是典型的多元统计的作图展现方法。同样的可以使用PCA等降维方法展现不同处理之间的差异。

# -*- coding: utf-8 -*-"""@author: Jianshu Zhao, zhaojs16@mails.tsinghua.edu.cn    Department of Environmental Metagenomics, School of Environment, Tsinghua University"""import matplotlibimport matplotlib.pylab as pylabimport matplotlib.pyplot as pltimport pandas as pdimport numpy as npfrom biom import load_tablefrom biom import example_tablefrom biom import parse_tableimport skbiofrom sklearn import manifoldfrom sklearn.decomposition import PCAfrom scipy.spatial import distancefrom plotnine import *# Say, "the default sans-serif font is Helvetica"matplotlib.rcParams['font.sans-serif'] = "Helvetica"# Then, "ALWAYS use sans-serif fonts"matplotlib.rcParams['font.family'] = "sans-serif"params={
    'axes.labelsize': '14',
    'xtick.labelsize':'14',
    'ytick.labelsize':'14',
    'lines.linewidth':1,
    'legend.fontsize': '14',
    'figure.figsize'   : '6, 4.5'}pylab.rcParams.update(params)### using path to open data Users/jaysonchao/Documents/python/MacQiime_hawaii_bacteria/usearch61_open_ref_no_chimera/# with open('16S_otu_table_mc2_w_tax_no_singletons_nolowcoverage_sorted_rarefied_40000.biom') as f:
    # table = parse_table(f)## a function to extrac otu table from biom file without taxonomydef exploding_panda(_bt):
    """BIOM->Pandas dataframe converter    Parameters    ----------    _bt : biom.Table        BIOM table    Returns    -------    pandas.DataFrame        The BIOM table converted into a DataFrame        object.    References    ----------    Based on this answer on SO:    http://stackoverflow.com/a/17819427/379593    """
    m = _bt.matrix_data
    data = [pd.SparseSeries(m[i].toarray().ravel()) for i in np.arange(m.shape[0])]
    out = pd.SparseDataFrame(data, index=_bt.ids('observation'),
                             columns=_bt.ids('sample'))
    


    
return out### using path to open data Users/jaysonchao/Documents/python/MacQiime_hawaii_bacteria/usearch61_open_ref_no_chimera/# with open('16S_otu_table_mc2_w_tax_no_singletons_nolowcoverage_sorted_rarefied_40000.biom') as f:
    # table = parse_table(f)# direct load data from biom filetable_all = load_table('16S_otu_table_mc2_w_tax_no_singletons_nolowcoverage_sorted_rarefied_40000.biom')# table_all = load_table('otu_table.biom')otu_table = exploding_panda(table_all)otu_table_T = otu_table.T# print(otu_table_T.index)index = otu_table_T.index# center the data# otu_table_T -= otu_table_T.mean()## calculating bray-curtis distance matrixDM_dist = distance.squareform(distance.pdist(otu_table_T, metric="braycurtis")) # (n x n) distance measureDM_dist_DF = pd.DataFrame(DM_dist,index=index, columns=index)# print(DM_dist_DF)## transform to distance matrix for pcoa in skbio# DM_dist = skbio.stats.distance.DistanceMatrix(Ar_dist, ids=otu_table_T.index)# mdsseed = np.random.RandomState(seed=3)mds = manifold.MDS(n_components=2, max_iter=3000, eps=1e-9, random_state=seed,
                   dissimilarity="precomputed", n_jobs=1)pos = mds.fit(DM_dist).embedding_# nmdsnmds = manifold.MDS(n_components=2, metric=False, max_iter=3000, eps=1e-12,
                    dissimilarity="precomputed", random_state=seed, n_jobs=1,
                    n_init=1)npos = nmds.fit_transform(DM_dist, init=pos)## rotate data# npos *= np.sqrt((otu_table_T ** 2).sum()) / np.sqrt((npos ** 2).sum())clf = PCA(n_components=2)npos = clf.fit_transform(npos)# print(npos)column = ['NMDS1','NMDS2']npos_DF = pd.DataFrame(npos,index=index, columns=column)print(npos_DF)## read group data and join two dataframegroup = pd.read_csv('group.csv',header=0,index_col=0)nmds_data = pd.concat([npos_DF, group], axis=1)print(nmds_data.index)print(nmds_data.columns)print(nmds_data.head)print(nmds_data)p1 = (ggplot(nmds_data, aes(x='NMDS1', y='NMDS2'))+
	geom_point(aes(color='precipitation',shape='depth'),size=3


    
)+ theme_classic()+scale_shape_manual(['s','^','o','])+
	theme(axis_text=element_text(size=14,color='black',family='sans-serif'),axis_title=element_text(size=14,color='black',family='sans-serif'),legend_text= element_text(size=12,color='black',family='sans-serif'))+
	theme(legend_background=element_rect(alpha=0))+
	stat_ellipse(aes(group='precipitation',color='precipitation'),level = 0.75,size=0.8,geom='path',type='t')+
	scale_color_manual(values=['#56B4E9', '#E69F00', '#D55E00', '#009E73', '#F0E442', '#0072B2', '#000000', '#CC79A7','#BE5C2E','#3F725D','#2CE1C9','#944F4A','#666666','#999999'])+
	labs(color='',shape='')
	)p1.save(filename = 'nmds_bacteria.pdf', height=5.35, width=5, units = 'in')

不同于R语言中的mds函数，python机器学习包sklearn需要几步才能完成nmds的计算，plotnine对figure参数的调整和改变几乎和ggplot2相同。值得一提的是stat_ellipse()函数在一个月前刚刚更行，之前的plotnine版本中是无法使用的。plotnine包仍然在大幅度的更新，很多ggplot2具有的其他功能也正在完善，有需要的同学可以时刻关注plotnine的作者个人网页

A Grammar of Graphics for Pythonplotnine.readthedocs.io

或者github

has2k1/plotninegithub.com

另外一种方法就是可以先在R中计算出nmds的记过再导入到python中可视化，同样的在R中完成可视化也同样很方便。

在完成堆积柱状图的时候，在R中非常方便，只要将对应的数据格式转换，导入读取画图即可，plotnine同样如此，注意将所有的element.text()改成element_text()。

MyColors 

和在ggplot2中相同，bar距离x轴和y轴的距离可以通过scale_y_continuous(expand=[0,0,0.1,0])调整而不需要保留空白一看就能看出来是ggplot2画的图。当然我们也可以通过添加分面的方式，把不同的组别分开只需要加上facidwrap(~group,scale="free_x")，按照什么类型分面，可以根据自己的需求修改画图导入的metadata。

当说到分组的时候，我们通常对应不同的处理或者样地而言，有些情况下，我们不仅需要各个处理的情况，还需要知道整体，尤其是做相关性分析的时候。分组的同时也能观察并且可视化总体的情况，这是用plotnine也是很快就可以完成的。

    
"""
@author: Jianshu Zhao (zhaojs16@mails.tsinghua.edu.cn)
"""
import matplotlib
# Say, "the default sans-serif font is Helvetica"
matplotlib.rcParams['font.sans-serif'] = "Helvetica"
# Then, "ALWAYS use sans-serif fonts"
matplotlib.rcParams['font.family'] = "sans-serif"

import numpy as np
import matplotlib.pylab as pylab
import matplotlib.pyplot as plt
params={
    'axes.labelsize': '14',
    'xtick.labelsize':'14',
    'ytick.labelsize':'14',
    'lines.linewidth':1,
    'legend.fontsize': '14',
    'figure.figsize'   : '6, 5.5'

pylab.rcParams.update(params)

import math
import pandas as pd
import numpy as np
from plotnine import *
from plotnine.data import *

# read dataframe and check
bacteria_ddr_group = pd.read_csv('bacteria_abundant_rare_ddr_group.csv',header=0,index_col=0)
print(bacteria_ddr_group.head())
plot_abundant=(ggplot(bacteria_ddr_group,aes(x='geodis_log',y='bray_abundant_log'))+geom_point(aes(color='category'),size=0.5)+theme_classic()+
	stat_smooth(aes(color='category'),method='lm',se=False)+
	scale_color_manual(values=['#1AA68C', '#D55E00','#56B4E9'])+
stat_smooth(method='lm',color='grey',se=False)+
theme(legend_position=(0.3,0.3))+theme(legend_background=element_rect(alpha=0))+theme(axis_text=element_text(size=12,color='black',family='sans-serif'),axis_title=element_text(size=13,color='black',family='sans-serif'),
		legend_text= element_text(size=12,color='black',family='sans-serif'),legend_title= element_text(size=12,color='black',family='sans-serif')))
plot_abundant.save('bacter_ddr11.pdf',width=6, height=5.5)

我们可以同时按照分组和不分组的方式对x和y做相关性分析，这个和其他可视化软件比如sigmaplot相比就显得非常人性化，可以同时叠加多个图层，而在sigmaplot则可能需要添加新的数据分组。

前天的文章很不错，没看的朋友请阅读学习：

plotnine：python数据可视化版ggplot2

plotnine数据可视化手册的领取方式：

1.  将该文章转发到微信朋友圈或者微信群，获得10个以上的赞；

2. 添加微信号：EasyCharts，将转发与点赞情况截图发给该微信。

欢迎大家加入QQ群一起探讨学习

如需联系EasyCharts团队

请加微信：EasyCharts

【书籍推荐】《Excel 数据之美--科学图表与商业图表的绘制》

【手册获取】国内首款-数据可视化参考手册：专业绘图必备

【必备插件】  EasyCharts -- Excel图表插件

【网易云课堂】  Excel 商业图表修炼秘笈之基础篇