美文网首页rna_seq
基于R的统计化30题

基于R的统计化30题

作者: 晓颖_9b6f | 来源:发表于2020-10-22 09:53 被阅读0次

    一.基础

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

    data(iris)
head(iris)   ##Species是定性数据,其他的都是定量数据.

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

    mode <- function(x) {
  uniqv <- unique(x)
  uniqv[which.max(tabulate(match(x, uniqv)))]
}
apply(iris,2,mode)   #可以对文字求众数
summary(iris)         
apply(iris[1:4],2,mean)

Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
       "5.0"        "3.0"        "1.4"        "0.2"     "setosa"

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

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

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

    apply(iris[1:4],2,sd)
apply(iris[1:4],2,var)
range <-function(x){
  max(x)-min(x)
 }
apply(iris[1:4],2,range)
library(raster)
install.packages('raster')
apply(iris[1:4],2,raster::cv)

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

    > cor(iris$Sepal.Length,iris$Sepal.Width,method = "pearson")
[1] -0.1175698
> cor(iris$Sepal.Length,iris$Sepal.Width,method = "kendall")
[1] -0.07699679
> cor(iris$Sepal.Length,iris$Sepal.Width,method = "spearman")
[1] -0.1667777
##这三种方法的区别???仔细看看

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

    iris <- apply(iris[1:4],2,scale)
> apply(iris,2,mean)
 Sepal.Length   Sepal.Width  Petal.Length   Petal.Width 
-4.484318e-16  2.034094e-16 -2.895326e-17 -3.663049e-17 
> apply(iris,2,sd)
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
           1            1            1            1 
##标准化是有好几种方法的.zcore标准化好像是针对平均值和方差进行标准化的.

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

    iris <- as.data.frame(iris)
cor(iris$Sepal.Length,iris$Sepal.Width,method = "pearson")
cor(iris$Sepal.Length,iris$Sepal.Width,method = "kendall")
cor(iris$Sepal.Length,iris$Sepal.Width,method = "spearman")

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

    data("iris")
head(iris)
unique(iris$Species)
setosa <- iris[iris$Species=="setosa",]
head(setosa)
versicolor <- iris[iris$Species=="versicolor",]
virginica <- iris[iris$Species=="virginica",]

all_calculate <- function(data){
  mode <- apply(data[,1:4],2,mode)#众数
  fivenum <- summary(data)#分位数
  mean <- apply(data[,1:4],2,mean)#平均数
  frequence <- table(data$Species)
  sd <- apply(data[,1:4],2,sd)#标准差
  var <- apply(data[,1:4],2,var)#方差
  cor <- cor(data[,1],data[,2])
  
  data_scale=t(scale(t(data[,1:4])))
  mean_scale <- apply(data_scale,2,mean)#平均数
  sd_scale <- apply(data_scale,2,sd)#标准差
  cor_scale <- cor(data_scale[,1],data_scale[,2])
  
  return(c(mean=mean, fivenum=fivenum, sd=sd, var=var, frequence=frequence,
           mode=mode, cor=cor, mean_scale=mean_scale, sd_scale=sd_scale,
           cor_scale=cor_scale))
}

all_calculate(setosa)

all_calculate(versicolor)

all_calculate(virginica)

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

    data("mtcars")
all_calculate <- function(data){
  head(data)
  mode <- apply(data[,1:4],2,mode)#众数
  fivenum <- summary(data)#分位数
  mean <- apply(data[,1:4],2,mean)#平均数
  frequence <- table(data$Species)
  sd <- apply(data[,1:4],2,sd)#标准差

  var <- apply(data[,1:4],2,var)#方差
  cor <- cor(data[,1],data[,2])
  
  data_scale=t(scale(t(data[,1:4])))
  mean_scale <- apply(data_scale,2,mean)#平均数
  sd_scale <- apply(data_scale,2,sd)#标准差
  cor_scale <- cor(data_scale[,1],data_scale[,2])
  
  return(c(mean=mean, fivenum=fivenum, sd=sd, var=var, frequence=frequence,
           mode=mode, cor=cor, mean_scale=mean_scale, sd_scale=sd_scale,
           cor_scale=cor_scale))
}

