写在前面——本文是介绍了一些统计学的基础知识，以及如何用R语言实现。旨在熟练使用编程语言解决实际问题。统计学是一门复杂的学科，还需更加系统的学习。下方 生统基础是大佬写的StartQuest的读书笔记，能够帮助快速理解生物统计学。

题目原文：http://www.bio-info-trainee.com/4385.html
参考答案：https://www.jianshu.com/p/1c0cbf1bfb75
参考答案：https://hc1023.github.io/2019/09/27/R-statistic/
生统基础：http://rvdsd.top/categories/%E7%94%9F%E7%BB%9F%E4%B9%8BStatQuest/

基础概念

1. 载入R中自带的数据集 iris，指出其每列是定性还是定量数据

colnames(iris)
head(iris)

基础概念讲解：https://mp.weixin.qq.com/s/OtB2h6f00U2SRZLzveJKfQ
前四个为定量变量，最后一个为定性变量
​

2. 对数据集 iris的所有定量数据列计算集中趋势指标：众数、分位数和平均数

# 求前四个变量的最大值、最小值、分位数和平均数
summary(iris[,1:4])# R中无直接求众数的函数
# 法一
# tmp <- table(iris[,1])
# index <- which.max(tmp)
# tmp[index]# 法二
getmode <- function(x){return(as.numeric(names(table(x))[table(x) == max(table(x))]))
}# table(iris[,1]) == max(table(iris[,1]))
# names(table(iris[,1]))
# names(table(iris[,1]))[table(iris[,1]) == max(table(iris[,1]))]for(i in 1:4){print(getmode(iris[,i]))
}# 第三列有两个众数

​

3. 对数据集 iris的所有定性数据列计算水平及频次

summary(iris[,5])#output#
setosa versicolor  virginica 50         50         50

​

4. 对数据集 iris的所有定量数据列计算离散趋势指标：方差和标准差等

# 方差
apply(iris[,1:4], 2, var)
# 标准差
apply(iris[,1:4], 2, sd)# 通过psych包中的describe()计算描述性统计量
# 输出结果依次为：非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误
# install.packages("psych")
library(psych)
describe(iris[,1:4])

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

cor(iris[,1:2], method = "pearson")
cor(iris[,1:2], method = "spearman")
cor(iris[,1:2], method = "kendall")

Pearson积差相关系数衡量了两个定量变量之间的线性相关程度（默认值）
Spearman等级相关系数则衡量分级定序变量之间的相关程度
Kendall’sTau相关系数也是一种非参数的等级相关度量

​

6. 对数据集 iris的所有定量数据列内部z-score标准化，并计算标准化后每列的平均值和标准差

newdata <- apply(iris[1:4],2,scale)
apply(newdata,2,mean)
apply(newdata,2,sd)

图源：https://blog.csdn.net/Orange_Spotty_Cat/article/details/80312154
​

scale(x, center = TRUE, scale = TRUE)
center和scale默认为真，center为真表示数据中心化，scale为真表示数据标准化
数据中心化是指：变量减去它的均值
数据标准化是指：数值减去均值，再除以标准差

所以不难理解，对数据的每列调用scale就完成了z-score标准化
标准化之后的数据均值为0，标准差为1（上面的均值结果近似于0，也算为0吗？）
​

参考：https://zhuanlan.zhihu.com/p/32482328
​

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

cor(newdata[,1:2], method = "pearson")
cor(newdata[,1:2], method = "spearman")
cor(newdata[,1:2], method = "kendall")

与第5题结果相同，标准化后的数据相关性不变
​

8. 根据数据集 iris的第五列拆分数据集后重复上面的Q2到Q7问题

