10X单细胞(10X空间转录组)Seurat分析之QQplot的详细解释及绘制

下面这张图大家应该都不陌生吧

这个图是函数JackStrawPlot()画出来的,用于我们分析全多少个主成分(PC)用于下游分析,大家应该都知道,那么我们来看看这个图的解释。

Plots the results of the JackStraw analysis for PCA significance. For each PC, plots a QQ-plot comparing the distribution of p-values for all genes across each PC, compared with a uniform distribution. Also determines a p-value for the overall significance of each PC .

Significant PCs should show a p-value distribution (black curve) that is strongly skewed to the left compared to the null distribution (dashed line) The p-value for each PC is based on a proportion test comparing the number of genes with a p-value below a particular threshold (score.thresh), compared with the proportion of genes expected under a uniform distribution of p-values.

这是官方的解释,这里有一个很有意思的地方,QQ-plot,什么是QQ-plot????是用来做什么的???我们来详细看一下

QQ plot也就是Quantile-Quantile Plots。是通过比较两个概率分布的分位数对这两个概率分布进行比较的概率图方法。

其想法就是,如果现在有从某个类型的概率分布中抽取的N个数据,那么如果想确定这个概率分布是否接近normal distribution该怎么办呢?

一种做法就是先假设这些点是符合normal distribution的,那么现在取一个正态分布作为对比模板,将其pdf下的面积等分成N份,每一块面积代表的概率都是相等的,都是1/N。那么如果现在我们要从这个正态分布中抽一个随机数,理想情况下,这个数出现在其中任何一个面积相等的小区域的概率都是相等的,因为这些小区域代表的发生的概率相同。也就是说,如果该正态分布的概率分布范围被分为N个等份,而我们刚好要抽N个数出来,那么理想状况下,应该是刚好这N个小区域中都被抽出了一个数,并且我们会预期这些数的值应该是每个区域的中位数(或二分位数)。

这些中位数(或二分位数)的值分别应该是多少呢?应该是下面这个式子:

n=N, i代表是第几个概率区间的中位数, -1代表cdf的反函数,括号里填入你想找有x%概率小于某个数y,那么这个函数会给出这个数y是多少。例如,对于standard normal distribution,Φ-1(0.1) =-1.28就代表抽出来的数有10%概率小于-1.28. 而(i-0.5)/n就是(0,1/N), (1/N,2/N)...((N-1)/N,1)这些概率区间的中位点(或二分位点)。

举个例子,假如我们有5个样本点,N=5,要判断它们是否符合standard normal, 那么如上面所说,我们要从作为对比模板的standard normal中抽5个数出来,也就是要把standard normal的概率分布范围分成5等份.可以看到下图一共有5个点,其中有4个分位点,这些点左边的面积分别是0.2,0.4,0.6,0.8,或者说从这个standard normal抽出来的随机数小于这些点的概率分别是20%,40%,60%,80%。定义上最后的第5个点不能叫分位点,但是它表示的含义还是和上面一样,即小于这个点的概率是100%. 需要注意的是,0.2,0.4...这些数值并不是真正的这些点的x值,而是代表小于这些点的概率。所以,我们也可以从图上看出来0.2与0.4之间的距离和0.4与0.6之间的距离明显不一样。这些点的实际x值如下图中第三行数字所示。

我们看到,分成的5等份概率区间分别为(0,1/5),(1/5,2/5),(2/5,3/5))...(4/5,1),即(0,0.2),(0.2,0.4)...(0.8,1)。此时要找出这些概率区间的中位数(或二分位数),也就是分别将这些概率区间的概率二等分的点,如上图中红点所示。我们发现我们要找的点就是Φ-1(0.1) ,Φ-1(0.3) ,Φ-1(0.5) ,Φ-1(0.7) ,Φ-1(0.9) ,即随机数分别有0.1,0.3...0.9概率小于该点的该点x值,它们分别为-1.282,-0.524...1.282。而公式 (i-0.5)/n就是帮我们总结出了Φ-1()括号里的数。Φ-1(0.1)...Φ-1(0.9)这些值就应该是我们从这个standard normal distribution中抽取5个随机数时,理想状态下,预期抽到的值。即每个概率区间都恰好被抽到一个数,而被抽到的数恰好是这些概率区间的中位数。将实际抽到的值(也就是需要判断概率分布是什么样子的样本点的值)从小到大排列好,如果这些样本点符合standard normal的话,它们的值在理想状态下应该就分别等于前面求出的那些期望值,即-1.282,-0.524...1.282。如果将这些实际抽到的值作为纵轴,理论上期望抽到的值作为点的横轴,画出这些点,这就叫做QQplot。即我们认为理想状态下,从某个概率分布中抽随机数,这些随机数就应该是这个概率分布的分位数,那么比较样本概率分布的分位数和作为对比模板的standard normal distribution的分位数,看它们的值是否相似,就可以判断样本点的概率分布是不是接近standard normal distribution了。由于理想状态下y=x,所以这些点应该落在y=x这样一条过0点,斜率45度的直线上。如果实际的样本点的值符合的是一般的normal distribution,mean=μ,std= , QQplot上画出的点就应该落在y= x+μ这条直线上,因为此时每个实际的样本点的值,应该就是理论期望值经过 x+μ这样的线性变化得到的。如果这些样本点并不符合normal distribution,那么画出来的QQplot中的点就不呈一条直线分布

