查看: 838|回复: 1

[R语言] 如何用R来计算离散系数?

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

    2019.9.23 10:01
  • 签到天数: 101 天

    连续签到: 1 天

    [LV.6]常住居民II

    管理员

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

    主题
    392
    奥币
    1106
    积分
    5611
    注册时间
    2018.4.19
    在线时间
    893 小时

    推广达人宣传达人


    发表于 2019.8.8 09:13:28 | 显示全部楼层 |阅读模式
    在一些实验设计中,大家比较关注的是在不同品系间的内参基因。类比于18s和actin,内参基因首先的定义就是不受内外因素影响、表达量趋向恒定的一类基因。虽然这些基因一般都是一些非常重要以及“基础”的管家基因。
    一般做qPCR等实验,用18s、GAPDH、actin这些内参基因就足够了。但一般转录组数据因为其通量高,一下检测几万的基因表达量,如此丰富的表达量信息,就为根据表达量的情况来筛选内参基因提供了一个很好的基础。

    言归正传,如何根据基因的表达量来进行内参基因的筛选呢?

    目前,我的常用方法是计算Unigene的离散系数(CV)。简单讲就是CV=标准差/平均数,数值越小,这个基因在不同样本间的表达量的变化越小。

    这个定义可以看:

    百度百科https://baike.baidu.com/item/%E5%8F%98%E5%BC%82%E7%B3%BB%E6%95%B0
    维基百科:https://en.wikipedia.org/wiki/Coefficient_of_variation

    当然,用这个CV来计算,并得到内参基因是否有文献支持,说实话,因为我没有这方面的研究背景,所以并不清楚。只是从统计学上讲,CV是非常不错的筛选标准。

    如果大家有更好的统计参数来衡量和筛选,非常欢迎留言讨论!目前看, CV值是个很不错的标准, 甚至有人为此写了一个R包。大家感兴趣的话可以去看一下这篇文章。



    用R来计算CV

    今天的主要内容是如何利用R来计算CV。以前,我会选择用excel来进行计算CV。

    其实excel计算起来并不复杂,用函数STDEV()和AVERAGE()就可以了,只需要计算两个值的比值就可以。
    但是为什么用R来计算,因为R首先是专业的统计软件,(并不是绘图软件!!!),并且用R会显得有逼格嘛。
    其实主要是写一个脚本,以后遇到这样的分析,就可以自动跑出来了啦,这是最重要的。
    为什么大家会认为R是绘图软件,我觉得这锅要Haddly来背(ggplot2的作者)。

    开始用R之前

    因为,我的R主要是靠自学和看书,我看的第一本是Haddly的《ggplot2:图形绘画与艺术》,第二本是Haddly的《R for data science》,目前第三本估计还是Haddly的《advance R》。

    兄弟,发现问题了没。所以我的R命令很多都依赖Haddly的R包。但其实Haddly只是众多R开发者之一,虽然无疑Haddly是改变R的人,但是R是一个社区,很多其他人R包也是很优秀的。

    所以我的命令有的时候就很不协调和冗余,这是因为为了兼顾Haddly的包特有的格式。比如tibble和R自带的data.frame之间的行名之间的冲突。、
    今天我才稍微知道了colunm_to_rownames命令如何协调两者。

    所以各位看官,看我的命令一定要知道为什么这样,以及换成简单命令。或者你也学学Haddly的R包吧。

    命令行


    在开始之前,推荐大家一个IDE和文本编辑器。

    sublime文本编辑器

    在此之前,有位朋友向我推荐了sublime这个文本编辑器,这个编辑器真心牛逼,针对R的脚本,实现了非常美观的代码高亮,非常nice。



    另外,在安装一些插件后,就可以变成R的IDE,这个非常完美了。

    我既嫌弃Rgui的简陋界面,又嫌弃Rstudio的繁琐界面。终于在安装sublime后,既有Rgui的简洁界面,又有Rstudio的代码补全。

    用R计算离散系数

    注释看代码行:

    [AppleScript] 纯文本查看 复制代码
    # R 3.5.2 win 10
    require(tidyverse)
    require(matrixStats)#这个包有rowSds函数,计算一行的sd标准差
    path = "C:/Users/woney/Desktop" # 这是我的工作目录桌面
    inname = "all.sample.annot.xls"
    setwd(path)
    dat=read_tsv(inname, col_names = T)
    nms <- str_subset(colnames(dat),pattern = "pm|km")
    dim(dat)
    dat_m <- dat %>% filter(rowSums(.[nms]) != 0) %>%
    mutate(mean = rowMeans(.[nms]),sd = rowSds(as.matrix(.[nms]))) %>%
    mutate(CV = sd/mean)
    dat_cv <- dat_m
    write_tsv(dat_cv,"unigene_cv.xls",col_names = T)
    dat_quan <- quantile(dat_cv[["CV"]],prob = 0.01)#寻找排名前1%的最小CV值基因
    dat_cv_fit <- dat_cv %>% arrange(CV) %>% filter(mean >= 100, CV <= dat_quan)
    write_tsv(dat_cv_fit,"unigene_cv_filter.xls",col_names = T)


    到此,我们就得到了一个完整结果的文件,一个我们筛选后的表格,但是感觉还是缺一个图形。

    绘制图形

    对于这个目标结果的可视化,我没什么好的经验。所以暂且用散点图来绘制,但是感觉不太匹配。要是有好的建议,可以分享一下,一起学习和进步~

    [AppleScript] 纯文本查看 复制代码
    #绘制图形
    p1 <- ggplot(dat =dat_cv_fit, mapping = aes(x = mean, y = sd))+
    geom_point(dat = dat_cv,aes(x = mean, y = sd),fill = "skyblue",size = 0.5)+
    geom_point(dat = dat_cv_fit,aes(fill = id),shape = 21,color ="red",size = 3)+
    labs(x = "The mean of Unigene RPKM",y = "The stdev of Unigene RPKM" ,title = "The CV value")+
    geom_vline(xintercept = c(100,2000),linetype = 4,color = "grey")
    #进行一点调整
    p2 <- p1 + theme_bw()+ theme(legend.position = "none")
    





    检验结果

    毕竟是自己的脚本,R实现的优势是速度快,但是只要其中一个命令出错,就可能导致结果的差异。对于自己不熟悉的东西,尤其是计算类的,还是小心一点比较好。

    所以我在excel上检查了一下,我的结果还是一致的,可靠!




    今天的内容就到这里啦,如果大家有什么问题或者建议,欢迎留言讨论。
         

    本帖子中包含更多资源

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

    x
    回复

    使用道具 举报

  • TA的每日心情
    害羞
    2019.9.24 10:09
  • 签到天数: 33 天

    连续签到: 1 天

    [LV.5]常住居民I

    迅猛龙

    Rank: 8Rank: 8

    主题
    161
    奥币
    1930
    积分
    1895
    注册时间
    2015.11.14
    在线时间
    76 小时

    发表于 2019.8.8 11:13:48 | 显示全部楼层
    学习了
    我的意中人是一个盖世英雄,有一天他会踩着七彩祥云来娶我……
    回复

    使用道具 举报

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

    本版积分规则

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