# 划分数据集
table(iris$Species)
setosa <- iris[iris$Species=="setosa",]
head(setosa)
versicolor <- iris[iris$Species=="versicolor",]
virginica <- iris[iris$Species=="virginica",]# 写个无脑输出函数
getResult <- function(x){print("### 2 ###")print(apply(x[,1:4], 2, getmode))print(summary(x[,1:4]))print("### 3 ###")print(summary(x[,5]))print("### 4 ###")print(apply(x[,1:4], 2, var))print(apply(x[,1:4], 2, sd))describe(x[,1:4])print("### 5 ###")print(cor(x[,1:2], method = "pearson"))print(cor(x[,1:2], method = "spearman"))print(cor(x[,1:2], method = "kendall"))print("### 6 ###")newx <- apply(x[1:4],2,scale)print(apply(newx,2,mean))print(apply(newx,2,sd))print("### 7 ###")print(cor(newx[,1:2], method = "pearson"))print(cor(newx[,1:2], method = "spearman"))print(cor(newx[,1:2], method = "kendall"))
}getResult(setosa)
getResult(versicolor)
getResult(virginica)

​

9. 载入R中自带的数据集 mtcars，重复上面的Q1到Q7个问题

# q1
> colnames(mtcars)[1] "mpg"  "cyl"  "disp" "hp"   "drat" "wt"   "qsec" "vs"   "am"   "gear" "carb"# mtcars每列数据的含义
Format
A data frame with 32 observations on 11 (numeric) variables.[, 1]   mpg Miles/(US) gallon
[, 2]   cyl Number of cylinders
[, 3]   disp    Displacement (cu.in.)
[, 4]   hp  Gross horsepower
[, 5]   drat    Rear axle ratio
[, 6]   wt  Weight (1000 lbs)
[, 7]   qsec    1/4 mile time
[, 8]   vs  Engine (0 = V-shaped, 1 = straight)
[, 9]   am  Transmission (0 = automatic, 1 = manual)
[,10]   gear    Number of forward gears
[,11]   carb    Number of carburetors# 定量变量一般使用具体的数字表示，定性变量也称分类变量，可以用来分类
# 在这里，8和9明显为定性变量
# 其余列为定量变量

# q2-q7
# 将定量变量取出，放在一起
newmtcars <- cbind(mtcars[,1:7],mtcars[10:11])# 没设置变量，直接暴力输出
getMtcars <- function(){print("### 2 ###")print(apply(newmtcars, 2, getmode))print(summary(newmtcars))print("### 3 ###")print(apply(mtcars[8:9], 2, table))print("### 4 ###")print(apply(newmtcars, 2, var))print(apply(newmtcars, 2, sd))describe(newmtcars)print("### 5 ###")print(cor(mtcars[,1:2], method = "pearson"))print(cor(mtcars[,1:2], method = "spearman"))print(cor(mtcars[,1:2], method = "kendall"))print("### 6 ###")newx <- apply(newmtcars,2,scale)print(apply(newx,2,mean))print(apply(newx,2,sd))print("### 7 ###")print(cor(newx[,1:2], method = "pearson"))print(cor(newx[,1:2], method = "spearman"))print(cor(newx[,1:2], method = "kendall"))
}getMtcars()# mpg和cyl的负相关性还是很强的

​

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

summarizedexperiment详细说明：https://bioconductor.org/packages/release/bioc/vignettes/SummarizedExperiment/inst/doc/SummarizedExperiment.html#anatomy-of-a-summarizedexperiment

其中关于airway的描述：
该软件包包含来自气道平滑肌每个基因的读取计数的 RNA-Seq实验的示例数据集。这些数据存储在一个对象中，该对象包含 8 个不同的实验和分析 64,102 个基因转录本。

人性化解读：https://www.dazhuanlan.com/nby/topics/1127479

rm(list=ls())
options(stringsAsFactors = F)# BiocManager::install("airway")
library(airway)
data("airway")
airway # 获取实验数据
RNAseq_expr <- assay(airway)
colnames(RNAseq_expr)
RNAseq_expr[1:4,1:4] # RNAseq_gl <- colData(airway)[,3]
# table(RNAseq_gl)# 计算相关性
M <- cor(RNAseq_expr)
pheatmap::pheatmap(M)

表达矩阵相关

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