另外,我们会发现如果样本点确实取自一个normal distribution, 那么样本点数量越大,QQplot就越接近一条直线,反之,当样本点数量过少时,QQplot中的点很容易偏离直线,类似大数定律的原理。

补充一下,Normal distribution中常说的的几倍 也是分位数,当它为standard normal时,μ=0, =1,这些分位数的值就刚好是1,2,3...而对于一般的normal distribution,这些分位数就为μ+ x,即μ+ ,μ+2 ,μ+3 ...

需要注意的是,在找前面所说的,抽取N个随机数应该期望得到的值的这一步,其实有很多不同版本的做法。有像前面所说,将作为对比模板的概率分布分成面积相等的N等份,然后取每个等份的中位数作为期望值的。还有的是将模板的概率分布分成面积相等的N+1等份,这样就有N个分位点,然后直接取这些分位数作为抽取N个随机数期望得到的值。可以看出,后面这种方法会和前一种方法算出的期望的值有些许不同,但是其实这对概率分布类型的判断并没有什么影响,因为我们只在意最后画出的点是否在一条直线上,这些期望的值的不同取法,只会使得直线上下平移而已,就像一般的normal distribution和standard normal distribution的分位点的值只是有整体μ的平移一样,这些分位点的间距是不变的

QQplot可以用来判断样本点是否符合任意一种概率分布,其方法就是和前面所讲normal distribution的例子一样,如果怀疑样本点符合某种分布,那么就找出从这样一个分布中抽取相同数目随机数时期望得到的数值,将这些数值与实际的样本点的值放在QQplot中进行比较,如果点的分布接近一条直线就说明,样本点的分布和猜测的分布相同,反之则不同。

这部分理论真的是有点费脑子

我们接下来实现一下这个图(R+Python)

我们还是简单理解一下QQ图吧,QQ图是用于验证一组数据是否符合正态分布,或者验证某两组数据是否来自同一分布情况,是一种散点图,通常情况下,其横坐标为标准正态分布的分位数,纵坐标为样本值。要利用QQ图判定测试样本数据是否近似于正态分布,只需看QQ图上的点是否近似地在一条直线附近,更多关于QQ图的含义理解,小伙伴们可自行搜索哈~~。QQ图样例如下:

QQ图(Quantile-Quantile Plots)绘制方法(R版本,ggplot2)

library(tidyverse)
library(ggtext)
library(hrbrthemes)
library(wesanderson)
library(LaCroixColoR)
library(RColorBrewer)
library(qqplotr)
library(ggsci)
#构造数据
df <- data.frame(y = rt(200, df = 5))
#可视化绘制
ggplot(df, aes(sample=y))+
  stat_qq(shape=21,fill="gray",size=4,colour="black",stroke=.5) +
  stat_qq_line(color="red") +
  labs(x = "Theoretical Quantiles", y = "Sample Quantiles",
    title = "Example of <span style='color:#D20F26'>ggplot2::stat_qq function</span>",
    subtitle = "processed charts with <span style='color:#1A73E8'>stat_qq()</span>",
    caption = "Visualization by <span style='color:#0057FF'>DataCharm</span>") +
  hrbrthemes::theme_ipsum(base_family = "Roboto Condensed") +
  theme(
    plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
                                  size = 20, margin = margin(t = 1, b = 12)),
    plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
    plot.caption = element_markdown(face = 'bold',size = 12)
      )

QQ图(Quantile-Quantile Plots)绘制方法(R版本,qqplotr包)

library(qqplotr)
set.seed(0)
smp <- data.frame(norm = rnorm(100))
ggplot(data = smp, mapping = aes(sample = norm)) +
    stat_qq_band(alpha=.4) +
    stat_qq_point(shape=21,size=4,fill="gray60",colour="black",stroke=.5) +
    stat_qq_line(color="red") +
    labs(x = "Theoretical Quantiles", y = "Sample Quantiles",
         title = "Example of <span style='color:#D20F26'>qqplotr::stat_qq_point__* function</span>",
         subtitle = "processed charts with <span style='color:#1A73E8'>stat_qq_point__*()</span>",
         caption = "Visualization by <span style='color:#0057FF'>DataCharm</span>"
        ) +
    hrbrthemes::theme_ipsum(base_family = "Roboto Condensed") +
    theme(
    plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
                                  size = 20, margin = margin(t = 1, b = 12)),
    plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
    plot.caption = element_markdown(face = 'bold',size = 12)
      )

注意:qqplotr包提供了stat_qq_band() 函数用于绘制区间,使QQ图更加易于理解。

多类别,R版本