all_calculate(mtcars)

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

    library(airway)
data(airway)
RNAseq_expr=assay(airway)

M <- cor(RNAseq_expr)
pheatmap::pheatmap(M)

    二.表达矩阵相关

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

    head(RNAseq_expr)
RNAseq_gl=colData(airway)[,3]  # ???
print(RNAseq_gl)
table(RNAseq_gl)
class(RNAseq_expr)
log <- log(RNAseq_expr[,1]+1)
head(log)
mean(log)
sd(log)

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

    length(log)
rnorm <- rnorm(64102,mean = 1.560774,sd=2.527118)

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

    RNAseq_expr <- RNAseq_expr[RNAseq_expr[,1]>=5,]
log <- log2(RNAseq_expr[,1]+1)
tmp <- RNAseq_expr[,1]
rnorm <-rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))

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

    head(RNAseq_expr)
RNAseq_expr <- RNAseq_expr[RNAseq_expr[,1]>=5,]
RNAseq_expr <- RNAseq_expr[RNAseq_expr[,2]>=5,]
t.test(RNAseq_expr[,1],RNAseq_expr[,2])

library(ggpubr)
library(ggplot2)

t <- data.frame(value=c(RNAseq_expr[,1],RNAseq_expr[,2]))
t$group <- c(rep('SRR1039508',length(RNAseq_expr[,1])),rep('SRR1039509',length(RNAseq_expr[,2])))
ggboxplot(t, y = "value", x = "group")

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

    RNAseq_expr <- assay(airway)
max <- which.max(apply(RNAseq_expr,1,sum))
RNAseq_expr[max,]

RNAseq_gl
class(RNAseq_gl)
t.test(RNAseq_expr[max,]~RNAseq_gl)

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

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

    Q7: 对RNAseq_expr全部加1后取log2后重复Q5和Q6

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

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

    a <- apply(RNAseq_expr,1,mad)
high100 <- names(tail(sort(a),100))
high100 <- RNAseq_expr[high100,]
dim(high100)
head(high100)
plot(hclust(dist(high100)))
high100 <- t(high100)
plot(hclust(dist(high100)))

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

    a <- apply(RNAseq_expr,1,sd)
high100 <- names(tail(sort(a),100))
high100 <- RNAseq_expr[high100,]
dim(high100)
head(high100)
plot(hclust(dist(high100)))
high100 <- t(high100)
plot(hclust(dist(high100)))

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

    pheatmap::pheatmap(scale(high100))
pheatmap::pheatmap(t(scale(t(high100))))

    三.统计检验相关

    先过滤,过滤那些每一列都为0的行。

    rm(list = ls())
library(airway)
data(airway)
RNAseq_expr <- assay(airway)
a <- RNAseq_expr[apply(RNAseq_expr,1,function(x) sum(x>0)>1),] 
dim(a)
dim(RNAseq_expr)
RNAseq_gl <- colData(airway)[,3]
table(RNAseq_gl)

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

    head(a)
t <- function(x){
  t.test(x~RNAseq_gl)$P.value
}
apply(a,1,t)
##问题在于t.test无法比较各项相等的两个向量

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

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

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

      e2=log2(e1+1)
  apply(e2,1,function(x) t.test(x~RNAseq_gl)$p.value)

    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)
  cor(p1,p2)

    相关文章

      网友评论

        本文标题：基于R的统计化30题

        本文链接：https://www.haomeiwen.com/subject/ijhsmktx.html

        延伸阅读

        深度阅读

        • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

        栏目导航

        热点阅读

        rna_seq

        关于我们|服务条款|联系我们|基于R的统计化30题|投稿指南|网站地图|RSS订阅|排版工具|手机版

        提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

        Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

        备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

        本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

        所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！