# 清空环境，重新取一下数据
rm(list=ls())
options(stringsAsFactors = F)library(airway)
data("airway")RNAseq_expr <- assay(airway)
head(RNAseq_expr)# +1是因为0没法取对数
tmp <- log2(RNAseq_expr[,1]+1)
head(tmp)mean(tmp)
sd(tmp)

​

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

# dnorm给出密度，pnorm给出分布函数，qnorm给出分位数函数，rnorm产生随机偏差
a <- rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
# 升序排序
a <- sort(a)
plot(a)# 跟处理前的数据分布对比一下
points(sort(tmp), col = "red")

​

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

# 删除对应数据
NewRNAseq_expr <- RNAseq_expr[RNAseq_expr[,1]>=5, ]# Q1
# 没有+1是因为数据不存在0了
tmp <- log2(NewRNAseq_expr[,1])
head(tmp)
mean(tmp)
sd(tmp)# Q2
b <- rnorm(length(tmp), mean = mean(tmp), sd = sd(tmp))
b <- sort(b)
plot(b)
points(sort(tmp), col = "red")# 可以看到，筛选过的数据就很接近正态分布了

​

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

# 选择数据进行T检验
x <- RNAseq_expr[RNAseq_expr[,1]>=5, ]
y <- RNAseq_expr[RNAseq_expr[,2]>=5, ]
x <- log2(x[,1])
y <- log2(y[,2])t.test(x,y)# 绘制箱图，查看数据分布情况
library(ggplot2)
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")

​

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

# 找到和最大的那一行
sumMax <- which.max(rowSums(RNAseq_expr))
sumMax
RNAseq_expr[sumMax, ]# 分组情况
RNAseq_gl <- colData(airway)[,3]
table(RNAseq_gl)# “~”左边的是一个数值型变量，右边是一个二分变量
t.test(RNAseq_expr[sumMax, ]~RNAseq_gl)

​

6. 取RNAseq_expr的MAD（绝对中位差）最大的那一行根据分组矩阵进行T检验

绝对中位差实际求法是用原数据减去中位数后得到的新数据的绝对值的中位数。但绝对中位差常用来估计标准差，估计标准差=1.4826*绝对中位差。R语言中返回的是估计的标准差。

madMax <- which.max(apply(RNAseq_expr, 1, mad))
madMax
RNAseq_expr[madMax, ]t.test(RNAseq_expr[madMax, ]~RNAseq_gl)

​

7. 对RNAseq_expr全部加1后取log2后重复Q5和Q6

RNAseq_expr_log2 <- log2(RNAseq_expr+1)# Q5
# 基因没变，P值变化不大
sumMax <- which.max(rowSums(RNAseq_expr_log2))
sumMax
RNAseq_expr_log2[sumMax, ]t.test(RNAseq_expr_log2[sumMax, ]~RNAseq_gl)# Q6
# 基因变了，P值变化很大
# mad可以反映数据的波动情况，只是粗略的一个估计值，辅助判断
madMax <- which.max(apply(RNAseq_expr_log2, 1, mad))
madMax
RNAseq_expr_log2[madMax, ]t.test(RNAseq_expr_log2[madMax, ]~RNAseq_gl)

​

8. 取RNAseq_expr矩阵的MAD最高的100行，对列和行分别进行层次聚类

# 取log2之后的数据
cg1 <- names(tail(sort(apply(RNAseq_expr_log2, 1, mad)), 100))
dat1 <- RNAseq_expr_log2[cg1, ]# 原始数据
cg2 <- names(tail(sort(apply(RNAseq_expr, 1, mad)), 100))
dat2 <- RNAseq_expr[cg2, ]# 对行进行聚类
plot(hclust(dist(t(dat1))))
plot(hclust(dist(t(dat2))))# 对列进行聚类
plot(hclust(dist(dat1)))
plot(hclust(dist(dat2)))# 聚类之后的图不知道怎么解读...挖坑待填...

9. 取RNAseq_expr矩阵的SD最高的100行，对列和行分别进行层次聚类

cg <- names(tail(sort(apply(RNAseq_expr_log2, 1, sd)), 100))
dat <- RNAseq_expr_log2[cg, ]plot(hclust(dist(t(dat))))plot(hclust(dist(dat)))

