R語言作業(yè)-統(tǒng)計30題

題目

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

我做題的時候主要翻閱學習了《R語言實戰(zhàn)》里統(tǒng)計相關內容。

基礎概念

需要掌握R內置數據集及R包數據集

定性/定量數據

鳶尾花(iris)數據集,包含150個鳶尾花的信息,共五列,分別為萼片長度(Sepal.Length)、萼片寬度(Sepal.Width)、花瓣長度(Petal.Length)、花瓣寬度(Petal.Width)和種類(Species)。前四列為定量數據,后一列種類為定性數據,是非連續(xù)的字符變量。

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


描述性數據

定量數據的集中趨勢指標主要是:眾數、分位數和平均數
定量數據的離散趨勢指標主要是:極差,方差和標準差,標準分數,相對離散系數(變異系數),偏態(tài)系數與峰態(tài)系數

# 平均數
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  



分組統(tǒng)計,最一開始是想將數據集分成三個數據框,重復之前函數,有幾個方法:

# 方法一: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()函數分開一個個計算統(tǒng)計數
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函數可以解決數據循環(huán)處理的問題,可以對矩陣、數據框、數組(二維、多維),按行或列進行循環(huán)計算,對子元素進行迭代,并把子元素以參數形式給自定義的FUN函數中,并返回計算結果。

函數定義:
apply(X,MARGIN,FUN,...)
參數列表:

  • X:數組、矩陣、數據框
  • MARGIN:1表示按行,2表示按列
  • FUN:自定義調用函數
  • ...
lapply函數

用來對list、data.frame進行循環(huán),并返回和X長度同樣的list結構作為結果集。

sapply函數

同lapply函數,多了2個參數simplify和USE.NAMES,返回值為向量,不是list對象。

vapply函數

類似sapply函數,提供了FUN.VALUE參數,用來控制返回值的行名。

mapply函數

類似sapply函數,第一個參數為FUN,可接受多個數據。

tapply函數

tapply函數用于分組的循環(huán)計算,相當于group by的操作。

函數定義:
tapply(X,INDEX,FUN,simplify,...)

參數列表:

  • X:向量
  • INDEX:用于分組的索引
  • FUN:自定義的調用函數
  • simplify:是否數組化
rapply函數

只處理list類型數據,對list的每個元素進行遞歸遍歷,如果list包括子元素則繼續(xù)遍歷。


相關性

R可以計算多種相關系數,包括Pearson相關系數、Spearman相關系數、Kendall相關系數、偏相關系數、多分格相關系數、多系列相關系數。cor()函數可以計算前三種相關系數,cov()函數可以計算協(xié)方差。

# 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列分組,統(tǒng)計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個樣本,分別是 trtuntrt 組。

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 是一個數值型矩陣,屬于連續(xù)性變量,可以探索眾數、分位數和平均數 ,極差,方差和標準差等統(tǒng)計學指標

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

# 根據上一步得到平均值和標準差生成同樣個數的隨機的正態(tài)分布數值
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統(tǒng)計值;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ǒng)計特征(比如上面研究的平均值、眾數、方差等)做判斷的方法為假設檢驗。

獨立樣本t檢驗

一個針對兩組的獨立樣本t檢驗可以用于檢驗兩個總體的均值相等的假設,檢驗調用格式為:
t.test( y ~ x, data )
其中y是一個數值型變量,x是一個二分變量。

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

現在還不能用自己的語言解釋清楚,整合幾篇寫的比較詳細的教程:
http://www.biye5u.com/article/R/2019/6399.html
http://www.itdecent.cn/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

統(tǒng)計這部分還是似懂非懂的狀態(tài),繼續(xù)學習~

更多學習資源:
生信技能樹公益視頻合輯
生信技能樹簡書賬號
生信工程師入門最佳指南
生信技能樹全球公益巡講
招學徒
...
你的宣傳能讓數以萬計的初學者找到他們的家,技能樹平臺一定不會辜負每一個熱愛學習和分享的同道中人 ??

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

友情鏈接更多精彩內容