ggplot(data = smp, mapping = aes(sample = norm)) +
    geom_qq_band(bandType = "ks", mapping = aes(fill = "KS"), alpha = 0.5) +
    geom_qq_band(bandType = "ts", mapping = aes(fill = "TS"), alpha = 0.5) +
    geom_qq_band(bandType = "pointwise", mapping = aes(fill = "Normal"), alpha = 0.3) +
    geom_qq_band(bandType = "boot", mapping = aes(fill = "Bootstrap"), alpha = 0.5) +
    stat_qq_point(shape=21,size=4,fill="gray60",colour="black",stroke=.5) +
    stat_qq_line() +
    labs(x = "Theoretical Quantiles", y = "Sample Quantiles",
         title = "Example of <span style='color:#D20F26'>qqplotr::stat__qq_point__* function</span>",
         subtitle = "processed charts with <span style='color:#1A73E8'>stat__qq_point__*()</span>",
         caption = "Visualization by <span style='color:#0057FF'>DataCharm</span>") +
    ggsci::scale_fill_jco(name="Bandtype") +
    hrbrthemes::theme_ipsum(base_family = "Roboto Condensed") +
    theme(
        plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
                                      size = 20, margin = margin(t = 1, b = 12)),
        plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
        plot.caption = element_markdown(face = 'bold',size = 12),
      )

分面绘制,R版本

data("barley", package = "lattice")
ggplot(data = barley, mapping = aes(sample = yield, color = site, fill = site)) +
    stat_qq_band(alpha=0.5) +
    stat_qq_line() +
    stat_qq_point() +
    ggsci::scale_color_nejm() +
    ggsci::scale_fill_nejm() +
    facet_wrap(~ site) +
    labs(x = "Theoretical Quantiles", y = "Sample Quantiles", 
         title = "Example of <span style='color:#D20F26'>qqplotr::stat_qq_point__* function</span>",
         subtitle = "processed charts with <span style='color:#1A73E8'>stat__qq_point__*()</span>",
         caption = "Visualization by <span style='color:#0057FF'>DataCharm</span>") +
    hrbrthemes::theme_ipsum(base_family = "Roboto Condensed") +
    theme(
        plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
                                      size = 20, margin = margin(t = 1, b = 12)),
        plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
        plot.caption = element_markdown(face = 'bold',size = 12),
      )

更多详细例子和绘图参数,小伙伴们可去qqplotr官网进行查询。

QQ图Python绘制教程

Python 绘制QQ图主要借助其用于统计分析的statsmodels库和scipy库,样例如下:

statsmodels范例:

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = "Times New Roman"

# 构建数据
test = np.random.normal(0,1, 200)
# 开始绘图
fig,ax = plt.subplots(figsize=(6,5),dpi=100,)
ax.set_facecolor("white")
sm.qqplot(test, line='45',marker='.', markerfacecolor='gray', markeredgecolor='k',markeredgewidth=.5,
          markersize=14,ax=ax)
ax.set_xlim((-3, 3.0))
ax.set_ylim((-3, 3.0))
ax.grid(linestyle="--",linewidth=.5,color="black")
ax.tick_params(direction='in',labelsize=13)
titlefontdict = {"size":17,"color":"k",'family':'Times New Roman'}
ax.set_title('Example Of QQ Plot In Python-Statsmodels Make',titlefontdict,pad=20)
ax.text(.82,.056,'\nVisualization by DataCharm',transform = ax.transAxes,
        ha='center', va='center',fontsize = 10,color='black')

注意:sm.qqplot()在绘制QQ图时,其定制化绘制的灵活性较大,可以对散点样式、颜色、大小、粗细等属性进行设置。更多内容可参考:gofplots.qqplot()

scipy版本

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
fig,ax = plt.subplots(figsize=(6,5),dpi=100)
stats.probplot(test, plot=ax)
ax.set_xlim((-3, 3.0))
ax.set_ylim((-3, 3.0))
ax.grid(linestyle="--",linewidth=.5,color="black")
ax.tick_params(direction='in',labelsize=13)
titlefontdict = {"size":17,"color":"k",'family':'Times New Roman'}
ax.set_title('Example Of QQ Plot In Python-scipy Made',titlefontdict,pad=20)
ax.text(.82,.056,'\nVisualization by DataCharm',transform = ax.transAxes,
        ha='center', va='center',fontsize = 10,color='black')
plt.savefig(r'G:\DataCharm\可视化包介绍(绘制)\QQ图绘制\QQ charts in Python-scipy.png',width=6,height=5.5,
            dpi=900,bbox_inches='tight',facecolor='white')

可以看出:stats.probplot()绘制QQ图,其在定制化操作上有所不足,无法较灵活的设置颜色、边框宽度等属性。更多详细内容可参考:scipy.stats.probplot()

不知道大家喜欢那个语言绘图,个人非常喜欢python绘图,灵活度很高

一些资源

* R-qqplotr包官网: https://aloy.github.io/qqplotr/
* statsmodels.graphics.gofplots.qqplot: https://www.statsmodels.org/stable/generated/statsmodels.graphics.gofplots.qqplot.html
* scipy.stats.probplot(): https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.probplot.html

前路多艰,生活很好,等你超越~~~

  • 27
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值