R统计练习题

基础统计

1、求众数（mode）：一组数据中出现次数最多的数值，叫众数，有时众数在一组数中有好几个。用M表示。 理性理解：简单的说，就是一组数据中占比例最多的那个数。
例如：2，3，3，3，4，5的众数是3。 但是，如果有两个或两个以上个数出现次数都是最多的，那么这几个数都是这组数据的众数。 例如：2，2，3，3，4，5的众数是2和3。 其次，如果所有数据出现的次数都一样，那么这组数据没有众数。 例如：2，3，4，5没有众数。
众数在R中没有特定的函数，需自己建立函数：（原文出自【易百教程】，原文链接：https://www.yiibai.com/r/r_mean_median_mode.html）

getmode <- function(v) {
   uniqv <- unique(v)
   uniqv[which.max(tabulate(match(v, uniqv)))]
}      #自定义函数
colnames(dat)
mode_Sepal.Length <- getmode(dat$Sepal.Length);mode_Sepal.Length
mod_Sepal.Width <- getmode(dat$Sepal.Width);mod_Sepal.Width
mod_Petal.Length <- getmode(dat$Petal.Length);mod_Petal.Length
mod_Petal.Width <- getmode(dat$Petal.Width);mod_Petal.Width

2、计算数据集 iris的前两列变量的相关性，提示cor函数可以选择3种methods

cor(dat$Sepal.Length, dat$Sepal.Width, use = "everything",
    method = "pearson")
[1] -0.1175698

3、对数据集 iris的所有定量数据列内部zcore标准化，并计算标准化后每列的平均值和标准差
（参考：1、https://zhuanlan.zhihu.com/p/33727799；2、https://www.cnblogs.com/ljhdo/archive/2018/08/09/4899086.html）

zcore标准化

