吨吨吨的生信技能树r语言习题（2）

统计专题30题

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

前四列都是定量最后一列是定性

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

注意：apply函数要求对象为matrix

iris2=as.matrix(iris[,c(1:4)])
apply(iris2,2,mean)
apply(iris2,2,quantile)
apply(iris2,2,function(x){as.numeric(names(table(x))[table(x) == max(table(x))])})

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

table(iris[,5])

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

apply(iris2,2,sd)
apply(iris2,2,var)

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

有Pearson，spearman和kendall三种可选。其中pearson对数据要求较高。
Pearson系数、Spearman系数、Kendall系数使用条件 | 挖数网 (ikddm.com)

cor(iris$Sepal.Length,iris$Sepal.Width,method = 'pearson')

Q6:对数据集 iris的所有定量数据列内部zcore标准化，并计算标准化后每列的平均值和标准差

zscore=as.matrix(scale(iris[,c(1:4)]))
apply(zscore,2,mean)
apply(zscore,2,sd)

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

cor(zscore[,1],zscore[,2],method = 'pearson')

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

library(dplyr)
datamean=group_by(iris[,c(1:4)], iris$Species) %>% summarize_each(funs(mean))
library(dplyr)
datamean=group_by(iris[,c(1:4)], iris$Species) %>% summarize_each(funs(mean))
dataquantile=group_by(iris[,c(1:4)], iris$Species) %>% summarize_each(funs(quantile))
datazhongshu=group_by(iris[,c(1:4)], iris$Species) %>% summarize_each(funs(zhongshu))
datasd=group_by(iris[,c(1:4)], iris$Species) %>% summarize_each(funs(sd))
datavar=group_by(iris[,c(1:4)], iris$Species) %>% summarize_each(funs(var))
cor(dat1$Sepal.Length,dat1$Sepal.Width,method = 'pearson')
cor(dat2$Sepal.Length,dat2$Sepal.Width,method = 'pearson')
cor(dat3$Sepal.Length,dat3$Sepal.Width,method = 'pearson')
datascale=group_by(iris[,c(1:4)], iris$Species) %>% summarize_each(funs(scale))
datasd=group_by(datascale, iris$Species) %>% summarize_each(funs(sd))
datamean=group_by(datascale, iris$Species) %>% summarize_each(funs(mean))
cor(dat1[,1],dat1[,2],method = 'pearson')
cor(dat2[,1],dat2[,2],method = 'pearson')
cor(dat3[,1],dat3[,2],method = 'pearson')

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

基本同上

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

options(stringsAsFactors = F)
library(airway)
data(airway)
airway 
RNAseq_expr=assay(airway)
colnames(RNAseq_expr) 
RNAseq_expr[1:4,1:4] #行是基因，列是样本
# RNAseq_expr 是一个数值型矩阵，属于连续性变量，可以探索众数、分位数和平均数 ，极差，方差和标准差等统计学指标
RNAseq_gl=colData(airway)[,3]
table(RNAseq_gl)#分组信息
cor(RNAseq_expr)

是 8个样本的RNA-seq数据的counts矩阵，这8个样本分成2组，每组是4个样本， 分别是 trt 和 untrt 组。
通过上面的代码，我们得到了对airway数据集的RNA-seq数据的counts矩阵，命名为 RNAseq_expr 下面会用得到。

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

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

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

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检验

max=which.max(rowSums(RNAseq_expr))
t.test(RNAseq_expr[4464,]~RNAseq_gl)
pos

不显著

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

pos=max=which.max(rowMads(RNAseq_expr))
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)

525eeac375d0bb153a4b15a4f5fb278.png

看看是不是基因变化了，统计结果也变化了:是的

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

cg=names(tail(sort(apply(RNAseq_expr,1,mad)),100))
#这里一定要用apply，apply得到的结果会加上name，其余没这么方便
dat=RNAseq_expr[cg,]
plot(hclust(dist(t(dat))))
colnames(dat)
RNAseq_gl
plot(hclust(dist( dat )))

Q9: 取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 )))

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

cg=names(tail(sort(apply(RNAseq_expr,1,mad)),100))
dat=RNAseq_expr[cg,]
pheatmap::pheatmap(t(scale(t(dat))))

统计检验相关

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

options(stringsAsFactors = F)
rm(list=ls())
library(airway) 
RNAseq_expr=assay(airway)
e1=RNAseq_expr[apply(RNAseq_expr,1,function(x) sum(x>0)>1),] #这里就是过滤
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矩阵剔除掉

t.test算不出来是因为标准差两个都是0

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,]#|同时满足在ab中

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

可以算出来了

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

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

e2=log(e1+1)

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

p1=apply(e1, 1, function(x){
 t.test(x~RNAseq_gl)$p.value
}) 
p2=apply(e2, 1, function(x){
 t.test(x~RNAseq_gl)$p.value
}) 
plot(p1,p2)
cor(p1,p2)