​

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

# 紧接Q8
pheatmap::pheatmap(scale(dat1))
pheatmap::pheatmap(t(scale(t(dat1))))# 图很乱，想看清楚还需再处理一下

统计检验相关

这里需要对前面的RNAseq_expr矩阵进行一定程度的过滤

rm(list=ls())
options(stringsAsFactors = F)library(airway)
data("airway")
RNAseq_expr <- assay(airway)
colnames(RNAseq_expr)
dim(RNAseq_expr)
RNAseq_gl <- colData(airway)[,3]
table(RNAseq_gl)# 删除 每列相加，和小于等于1 的数据
e1 <- RNAseq_expr[apply(RNAseq_expr,1,function(x) sum(x>0)>1),]
dim(e1)#output#
> dim(RNAseq_expr)
[1] 64102     8
> dim(e1)
[1] 28877     8

上面的 e1 矩阵下面的习题就用得到。

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

apply(e1, 1, function(x){t.test(x ~ RNAseq_gl)$p.value
})#output#Error in t.test.default(x = c(SRR1039509 = 1L, SRR1039513 = 1L, SRR1039517 = 1L,  : data are essentially constant

# 使用二分法手动查找出错的位置
> t.test(e1[4257,] ~ RNAseq_gl)$p.value
Error in t.test.default(x = c(SRR1039509 = 1L, SRR1039513 = 1L, SRR1039517 = 1L,  : data are essentially constant
> rownames(e1)[4257]
[1] "ENSG00000116852"
> e1[4257,]
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 SRR1039521 0          1          0          1          0          1          0          1 # 通过观察数据可以发现，报错应该是因为相同组的数据都一样，没法进行T检验

​

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

# 分类
e1_a <- e1[,RNAseq_gl=='trt']
e1_b <- e1[,RNAseq_gl=='untrt']# 筛选，如果数据相同的话，行标准差为0
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)
table(b_filter)
table(a_filter,b_filter)# 去除同时为FALSE的值
e1 <- e1[a_filter | b_filter, ]

​

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

p1 <- apply(e1, 1, function(x){t.test(x~RNAseq_gl)$p.value
})
head(p1)#output#
ENSG00000000003 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000938 0.2208266       0.8275543       0.6747468       0.3841246       0.2151699
ENSG00000000971 0.4949860

​

4. 对e1矩阵进行加1后log2的归一化命名为e2再对每一行独立根据分组矩阵进行T检验

e2 <- log2(e1+1)
p2 <- apply(e2, 1, function(x){t.test(x ~ RNAseq_gl)$p.value
})
head(p2)#output#
ENSG00000000003 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000938 0.1705491       0.9265533       0.5472763       0.4442560       0.1975805
ENSG00000000971 0.6024837

​

5. 对e1,e2的T检验P值做相关性分析

plot(p1, p2, cex = 0.5, pch = 21, col = "red", bg = "blue")# 画出最小二乘拟合直线
abline(lsfit(p1, p2), col = "yellow", lty = 3)
cor(p1,p2) #output#
[1] 0.9470391

生信学习——基于R的统计习题（附详细答案解读）相关推荐

1. 生信学习——基于R的可视化习题30个（附详细答案解读）

   题目目录 一.基础绘图 1. 对RNAseq_expr的每一列绘制boxplot图 2. 对RNAseq_expr的每一列绘制density图 3. 对RNAseq_expr的每一列绘制条形图 4. ...

2. 生信学习——R语言练习题-初级（附详细答案解读）

   题目目录 1. 打开 Rstudio 告诉我它的工作目录. 2. 新建6个向量,基于不同的数据类型.(重点是字符串,数值,逻辑值) 3. 告诉我在你打开的rstudio里面 getwd() 代码运行后 ...

3. 生信学习——生信人的20个R语言习题（上）（附详细答案解读）

   题目目录 1. 安装一些R包. 2. 了解ExpressionSet对象,比如CLL包里面就有data(sCLLex),找到它包含的元素,提取其表达矩阵(使用exprs函数),查看其大小. 3. 了解 ...

4. 生信学习——生信人的20个R语言习题（下）（附详细答案解读）

   题目目录 12. 理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标,最后按照mad值排序,取top 50 mad值的基因,得到列表. ...

5. 生信学习——R语言小作业-中级（附详细答案解读）

   题目目录 1. 请根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol). 2. 根据R包hgu133a.db找到下面探针对应的基因名(symbol). 3. ...