主要是scale(）函数：

scale(x, center = TRUE, scale = TRUE)

参数注释：x是数值矩阵，当center=TRUE时，为数据对象x按列进行中心化；当scale=TRUE时，为数据对象x按列进行标准化。
中心化是把每个数据减去均值；
标准化是在中心化后的数据基础上再除以数据的标准差；
默认情况下，center=TRUE，scale=TRUE，scale()函数首先把一组数的每个数都减去这组数的平均值，然后除以这组数的均方根。
如果scale=TRUE，而center=FALSE，那么，scale()函数不会把一组数中的每个数减去平均值，而直接处以这组数据的均方根。

#method1
scale(data,center=T,scale=F) #center表示数据中心化，scale表示数据标准化
scale(data,center=T,scale=T) 或scale(data)
a<-round(scale(data),2)  #约取2位有效数字
a
        V1    V2    V3    V4
[1,] -1.16 -1.16 -1.16 -1.16
[2,] -0.39 -0.39 -0.39 -0.39
[3,]  0.39  0.39  0.39  0.39
[4,]  1.16  1.16  1.16  1.16
#看一下正态分布图片
d<-density(a)
plot(d)
#method2
#或者调用BBmisc包中的normalize函数
normalize(data) #be equal to scale(data)
#method3
#调用limma包中normalizeBetweenArrays()函数
library(limma) 
exprSet=normalizeBetweenArrays(exprSet)  ##表达矩阵的标准化

4、计算列内部zcore标准化后 iris的前两列变量的相关性

> dat <- as.data.frame(dat)
> attach(dat)
> cor(Sepal.Length, Sepal.Width, use = "everything",
+     method = "pearson")
[1] -0.1175698
> detach(dat)

5、载入r包airway并且通过assay函数拿到其表达矩阵后计算每列之间的相关性

options(stringsAsFactors = F)
library(airway)
data(airway)
# 这里需要自行学习bioconductor里面的RangedSummarizedExperiment对象
airway 
RNAseq_expr=assay(airway)
colnames(RNAseq_expr) 
RNAseq_expr[1:4,1:4] 
# RNAseq_expr 是一个数值型矩阵，属于连续性变量，可以探索众数、分位数和平均数 ，极差，方差和标准差等统计学指标
RNAseq_expr <- log2(RNAseq_expr+1)  #进行log2()转化
RNAseq_gl=colData(airway)[,3]
table(RNAseq_gl)
RNAseq_gl <- relevel(RNAseq_gl,ref = "untrt")
annotation_col=data.frame(group=RNAseq_gl)
rownames(annotation_col)=colnames(RNAseq_expr) 
M <- cor(RNAseq_expr);M
pheatmap::pheatmap(M,show_rownames = T, show_colnames = T, annotation_col = annotation_col)
相关性热图

表达矩阵相关

Q1: 把RNAseq_expr第一列全部加1后取log2后计算平均值和标准差

tmp=log2(RNAseq_expr[,1]+1)
mean(tmp)
sd(tmp)

Q2: 根据上一步得到平均值和标准差生成同样个数的随机的正态分布数值

a=rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
a=sort(a)
plot(a)
points(sort(tmp))

Q3: 删除RNAseq_expr第一列低于5的数据后，重复Q1和Q2

tmp=RNAseq_expr[,1]
tmp=tmp[tmp>5]
tmp=log2(tmp)
a=rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
a=sort(a)
plot(a)
points(sort(tmp))

Q4: 基于Q3对RNAseq_expr的第一列和第二列进行T检验

x=RNAseq_expr[,1]
x=x[x>5]
x=log2(x)   #视情况进行转化

y=RNAseq_expr[,2]
y=y[y>5]
y=log2(y) #视情况进行转化

t.test(x,y)
library(ggpubr)
df=data.frame(value=c(x,y),
 group=c(rep('x',length(x)),rep('y',length(y))))
ggboxplot(df, y = "value", x = "group")

Q5: 取RNAseq_expr行之和最大的那一行根据分组矩阵进行T检验

#法一：（没有标准化）
pos=which.max(rowSums(RNAseq_expr))
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos
#法二：（log2转化后的数据）
df <- as.data.frame(RNAseq_expr)
df$sum <- apply(df,1,sum)
df <- df[order(df$sum,decreasing = T),]
dat <- df[1,1:ncol(df)-1]
library(reshape2)
dat_m <- melt(dat)
group <- colData(airway)[,3]
group <- relevel(group,ref = "untrt");group
dat_m$group <- group
t.test(dat_m$value~dat_m$group)

Q6: 取RNAseq_expr的MAD最大的那一行根据分组矩阵进行T检验

pos=which.max(apply(RNAseq_expr,1,mad))
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos

Q7: 对RNAseq_expr全部加1后取log2后重复Q5和Q6,看看是不是基因变化了，统计结果也变化了

RNAseq_expr=log2(RNAseq_expr+1)
pos=which.max(rowSums(RNAseq_expr))
pos
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos=which.max(apply(RNAseq_expr,1,mad))
pos
t.test(RNAseq_expr[pos,]~RNAseq_gl)

Q8: 取RNAseq_expr矩阵的MAD最高的100行，对列和行分别进行层次聚类（画聚类树状图）

cg=names(tail(sort(apply(RNAseq_expr,1,mad)),100))
dat=RNAseq_expr[cg,]
plot(hclust(dist(t(dat))))   #主要对样本进行聚类
colnames(dat)
RNAseq_gl
plot(hclust(dist( dat )))
#检查一下聚类结果跟样本的处理信息是否能对应
mad100的样本聚类（偶数为对照，奇数为处理组）

Q9: 取RNAseq_expr矩阵的SD最高的100行，对列和行分别进行层次聚类（画聚类树状图）

dat <- tail(sort(apply(RNAseq_expr,1,sd)),100)
df <- RNAseq_expr[names(dat),]
plot(hclust(dist(t(df))))
colnames(df)
group
plot(hclust(dist(df)))
#不同的标准（mad、sd）得到的聚类结果不一样，再从主成分分析方面进行验证
dat=as.data.frame(t(RNAseq_expr))
dat[1:4,1:4]
library(FactoMineR)#画主成分分析图需要加载这两个包
library(factoextra) 

#pca的统一操作走起
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = group, # color by groups
                         #palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups"
)
pca_plot  #浓度椭圆

pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = group, # color by groups
                         #palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, 
                         ellipse.type = "confidence",
                         legend.title = "Groups"
)
pca_plot  #置信椭圆

  • 感想不同的标准（mad、sd）得到的聚类结果不一样。


    sd100的样本聚类（偶数为对照，奇数为处理组）

    PCA主成分分析

Q10: 对Q8矩阵按照行和列分别归一化并且热图可视化

cg=names(tail(sort(apply(RNAseq_expr,1,mad)),100))
dat=RNAseq_expr[cg,]
pheatmap::pheatmap(scale(dat))  #scale默认对列进行标准化
pheatmap::pheatmap(t(scale(t(dat))))
#使用mad前100个基因并不能画出有明显差异的热图，还需要进行差异分析，效果可能更好。

