系列文章的第二讲主要关注更加plotnine作图更加细节的地方。首先作为科研工作者，我们很大一部分工作，无论是某个算法相比于之前的算法更加节省时间，稳定性更高，还是生命科学里面一个基因敲除之后蛋白表达水平降低，这些数据结果的呈现方式都需要按照不同的类别分组，然后再进一步呈现，如果用生命科学中的术语来说，我们需要和control进行比较，从而得出我们的结论。这就免不了需要分组做图，这样更加直观，说服力更强。The gramma of graphics这门书，也就是ggplot2的灵感来源，实际上就是告诉我们通过什么样的数据格式去展现不同情况下的数据。
首先，传统的参数统计以及方差分析等作图的方式相对简单，数据量较小，而且使用柱状图和折线图即可，简单的分组和事后比较。然而，非参数统计，或者多元统计方面，则需要一些复杂的方法，比如降维，排序等。下图以python非度量多维尺度为例，详细展示了如何使用sklearn中的MDS函数进行nmds分析并且将数据转化为dataframe，用于plotnine作图。在这个案例中，不同的样品点之间，或者不同的处理之间，可以通通过不同的形状和颜色区别开来，对于嵌套的实验设计同样可以这样做。而且同一个处理之间通过置信椭圆区别开来。这是典型的多元统计的作图展现方法。同样的可以使用PCA等降维方法展现不同处理之间的差异。
# -*- coding: utf-8 -*-
"""@author: Jianshu Zhao, zhaojs16@mails.tsinghua.edu.cnDepartment of Environmental Metagenomics, School of Environment, Tsinghua University"""
import matplotlib
import matplotlib.pylab as pylab
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from biom import load_table
from biom import example_table
from biom import parse_table
import skbio
from sklearn import manifold
from sklearn.decomposition import PCA
from scipy.spatial import distance
from 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 taxonomy
def exploding_panda(_bt):
"""BIOM->Pandas dataframe converterParameters----------_bt : biom.TableBIOM tableReturns-------pandas.DataFrameThe BIOM table converted into a DataFrameobject.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 file
table_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 matrix
DM_dist = distance.squareform(distance.pdist(otu_table_T, metric="braycurtis")) # (n x n) distance measure
DM_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)
# mds
seed = 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_
# nmds
nmds = 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 dataframe
group = 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
或者githubhas2k1/plotninegithub.com
另外一种方法就是可以先在R中计算出nmds的记过再导入到python中可视化，同样的在R中完成可视化也同样很方便。
在完成堆积柱状图的时候，在R中非常方便，只要将对应的数据格式转换，导入读取画图即可，plotnine同样如此，注意将所有的element.text()改成element_text()。
MyColors <- ["#E6478A","#56B4E9", "#E33D36","#7570B0", "#E69F00", "#B3DB70", "#0072B2","#DBDBD9","#944F4A","#5CB5E8","#E39E26","#BFBAD9","#8FD4C7","#DBBA29","#47A133","#945457","#F07D3D","#F5E369","#999999","#ED8CC2"]
#bar
a=(ggplot(data=ups.melt, aes(x=variable, y=value, fill=phylum)) +
geom_bar(stat="identity",width=0.94)+theme_classic()+theme(legend.title = element_blank())+theme(legend_position="right")+scale_fill_manual(values=MyColors)+theme(axis_ticks_x = element_blank(),panel_grid_major=element_blank(),panel_grid_minor=element_blank(),panel_background = element_blank())+theme(legend_text = element_text(size = 12))+theme(legend_key_size = unit(0.5, "cm"))+scale_y_continuous(expand=[0, 0, .1, 0])+labs(y="Relative Abundance",x=""+theme(axis_text_x = element_blank()）相对丰度堆积柱状图
和在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则可能需要添加新的数据分组。
以上主要介绍简单的分组方式。最近太忙了，乘着坐飞机的空隙写的。
扫一扫
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
