R语言作业-统计30题

题目

链接：http://www.bio-info-trainee.com/4385.html

我做题的时候主要翻阅学习了《R语言实战》里统计相关内容。

基础概念

需要掌握R内置数据集及R包数据集

定性/定量数据

鸢尾花(iris)数据集，包含150个鸢尾花的信息，共五列，分别为萼片长度(Sepal.Length)、萼片宽度(Sepal.Width)、花瓣长度(Petal.Length)、花瓣宽度(Petal.Width)和种类(Species)。前四列为定量数据，后一列种类为定性数据，是非连续的字符变量。

data("iris")
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

colnames(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
## [5] "Species"

table(iris$Species)
##     setosa versicolor  virginica 
##         50         50         50


描述性数据

定量数据的集中趋势指标主要是：众数、分位数和平均数
定量数据的离散趋势指标主要是：极差，方差和标准差，标准分数，相对离散系数（变异系数），偏态系数与峰态系数

# 平均数
apply(iris[,1:4],2,mean)
# 中位数
apply(iris[,1:4],2,median)
# 标准差
apply(iris[,1:4],2,sd)
# 方差
apply(iris[,1:4],2,var)
# 绝对中位数(median absolute deviation)
apply(iris[,1:4],2,mad)
# 众数
mode <- function(x) {
  uniq <- unique(x)
  tab <- tabulate(match(x,uniq))
  uniq[tab == max(tab)]
}
apply(iris[,1:4],2,mode)
## $Sepal.Length
##      mean    median        sd       var       mad      mode 
## 5.8433333 5.8000000 0.8280661 0.6856935 1.0378200 5.0000000 
## 
## $Sepal.Width
##      mean    median        sd       var       mad      mode 
## 3.0573333 3.0000000 0.4358663 0.1899794 0.4447800 3.0000000 
## 
## $Petal.Length
##     mean   median       sd      var      mad    mode1    mode2 
## 3.758000 4.350000 1.765298 3.116278 1.853250 1.400000 1.500000 
## 
## $Petal.Width
##      mean    median        sd       var       mad      mode 
## 1.1993333 1.3000000 0.7622377 0.5810063 1.0378200 0.2000000

# 可以写个函数一起算
iris_count <- function(x){
  mean <- mean(x)
  median <- median(x)
  sd <- sd(x)
  var <- var(x)
  mad <- mad(x)
  uniq <- unique(x)
  tab <- tabulate(match(x,uniq))
  mode <- uniq[tab == max(tab)]
  return(c(mean=mean,median=median,sd=sd,var=var,mad=mad,mode=mode))
}

apply(iris[,1:4],2,iris_count)
# pastecs包中的stat.desc()函数也可以计算等等

# 定性数据频次
table(iris[,5])  
##     setosa versicolor  virginica 
##         50         50         50

# summary()函数提供了最大值、最小值、四分位数和数值变量的均值
summary(iris)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50



分组统计，最一开始是想将数据集分成三个数据框，重复之前函数，有几个方法：

# 方法一：which()定位函数 cbind()按列合并对象
setosa <- iris[which(iris$Species == 'setosa'),]

# 方法二：subset()
setosa2 <- subset(iris,iris$Species == 'setosa')

# 方法三：split()--> sapply()--summary()
iris_s <- split(iris,iris$Species)
setosa3 <- iris_s$setosa



或者不分开，之前对原数据集计算：

# aggregate()函数分开一个个计算统计数
iris_mean <- aggregate(iris[,1:4],by = list(iris$Species),FUN = mean)
iris_med <- aggregate(iris[,1:4],by = list(iris$Species),FUN = median)
iris_sd <- aggregate(iris[,1:4],by = list(iris$Species),FUN = sd)
iris_var <- aggregate(iris[,1:4],by = list(iris$Species),FUN = var)

# 整体计算，利用之前的函数
# 分组计算
iris_count2 <- function(x)sapply(x,iris_count)
by(iris[,1:4], iris$Species,iris_count2)


apply函数家族

apply函数可以解决数据循环处理的问题，可以对矩阵、数据框、数组（二维、多维），按行或列进行循环计算，对子元素进行迭代，并把子元素以参数形式给自定义的FUN函数中，并返回计算结果。

函数定义:
apply(X,MARGIN,FUN,...)
参数列表：

  • X:数组、矩阵、数据框
  • MARGIN：1表示按行，2表示按列
  • FUN：自定义调用函数
  • ...
lapply函数

用来对list、data.frame进行循环，并返回和X长度同样的list结构作为结果集。

sapply函数

同lapply函数，多了2个参数simplify和USE.NAMES，返回值为向量，不是list对象。

vapply函数

类似sapply函数，提供了FUN.VALUE参数，用来控制返回值的行名。

mapply函数

类似sapply函数，第一个参数为FUN，可接受多个数据。

tapply函数

tapply函数用于分组的循环计算，相当于group by的操作。

函数定义：
tapply(X,INDEX,FUN,simplify,...)

参数列表：

  • X：向量
  • INDEX：用于分组的索引
  • FUN：自定义的调用函数
  • simplify：是否数组化
rapply函数

只处理list类型数据，对list的每个元素进行递归遍历，如果list包括子元素则继续遍历。


相关性

R可以计算多种相关系数，包括Pearson相关系数、Spearman相关系数、Kendall相关系数、偏相关系数、多分格相关系数、多系列相关系数。cor()函数可以计算前三种相关系数，cov()函数可以计算协方差。

# Pearson积差相关系数衡量了两个定量变量之间的线性相关程度。
# Spearman等级相关系数则衡量分级定序变量之间的相关程度。
# Kendall’s Tau相关系数也是一种非参数的等级相关度量。

# 计算数据集 iris的前两列变量的相关性
cor(iris[,1:2],use = "everything", method = "pearson")
##              Sepal.Length Sepal.Width
## Sepal.Length    1.0000000  -0.1175698
## Sepal.Width    -0.1175698   1.0000000

cor(iris[,1:2],use = "everything", method = "kendall")
##              Sepal.Length Sepal.Width
## Sepal.Length   1.00000000 -0.07699679
## Sepal.Width   -0.07699679  1.00000000

cor(iris[,1:2],use = "everything", method = "spearman")
##              Sepal.Length Sepal.Width
## Sepal.Length    1.0000000  -0.1667777
## Sepal.Width    -0.1667777   1.0000000

cov(iris[,1:2])
##              Sepal.Length Sepal.Width
## Sepal.Length    0.6856935  -0.0424340
## Sepal.Width    -0.0424340   0.1899794

# 不指定默认是Person系数
by(iris[,1:2], iris$Species,cor)
## iris$Species: setosa
##              Sepal.Length Sepal.Width
## Sepal.Length    1.0000000   0.7425467
## Sepal.Width     0.7425467   1.0000000
## -------------------------------------------------------- 
## iris$Species: versicolor
##              Sepal.Length Sepal.Width
## Sepal.Length    1.0000000   0.5259107
## Sepal.Width     0.5259107   1.0000000
## -------------------------------------------------------- 
## iris$Species: virginica
##              Sepal.Length Sepal.Width
## Sepal.Length    1.0000000   0.4572278
## Sepal.Width     0.4572278   1.0000000


标准化

数据的标准化是指中心化之后的数据在除以数据集的标准差，即数据集中的各项数据减去数据集的均值再除以数据集的标准差。scale()函数可以完成标准化。

iris_z <- apply(iris[,1:4], 2, scale)
apply(iris_z[,1:4],2,iris_count)
## $Sepal.Length
##          mean        median            sd           var           mad 
## -4.484318e-16 -5.233076e-02  1.000000e+00  1.000000e+00  1.253306e+00 
##          mode 
## -1.018437e+00 
## 
## $Sepal.Width
##          mean        median            sd           var           mad 
##  2.034094e-16 -1.315388e-01  1.000000e+00  1.000000e+00  1.020451e+00 
##          mode 
## -1.315388e-01 
## 
## $Petal.Length
##          mean        median            sd           var           mad 
## -2.895326e-17  3.353541e-01  1.000000e+00  1.000000e+00  1.049823e+00 
##         mode1         mode2 
## -1.335752e+00 -1.279104e+00 
## 
## $Petal.Width
##          mean        median            sd           var           mad 
## -3.663049e-17  1.320673e-01  1.000000e+00  1.000000e+00  1.361544e+00 
##          mode 
## -1.311052e+00

# 分组
by(iris_z[,1:4], iris$Species,iris_count2)

# 标准化后的相关性和原数据结果相同
cor(iris_z[,1:2],use = "everything", method = "pearson")
##              Sepal.Length Sepal.Width
## Sepal.Length    1.0000000  -0.1175698
## Sepal.Width    -0.1175698   1.0000000

mtcars

mtcars数据集是32辆汽车在11个指标上的数据。

data("mtcars")
head(mtcars) #32辆汽车在11个指标上的数据
summary(mtcars)

# tapply 根据gear列分组，统计mpg值
tapply(mtcars$mpg,mtcars$gear, summary)
# 求标准差
apply(mtcars,2,sd)
# 前两列的相关性
cor(mtcars[,1:2],use = "everything", method = "pearson")
# 标准化
mtcars_z <- apply(mtcars, 2, scale)
# 求标准化后的每列的平均值
apply(mtcars_z,2,mean)


表达矩阵相关

airway包是8个样本的RNA-seq数据的counts矩阵，这8个样本分成2组，每组是4个样本，分别是 trt 和 untrt 组。

options(stringsAsFactors = F)
suppressMessages(library(airway))
data(airway)
# 这里需要自行学习bioconductor里面的RangedSummarizedExperiment对象
RNAseq_expr=assay(airway)
colnames(RNAseq_expr) 
## [1] "SRR1039508" "SRR1039509" "SRR1039512" "SRR1039513" "SRR1039516"
## [6] "SRR1039517" "SRR1039520" "SRR1039521"

RNAseq_expr[1:4,1:4] 
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513
## ENSG00000000003        679        448        873        408
## ENSG00000000005          0          0          0          0
## ENSG00000000419        467        515        621        365
## ENSG00000000457        260        211        263        164
# RNAseq_expr 是一个数值型矩阵，属于连续性变量，可以探索众数、分位数和平均数 ，极差，方差和标准差等统计学指标

RNAseq_gl=colData(airway)[,3]
table(RNAseq_gl)
## RNAseq_gl
##   trt untrt 
##     4     4

# 把RNAseq_expr第一列全部加1后取log2后计算平均值和标准差
tmp=log2(RNAseq_expr[,1]+1)
mean(tmp)
## [1] 2.251721
sd(tmp)
## [1] 3.64586

# 根据上一步得到平均值和标准差生成同样个数的随机的正态分布数值
a=rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
a=sort(a) # 排序

# 删除RNAseq_expr第一列低于5的数据后，重复Q1和Q2
tmp = RNAseq_expr[,1][which(RNAseq_expr[,1]>5)]
tmp=log2(tmp)
a=rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
a=sort(a)
plot(a)
points(sort(tmp))

# 基于Q3对RNAseq_expr的第一列和第二列进行T检验
x = tmp
y = log2(RNAseq_expr[,2][which(RNAseq_expr[,2]>5)])
t.test(x,y)
##  Welch Two Sample t-test
## 
## data:  x and y
## t = 2.8841, df = 34335, p-value = 0.003928
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.02800383 0.14680584
## sample estimates:
## mean of x mean of y 
##  7.651705  7.564300
# t  t统计值；df 自由度；p H0假设成立的概率
# 基因表达分析的零假设是： 基因在不同处理下的表达量相同

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")
x ~ y.png
# 取RNAseq_expr行之和最大的那一行根据分组矩阵进行T检验
pos=which.max(rowSums(RNAseq_expr))
pos
t.test(RNAseq_expr[pos,]~RNAseq_gl)


# 取RNAseq_expr的MAD最大的那一行根据分组矩阵进行T检验
pos=which.max(apply(RNAseq_expr,1,mad))
t.test(RNAseq_expr[pos,]~RNAseq_gl)

# 对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)


# 取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 )))
聚类-列-mad.png
# 取RNAseq_expr矩阵的SD最高的100行，对列和行分别进行层次聚类
cg=names(tail(sort(apply(RNAseq_expr,1,sd)),100))
dat=RNAseq_expr[cg,]
plot(hclust(dist(t(dat))))
colnames(dat)
RNAseq_gl
plot(hclust(dist( dat )))
聚类-列-sd.png

t检验

t检验是一种可用于比较的假设检验。

理解t检验：一个年纪共有好多学生，需要研究他们的平均身高。这时，这批学生是我们要研究的对象，即总体。从这个年纪中每个班级随机挑选10名同学，这部分同学则为样本，通过样本来对总体的某个统计特征（比如上面研究的平均值、众数、方差等）做判断的方法为假设检验。

独立样本t检验

一个针对两组的独立样本t检验可以用于检验两个总体的均值相等的假设，检验调用格式为：
t.test( y ~ x, data )
其中y是一个数值型变量，x是一个二分变量。

t.test(y1,y2)
其中y1、y2为数值型向量。

现在还不能用自己的语言解释清楚，整合几篇写的比较详细的教程：
http://www.biye5u.com/article/R/2019/6399.html
//www.greatytc.com/p/67be9b3806cd

options(stringsAsFactors = F)
suppressMessages(library(airway))
data(airway)
# 这里需要自行学习bioconductor里面的RangedSummarizedExperiment对象
RNAseq_expr=assay(airway)
RNAseq_gl=colData(airway)[,3]
e1=RNAseq_expr[apply(RNAseq_expr,1,function(x) sum(x>0)>1),] 

# 对e1每一行独立根据分组矩阵进行T检验，检查为什么有些行T检验失败
apply(e1, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
})

# 报错：data are essentially constant

# 找出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)
##         b_filter
## a_filter FALSE  TRUE
##    FALSE    10   948
##    TRUE    603 27316

e1=e1[a_filter | b_filter,]

# 对过滤后的e1矩阵进行每一行独立根据分组矩阵进行T检验
p1=apply(e1, 1, function(x){
  t.test(x~RNAseq_gl)$p.value
})

# 对e1矩阵进行加1后log2的归一化命名为e2再对每一行独立根据分组矩阵进行T检验
e2=log(e1+1)
p2=apply(e2, 1, function(x){
 t.test(x~RNAseq_gl)$p.value
}) 

# 对e1,e2的T检验P值做相关性分析
plot(p1,p2)
cor(p1,p2)
## [1] 0.9470391
相关性.png

统计这部分还是似懂非懂的状态，继续学习~

更多学习资源：
生信技能树公益视频合辑
生信技能树简书账号
生信工程师入门最佳指南
生信技能树全球公益巡讲
招学徒
...
你的宣传能让数以万计的初学者找到他们的家，技能树平台一定不会辜负每一个热爱学习和分享的同道中人 😎

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