6. 生信学习——R语言学习总结

   写在前面--经过了四十天断断续续的学习,算是对R语言有了初步的了解.其实使用R语言,无非就是对数据进行处理分析,然后把结果可视化.但是数据的千变万化,还有数以万计的函数.数据格式,使得这个过程变得很复 ...

7. 送书 | 知乎阅读300w+的生信学习指南（更新版）

   先送书 在上周的留言送书活动中,恭喜下面这位读者获得书籍"Oracle高性能系统架构实战大全",请及时与生信宝典编辑(shengxinbaodian)联系. 2020过去三分之一了 ...

8. 生信学习学的是什么？常识！

   生物信息学学的是什么?常识! 学习的是基本生物学概念的常识! 学习的是计算机基础的常识! 学习的是图形解读的常识! 学习的是统计的常识! 拦住生信学习脚步的不是技术有多难,而是有些常识你还不知道. 这 ...

9. 知乎阅读三百万的生信学习指南

   作为本科学生物,硕博转行生物信息的人,经常会被人问起,为啥学习生物信息了呢?这背后通常会带着一些困惑,生物信息分析好不好学? 生信的作用越来越大,想学的人越来越多,不管是为了以后发展,还是为了解决眼下 ...

10. 《机器学习与数据科学（基于R的统计学习方法）》——1.2　机器学习的实际案例...

    本节书摘来异步社区<机器学习与数据科学(基于R的统计学习方法)>一书中的第1章,第1.2节,作者:[美]Daniel D. Gutierrez(古铁雷斯),更多章节内容可以访问云栖社区&q ...

最新文章

1. Full Gc经历分析
2. 批量插入以及数据存在重复就进行更新操作
3. Java程序员需要掌握的计算机底层知识（三）：进程、线程、纤程、中断
4. Angular Component代码和编译后生成的JavaScript代码
5. android webviewclient 点击事件,Android Api WebViewClient 详细解析
6. Android 8.0 学习（24）---Android8.0 WiFi热点适配
7. docker安装软件（vim,service）
8. 阿里云弹性计算负责人蒋林泉：亿级场景驱动的技术自研之路 | 问底中国 IT 技术演变...
9. Linux连接锐捷校园网客户端
10. 服务器运维软硬件维护月报,运维月报ppt
11. 【Excel】在单元格中插入换行符
12. STEAM上的一款电路模拟神器 — CRUMB Circuit Simulator
13. 文件(图片)上传保存与展示
14. matlab求平均聚集系数,复杂网络聚类系数和平均路径长度计算的MATLAB源代码
15. CLOCs：一种相机-激光雷达3D目标检测后融合方法
16. UCI计算机工程必修专业课,想问问加州大学欧文分校计算机工程专业怎么样？
17. 转：Metalink 账户
18. 两个整数相除，不使用乘法，除法和取余
19. 一面准备：项目经历、开放问题回答、基础
20. 《谋圣鬼谷子》曝概念海报

热门文章

1. 第七章 浮动（制作京东登录页面）
2. 深度学习图像处理目标检测图像分割计算机视觉 10--通用场景图像分割
3. JAVA基础之单例模式
4. 为什么入侵地球的总是火星人，是别的星球不够努力吗？
5. 贪婪的大脑:为何人类会无止境地寻求意义 目录
6. 飞桨AI Studio之加州房价预测——机器学习的Hello world
7. Oracle AutoVue 安装与配置说明
8. 我总结了程序员转行得最有出路5个方向
9. 苹果MAC电脑双系统教程——MAC安装Windows双系统教程
10. String spilt的使用和常见的坑以及平替方案性能比较