統(tǒng)計學(xué)是一門很深的學(xué)問赫段，這里僅僅是出題幫助大家熟練使用R語言來學(xué)習(xí)統(tǒng)計學(xué)知識糯笙，具體知識點需要更深入閱讀書籍或者教程：
- 推薦一下 統(tǒng)計學(xué)基礎(chǔ)：[link]https://mp.weixin.qq.com/s/OtB2h6f00U2SRZLzveJKfQ
- 統(tǒng)計學(xué)精華-statQuest教學(xué)視頻：[link]https://mp.weixin.qq.com/s/X0PE9S0BgSuCcAV9zeY1jQ
基礎(chǔ)概念
需要掌握R內(nèi)置數(shù)據(jù)集及R包數(shù)據(jù)集
- 內(nèi)置數(shù)據(jù)集： [link]https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g
- airway 數(shù)據(jù)集 [link] https://mp.weixin.qq.com/s/FTdOyxvlS4L-zL3inHTtZg
理解 定性變量（qualitative variable） 和 定量變量(quantitative variable)
定量數(shù)據(jù)的集中趨勢指標(biāo)主要是：眾數(shù)豺憔、分位數(shù)和平均數(shù)
定量數(shù)據(jù)的離散趨勢指標(biāo)主要是：極差，方差和標(biāo)準(zhǔn)差焕阿，標(biāo)準(zhǔn)分?jǐn)?shù)暮屡，相對離散系數(shù)（變異系數(shù)）褒纲，偏態(tài)系數(shù)與峰態(tài)系數(shù)
Q1: 載入R中自帶的數(shù)據(jù)集 iris钥飞，指出其每列是定性還是定量數(shù)據(jù)
Q2: 對數(shù)據(jù)集 iris的所有定量數(shù)據(jù)列計算集中趨勢指標(biāo)：眾數(shù)读宙、分位數(shù)和平均數(shù)
Q3:對數(shù)據(jù)集 iris的所有定性數(shù)據(jù)列計算水平及頻次
Q4:對數(shù)據(jù)集 iris的所有定量數(shù)據(jù)列計算離散趨勢指標(biāo)：方差和標(biāo)準(zhǔn)差等
Q5:計算數(shù)據(jù)集 iris的前兩列變量的相關(guān)性唇兑，提示cor函數(shù)可以選擇3種methods
Q6:對數(shù)據(jù)集 iris的所有定量數(shù)據(jù)列內(nèi)部zcore標(biāo)準(zhǔn)化桦锄，并計算標(biāo)準(zhǔn)化后每列的平均值和標(biāo)準(zhǔn)差
Q7:計算列內(nèi)部zcore標(biāo)準(zhǔn)化后 iris的前兩列變量的相關(guān)性
Q8: 根據(jù)數(shù)據(jù)集 iris的第五列拆分?jǐn)?shù)據(jù)集后重復(fù)上面的Q2到Q7問題
Q9:載入R中自帶的數(shù)據(jù)集 mtcars结耀，重復(fù)上面的Q1到Q7個問題
Q10: 載入r包airway并且通過assay函數(shù)拿到其表達(dá)矩陣后計算每列之間的相關(guān)性
關(guān)于 airway 代碼如下图甜，需要理解：
options(stringsAsFactors = F)
library(airway)
data(airway)
# 這里需要自行學(xué)習(xí)bioconductor里面的RangedSummarizedExperiment對象
airway
RNAseq_expr=assay(airway)
colnames(RNAseq_expr)
RNAseq_expr[1:4,1:4]
# RNAseq_expr 是一個數(shù)值型矩陣，屬于連續(xù)性變量嚼摩，可以探索眾數(shù)低斋、分位數(shù)和平均數(shù) 膊畴，極差，方差和標(biāo)準(zhǔn)差等統(tǒng)計學(xué)指標(biāo)
RNAseq_gl=colData(airway)[,3]
table(RNAseq_gl)
是 8個樣本的RNA-seq數(shù)據(jù)的counts矩陣病游，這8個樣本分成2組唇跨，每組是4個樣本稠通， 分別是 trt 和 untrt 組。
通過上面的代碼买猖，我們得到了對airway數(shù)據(jù)集的RNA-seq數(shù)據(jù)的counts矩陣改橘，命名為RNAseq_expr下面會用得到。
表達(dá)矩陣相關(guān)
首先了解各種統(tǒng)計分布： [link]https://mp.weixin.qq.com/s/uly4jlQomk9LZlHyknkNdg 在R語言的實現(xiàn)方式玉控。
Q1: 把RNAseq_expr第一列全部加1后取log2后計算平均值和標(biāo)準(zhǔn)差
tmp=log2(RNAseq_expr[,1]+1)
mean(tmp)
sd(tmp)
Q2: 根據(jù)上一步得到平均值和標(biāo)準(zhǔn)差生成同樣個數(shù)的隨機的正態(tài)分布數(shù)值
a=rnorm(length(tmp),mean = mean(tmp),sd = sd(tmp))
a=sort(a)
plot(a)
points(sort(tmp))
Q3: 刪除RNAseq_expr第一列低于5的數(shù)據(jù)后飞主，重復(fù)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的第一列和第二列進(jìn)行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行之和最大的那一行根據(jù)分組矩陣進(jìn)行T檢驗
pos=which.max(rowSums(RNAseq_expr))
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos
Q6: 取RNAseq_expr的MAD最大的那一行根據(jù)分組矩陣進(jìn)行T檢驗
pos=which.max(apply(RNAseq_expr,1,mad))
t.test(RNAseq_expr[pos,]~RNAseq_gl)
pos
Q7: 對RNAseq_expr全部加1后取log2后重復(fù)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)
看看是不是基因變化了，統(tǒng)計結(jié)果也變化了
Q8: 取RNAseq_expr矩陣的MAD最高的100行高诺，對列和行分別進(jìn)行層次聚類
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 )))
檢查一下聚類結(jié)果跟樣本的處理信息是否能對應(yīng)
Q9: 取RNAseq_expr矩陣的SD最高的100行碌识，對列和行分別進(jìn)行層次聚類
Q10: 對Q8矩陣按照行和列分別歸一化并且熱圖可視化
cg=names(tail(sort(apply(RNAseq_expr,1,mad)),100))
dat=RNAseq_expr[cg,]
pheatmap::pheatmap(scale(dat))
pheatmap::pheatmap(t(scale(t(dat))))
統(tǒng)計檢驗相關(guān)
這里需要對前面的RNAseq_expr矩陣進(jìn)行一定程度的過濾，主要是過濾那些每一列都為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)
上面的 e1 矩陣下面的習(xí)題就用得到惠呼。
Q1: 對e1每一行獨立根據(jù)分組矩陣進(jìn)行T檢驗趟畏，檢查為什么有些行T檢驗失敗
apply(e1, 1, function(x){
t.test(x~RNAseq_gl)$p.value
})
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矩陣進(jìn)行每一行獨立根據(jù)分組矩陣進(jìn)行T檢驗
Q4: 對e1矩陣進(jìn)行加1后log2的歸一化命名為e2再對每一行獨立根據(jù)分組矩陣進(jìn)行T檢驗
Q5: 對e1,e2的T檢驗P值做相關(guān)性分析
p1=apply(e1, 1, function(x){
t.test(x~RNAseq_gl)$p.value
})
e2=log(e1+1)
p2=apply(e2, 1, function(x){
t.test(x~RNAseq_gl)$p.value
})
plot(p1,p2)
cor(p1,p2)
寫到這里猎莲，我突然間感覺自己代碼很壯觀，也不知道為什么身笤，就覺得自己寫代碼很帥氣脱篙，后面的5個題目就不寫啦适刀，這些已經(jīng)夠大家用的笔喉。
我這里并沒有提到基因和樣本這樣的詞語吧彪，就是希望其他領(lǐng)域?qū)WR的朋友也可以看看怨酝，如果生物信息學(xué)領(lǐng)域赡艰，這樣的簡單T檢驗是有很多不合理的地方揍堕，比如文庫大小芹血，比如基因表達(dá)量分布等等幔烛。
