查看: 2872|回复: 19

[R语言] R语言给PCA加个圈圈

  [复制链接]
  • TA的每日心情

    6 小时前
  • 签到天数: 94 天

    连续签到: 1 天

    [LV.6]常住居民II

    管理员

    Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

    主题
    349
    奥币
    1079
    积分
    5063
    注册时间
    2018.4.19
    在线时间
    800 小时

    推广达人宣传达人


    发表于 2019.1.9 14:23:10 | 显示全部楼层 |阅读模式
    小伙伴们,在遇到组学实验数据分析得时候,是少不了绘制PCA图的,但是除了常规的PCA图以外,往往也会需要在我们的流程结果的PCA上展现组内样品的分布范围:

    像这样的


    这样的


    其实,本质上这些圆圈是在PCA图的基础上,加了一个置信区间。但是我之前都是在ppt上手动加的。这样总觉得哪里不太对,万一这样和老师说,到时候发文章被审稿人问起怎么办呢。所以这种分析要求,最好的方法是提交给分析。另外我们公司R语言培训中也有PCA的图形绘制。但是如果老师想自己动手,丰衣足食呢,别急,我们有简化版的PCA绘制教程。用R语言绘制完成后,再用PS和AI适当修改一下就可以大功告成了。

    PCA降维分析,不仅可以展现我们的样品重复性,还可以用来发现内在的生物学规律。所以这次我们就围绕这两个主题开展今天帖子的内容。

    一,PCA绘图及美化

    首先今天我们用的R包是ggord(https://github.com/fawda123/ggord
    因为这个包并不放在CRAN中,所以直接用install.packages(“ggord”)是不可行的。
    所以我们需要在github上下载。

    #首先安装devtools
    install.packages('devtools')
    library(devtools)
    #我们再安装ggord
    install_github('fawda123/ggord')
    library(ggord)
    好的,接下来就要进入正题了,以我们的流程中的表达量总表为例(不用预先处理)

    library(ggord)
    library(ggplot2)#我们需要一些ggplot2功能
    getwd()#检查一下路径
    setwd("C:/Users/woney/Desktop")#将工作目录改成桌面
    pca_data=read.table("all.txt",header=T,sep="\t",row=1)#读取表格
    pca_data=t(as.matrix(pca_data))#将表达量总表反转,PCA函数计算的是展现的第一列的ID
    pca_group=factor(c(rep('A',4),rep('B',4), rep('C', '4'), rep('D', '4')))#分组,后续可以自己设置对应的组名,rep('组名','样品数')
          pca1=prcomp(pca_data,center=TRUE,retx=T)#prcomp是R自带的pca分析函数
          p1 <- ggord(pca1, grp_in = pca_group, arrow=0, vec_ext =0,txt=NULL)
    哒哒哒,出图如下
                                 
         
    然后,我们输入?ggord就可以查看具体的参数了。根据帮助文档的参数就可以自行修改。
    我简单列举一下我做的参数修改。

    变成多边形
    p2 <-ggord(pca1, grp_in = pca_group, arrow=0,vec_ext = 0,txt=NULL,ellipse=F,ellipse_pro=0.68,hull=T)

    修改置信区间,以及将椭圆改成空心
    p3<-ggord(pca1, grp_in = pca_group, arrow=0,vec_ext =0,txt=NULL,ellipse=T,ellipse_pro=0.68,poly = F)

    不展示点
    p4 <-ggord(pca1, grp_in = pca_group, arrow=0,vec_ext =0,txt=NULL,ellipse=F,ellipse_pro=0.68,hull=T,obslab =TRUE,repel=TRUE)

    此外,ggord是基于ggplot2基础上绘图,那么理论上gggplot2的图层语法是可以用的。
    比如加上标签
    library(ggplot2)
    p5=p1+geom_text(aes(label=rownames(pca_data)),vjust=1.5,colour="red")

    或者theme()主题都是可以的。总之ggord设置了大量的参数,都可以我们修改图形,其中作者也写了一些例子((https://github.com/fawda123/ggord))或者参考帮助文档中的参数(?ggord),可以帮助我们达成我们的目的。

    二,另类的PCA图

      那么第二点,我们如果想做成有线条和箭头来展现我们数据中的一些变量之间的关系呢。其实呢,在PCA中,最原始的就是会把其中的一些变量给展现出来,但是平时的分析中,我们都不会画出来,大家可以想象一下,当几万的基因代表的线条展现出来的时候。图形长什么样子。通过某个变量所代表的线条在PC1和PC2上的投影,我们可以看出这个变量对样本分离的贡献度,线条越长,代表投影越大,影响越显著,同时两个线条之间的夹角,可以理解为两个变量之间的相关性,夹角小于90度,可以认为两个变量正相关,大于90度,可以认为两个变量负相关。

    那么具体又该如何去画呢,其实用ggord是很简单去实现的。(其实ggord可以画很多图,有CCA,RDA,MAC,LDA等等)
    我们以16s测序为对象,我们的项目中会有一个表格对各个样本的α-多样性指数值数汇总,存于 04.alpha_diversity\summary.alpha_diversity.xls
    如果像展现环境因子与多样性指数之间的关系,可以把环境的理化因子放入其中。如下表


    数据整理完,就可以开始画图了。
    library(ggord)
    library(ggplot2)
    getwd()#检查一下路径
    setwd("C:/Users/woney/Desktop")#将工作目录改成桌面
    pca_data=read.table("summary.alpha_diversity.xls",header=T,sep="\t",row= 1)#读取表格将表达量总表反转,PCA函数计算的是展现的第一列的ID
    pca_group=factor(c(rep('SN','3'),rep('SNP','3'), rep('SNPK','3')))#分组,后续可以自己设置对应的组名
    pca1=prcomp(pca_data,center=TRUE,retx=T)#prcomp是R自带的pca分析函数
    p1 <-ggord(pca1, grp_in = pca_group, arrow=0.3, vec_ext =700,txt=5,repel=T)
    参数也是自己可以修改,结果如下:

    我们简单的可以推测,其中PH和TN影响最大,盐度与α-多样性指数呈负相关,温度呈现正相关。

    本文作者:技术咖

    本帖子中包含更多资源

    您需要 登录 才可以下载或查看,没有帐号?立即注册

    x
    回复

    使用道具 举报

  • TA的每日心情
    yes!
    2019.1.28 16:13
  • 签到天数: 113 天

    连续签到: 1 天

    [LV.6]常住居民II

    钵水母

    Rank: 3Rank: 3

    主题
    6
    奥币
    797
    积分
    163
    注册时间
    2016.5.20
    在线时间
    55 小时

    发表于 2019.1.9 16:45:46 | 显示全部楼层
    好看
    回复

    使用道具 举报

  • TA的每日心情
    害羞
    昨天 11:09
  • 签到天数: 599 天

    连续签到: 2 天

    [LV.9]以坛为家II

    迅猛龙

    Rank: 8Rank: 8

    主题
    26
    奥币
    3589
    积分
    1653
    注册时间
    2016.1.8
    在线时间
    379 小时

    发表于 2019.1.9 17:43:02 | 显示全部楼层
    回复

    使用道具 举报

  • TA的每日心情
    忙~
    6 小时前
  • 签到天数: 231 天

    连续签到: 4 天

    [LV.7]常住居民III

    帝王蝶

    Rank: 4

    主题
    10
    奥币
    971
    积分
    403
    注册时间
    2018.1.3
    在线时间
    158 小时

    热心会员活跃会员


    发表于 2019.1.10 10:02:17 | 显示全部楼层
    以后也会画圈圈啦
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    吃饭
    2019.7.22 10:38
  • 签到天数: 450 天

    连续签到: 1 天

    [LV.9]以坛为家II

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    21
    奥币
    4680
    积分
    1837
    注册时间
    2015.12.29
    在线时间
    427 小时

    突出贡献优秀版主热心会员


    发表于 2019.1.10 18:22:31 | 显示全部楼层
    加多一个标题
    p6=p5+abs(title = "PCA Figure")+ theme(plot.title = element_text(hjust = 0.5))

    评分

    参与人数 1奥币 +5 收起 理由
    基迪奥-李泽标 + 5 楼主V5!

    查看全部评分

    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    害羞
    2019.8.7 09:16
  • 签到天数: 312 天

    连续签到: 1 天

    [LV.8]以坛为家I

    帝王蝶

    Rank: 4

    主题
    2
    奥币
    789
    积分
    450
    注册时间
    2016.3.1
    在线时间
    113 小时

    发表于 2019.1.11 13:05:00 | 显示全部楼层
    回复

    使用道具 举报

  • TA的每日心情

    2019.1.15 08:44
  • 签到天数: 385 天

    连续签到: 1 天

    [LV.9]以坛为家II

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    3
    奥币
    2564
    积分
    544
    注册时间
    2016.3.23
    在线时间
    230 小时

    突出贡献优秀版主


    发表于 2019.1.11 14:37:30 | 显示全部楼层
    楼 主。来电数据练练手撒
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    吃饭
    2019.8.6 11:03
  • 签到天数: 477 天

    连续签到: 3 天

    [LV.9]以坛为家II

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    62
    奥币
    8320
    积分
    3223
    注册时间
    2017.9.21
    在线时间
    305 小时

    突出贡献优秀版主论坛元老


    发表于 2019.1.11 15:14:37 | 显示全部楼层
    优秀~
    回复

    使用道具 举报

  • TA的每日心情

    6 小时前
  • 签到天数: 94 天

    连续签到: 1 天

    [LV.6]常住居民II

    管理员

    Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

    主题
    349
    奥币
    1079
    积分
    5063
    注册时间
    2018.4.19
    在线时间
    800 小时

    推广达人宣传达人


     楼主| 发表于 2019.1.11 15:28:30 | 显示全部楼层
    gshwq0001 发表于 2019.1.11 14:37
    楼 主。来电数据练练手撒

    数据用老师的,所以不能提供给大家。大家参考一下思路就可以啦
    回复 支持 反对

    使用道具 举报

    该用户从未签到

    草履虫

    Rank: 2

    主题
    0
    奥币
    30
    积分
    5
    注册时间
    2018.11.28
    在线时间
    3 小时

    发表于 2019.1.11 22:50:01 | 显示全部楼层
    包没装上,尴尬
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    忙~
    6 天前
  • 签到天数: 70 天

    连续签到: 1 天

    [LV.6]常住居民II

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    175
    积分
    67
    注册时间
    2018.2.22
    在线时间
    23 小时

    发表于 2019.1.13 19:04:57 | 显示全部楼层
    回复

    使用道具 举报

  • TA的每日心情
    忙~
    2019.3.13 16:11
  • 签到天数: 3 天

    连续签到: 1 天

    [LV.2]偶尔看看I

    草履虫

    Rank: 2

    主题
    0
    奥币
    50
    积分
    22
    注册时间
    2019.1.11
    在线时间
    4 小时

    发表于 2019.1.14 10:50:01 | 显示全部楼层
    优秀的人
    回复

    使用道具 举报

  • TA的每日心情
    害羞
    2019.7.15 11:34
  • 签到天数: 15 天

    连续签到: 1 天

    [LV.4]偶尔看看III

    中华鲟

    Rank: 5Rank: 5

    主题
    43
    奥币
    2011
    积分
    660
    注册时间
    2018.7.19
    在线时间
    53 小时

    热心会员活跃会员


    发表于 2019.1.14 17:06:22 | 显示全部楼层
    学习了
    回复

    使用道具 举报

  • TA的每日心情

    2018.12.10 10:36
  • 签到天数: 1 天

    连续签到: 1 天

    [LV.1]初来乍到

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    54
    积分
    73
    注册时间
    2016.12.11
    在线时间
    28 小时

    发表于 2019.1.16 05:19:04 | 显示全部楼层
    谢谢
    回复

    使用道具 举报

  • TA的每日心情
    害羞
    6 天前
  • 签到天数: 31 天

    连续签到: 1 天

    [LV.5]常住居民I

    钵水母

    Rank: 3Rank: 3

    主题
    2
    奥币
    105
    积分
    61
    注册时间
    2018.1.10
    在线时间
    15 小时

    发表于 2019.3.21 08:48:00 | 显示全部楼层
    请问它的group是自己设置数量吗?而不是自己根据聚类情况而自动分配的?
    回复 支持 反对

    使用道具 举报

  • TA的每日心情

    2 小时前
  • 签到天数: 215 天

    连续签到: 1 天

    [LV.7]常住居民III

    帝王蝶

    Rank: 4

    主题
    0
    奥币
    672
    积分
    210
    注册时间
    2018.8.28
    在线时间
    53 小时

    发表于 2019.3.22 14:20:57 | 显示全部楼层
    收藏
    回复

    使用道具 举报

  • TA的每日心情
    吃饭
    2019.7.29 11:42
  • 签到天数: 43 天

    连续签到: 1 天

    [LV.5]常住居民I

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    185
    积分
    88
    注册时间
    2017.10.22
    在线时间
    30 小时

    发表于 2019.3.28 20:02:24 | 显示全部楼层
    非常好
    回复

    使用道具 举报

  • TA的每日心情
    害羞
    2019.6.26 11:16
  • 签到天数: 26 天

    连续签到: 1 天

    [LV.4]偶尔看看III

    草履虫

    Rank: 2

    主题
    3
    奥币
    110
    积分
    46
    注册时间
    2018.4.10
    在线时间
    16 小时

    发表于 2019.5.15 16:55:32 | 显示全部楼层
    本帖最后由 ayi 于 2019.5.15 19:13 编辑

    请问怎么更改显示符号啊,不用圆点改成三角形,谢谢
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    好棒
    8 小时前
  • 签到天数: 602 天

    连续签到: 2 天

    [LV.9]以坛为家II

    帝王蝶

    Rank: 4

    主题
    2
    奥币
    1371
    积分
    320
    注册时间
    2016.12.1
    在线时间
    187 小时

    发表于 2019.5.28 16:34:52 | 显示全部楼层
    pca_data=read.table("summary.alpha_diversity.xls",header=T,sep="\t",row= 1)  这一行代码中,header=T和sep="\t"代表什么意思?  
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    好棒
    8 小时前
  • 签到天数: 602 天

    连续签到: 2 天

    [LV.9]以坛为家II

    帝王蝶

    Rank: 4

    主题
    2
    奥币
    1371
    积分
    320
    注册时间
    2016.12.1
    在线时间
    187 小时

    发表于 2019.5.28 16:37:16 | 显示全部楼层
    p1 <-ggord(pca1, grp_in = pca_group, arrow=0.3, vec_ext =700,txt=5,repel=T)  这行代码中grp_in = pca_group, arrow=0.3, vec_ext =700,txt=5,repel=T代表什么意思?
    回复 支持 反对

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    快速回复 返回顶部 返回列表