统计检验相关

这里需要对前面的RNAseq_expr矩阵进行一定程度的过滤，主要是过滤那些每一列都为0的行

options(stringsAsFactors = F)
rm(list=ls())
library(airway) 
data("airway")
RNAseq_expr=assay(airway)
e1=RNAseq_expr[apply(RNAseq_expr,1,function(x) sum(x>0)>1),] #过滤掉每列均为0的行
colnames(RNAseq_expr) 
RNAseq_gl=colData(airway)[,3]
table(RNAseq_gl)

Q1: 对e1每一行独立根据分组矩阵进行T检验，检查为什么有些行T检验失败

apply(e1, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
})

Q2: 找出T检验失败的行并且从e1矩阵剔除掉

e1_a=e1[,RNAseq_gl=='trt']
e1_b=e1[,RNAseq_gl=='untrt']
a_filter=apply(e1_a, 1,function(x) sd(x)>0)
b_filter=apply(e1_b, 1,function(x) sd(x)>0)
table(a_filter,b_filter)
e1=e1[a_filter | b_filter,]

Q3: 对过滤后的e1矩阵进行每一行独立根据分组矩阵进行T检验

dat <- apply(e1, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
})
class(dat)
df <- data.frame(expression = dat)

Q4:对e1矩阵进行加1后log2的归一化命名为e2再对每一行独立根据分组矩阵进行T检验,对e1,e2的T检验P值做相关性分析

p1=apply(e1, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
}) 
e2=log(e1+1)
p2=apply(e2, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
}) 
plot(p1,p2,cex = 0.5, pch = 21, col = "red", bg = "blue")
abline(lsfit(p1,p2),col = "yellow", lty = 3)
cor(p1,p2)
原始数据与log2转化后数据的每组基因t.test检验p值的相关性点图

统计检验中的类

  • 对于最后一道习题中出现t.test(x~RNAseq_gl)$p.value这样的代码，我想知道t.test(x~RNAseq_gl)返回的是什么对象，因为一般来说数据框和列表常用$符号取值。用习题中的数据进行了验证。
class(t.test(e1[1,]~RNAseq_gl))
[1] "htest"
> t.test(e1[1,]~RNAseq_gl)

    Welch Two Sample t-test

data:  e1[1, ] by RNAseq_gl
t = -1.3886, df = 5.2655, p-value = 0.2208
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -695.2741  202.7741
sample estimates:
  mean in group trt mean in group untrt 
             618.75              865.00
  • 可见t.test返回的结果为htest对象，进行搜索后有一大致的了解：
    在作统计检验的时候，输出的结果都属于一个叫htest的类。t.test的结果其实就是个列表，但是列表不是好的表现方法，于是通过htest这个类，换了一种输出到屏幕的方法。
    察看一下t.test的结果的类：class(y)
    得到 ：[1] "htest"。
    去掉y的类：unclass(y)
    得到：
 <br />
$statistic<br />
       t <br />
7.002526 <br />
<br />
$parameter<br />
df <br />
12 <br />
<br />
$p.value<br />
[1] 1.429288e-05<br />
<br />
$conf.int<br />
[1] 0.4604719 0.8764511<br />
attr(,"conf.level")<br />
[1] 0.95<br />
<br />
$estimate<br />
mean of x <br />
0.6684615 <br />
<br />
$null.value<br />
mean <br />
   0 <br />
<br />
$alternative<br />
[1] "two.sided"<br />
<br />
$method<br />
[1] "One Sample t-test"<br />
<br />
$data.name<br />
[1] "x - 1.44"<br />
  • 可见，要取出t.test返回对象中的值，可以用$的办法。

本练习题来自于生信技能树培训。
最后的htest类参考：https://d.cosx.org/d/8800-8800

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 214,172评论 6赞 493
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 91,346评论 3赞 389
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 159,788评论 0赞 349
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 57,299评论 1赞 288
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 66,409评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,467评论 1赞 292
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,476评论 3赞 412
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,262评论 0赞 269
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,699评论 1赞 307
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,994评论 2赞 328
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,167评论 1赞 343
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,827评论 4赞 337
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,499评论 3赞 322
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,149评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,387评论 1赞 267
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 47,028评论 2赞 365
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,055评论 2赞 352