统计学基本概念
四分位数（Quartile)
统计学中，把所有数值由小到大排列并分成四等份，处于三个分割点位置的数值就是四分位数。第一，四分位数 (Q1)，又称“较小四分位数”，等于该样本中所有数值由小到大排列后第25%的数字。第二，四分位数 (Q2)，又称“中位数”，等于该样本中所有数值由小到大排列后第50%的数字。第三，四分位数 (Q3)，又称“较大四分位数”，等于该样本中所有数值由小到大排列后第75%的数字。第三，四分位数与第一四分位数的差距又称四分位距（InterQuartile Range,IQR）。四分位距的优点在于，与全距（极差）相比，较少受异常值的影响。 R中用来显示四分位数的函数是quantile，另外用boxplot可以绘制出某个数据集的箱线图。
箱线图结果boxplot(w)如下所示：
image.png
线图解读：最下面是下界，最上面的圆圈是上界，上界的圆圈是异常值，中间矩形的底边是下分位数，上边是上四分位数，中间的粗线是中位数，箱体的高是四分位距。
方差
image.png
标准差
是方差的正平方根，其单位与原变量值的单位相同
image.png
p值
P 值即概率，反映某一事件发生的可能性大小。统计学根据显著性检验方法所得到的P 值，一般以P < 0.05 为有统计学差异， P<0.01 为有显著统计学差异，P<0.001为有极其显著的统计学差异。其含义是样本间的差异由抽样误差所致的概率小于0.05 、0.01、0.001。
数据解释
P>0.05
碰巧出现的可能性大于5%
不能否定无效假设
两组差别无显著意义
P<0.05
碰巧出现的可能性小于5%
可以否定无效假设
两组差别有显著意义
P <0.01
碰巧出现的可能性小于1%
可以否定无效假设
两者差别有非常显著意义
p值的意义
P的意义不表示两组差别的大小，P反映两组差别有无统计学意义，并不表示差别大小
t检验
T检验，亦称student t检验（Student's t test），主要用于样本含量较小（例如n < 30），总体标准差σ未知的正态分布。 [1] T检验是用t分布理论来推论差异发生的概率，从而比较两个平均数的差异是否显著。
GEO分析流程
image.png
数据下载
> rm(list = ls())
> options(stringsAsFactors = F)
> library(GEOquery)
> eSet <- getGEO("GSE78179",
+ destdir = '.',
+ getGPL = F)###eSet是一个对象（比列表更高级）
#（1）从eSet中提取表达矩阵exp
> eSet[[1]]@annotation###获取平台信息
[1] "GPL13497"
> #（1）从eSet中提取表达矩阵exp
> class(eSet)
[1] "list"
> length(eSet)
[1] 1
> class(eSet[[1]])
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
image.png
#(2)提取临床信息,用pData函数
> pd <- pData(eSet[[1]])
#保存文件
> save(pd,exp,file = "step1output.Rdata")
获得分组信息
> rm(list = ls()) ## 魔幻操作，一键清空~
> options(stringsAsFactors = F)
> load(file = "step1output.Rdata")
#此处需要耗时端详这个pd,找到里面的分组信息，生成一个grouplist向量。
> class(pd)
[1] "data.frame"
> dim(pd)
[1] 6 39
> colnames(pd)
[1] "title" "geo_accession" "status"
[4] "submission_date" "last_update_date" "type"
[7] "channel_count" "source_name_ch1" "organism_ch1"
[10] "characteristics_ch1" "characteristics_ch1.1" "characteristics_ch1.2"
[13] "treatment_protocol_ch1" "growth_protocol_ch1" "molecule_ch1"
[16] "extract_protocol_ch1" "label_ch1" "label_protocol_ch1"
[19] "taxid_ch1" "hyb_protocol" "scan_protocol"
[22] "description" "data_processing" "platform_id"
[25] "contact_name" "contact_email" "contact_phone"
[28] "contact_laboratory" "contact_department" "contact_institute"
[31] "contact_address" "contact_city" "contact_zip/postal_code"
[34] "contact_country" "supplementary_file" "data_row_count"
[37] "age:ch1" "cell line:ch1" "cell type:ch1"
#找到title列，查看其分组信息
> group_list=c(rep("carcinoma",times=3),rep("health",times=3))###对其进行分组
> group_list
[1] "carcinoma" "carcinoma" "carcinoma" "health" "health" "health"
> save(group_list,file = "step2output.Rdata")
检查数据（检查分组之间是否有差异）
> rm(list = ls()) ## 魔幻操作，一键清空~
> options(stringsAsFactors = F)
> load(file = "step1output.Rdata")
> load(file = "step2output.Rdata")
# 看看数据
> exp[1:4,1:4]
GSM2068840 GSM2068841 GSM2068842 GSM2068843
(+)E1A_r60_1 -0.10223055 -0.02216482 -0.04937220 0.02216482
(+)E1A_r60_3 -0.06468773 0.12406254 -0.41914177 0.01200771
(+)E1A_r60_a104 -0.10492134 -0.09622002 -0.27500343 0.15058613
(+)E1A_r60_a107 -0.32794570 -0.03319979 -0.02773619 0.06135511
#exp = log2(exp+1)
> dat=as.data.frame(t(exp))#转置
> dim(exp)
[1] 34183 6
> dim(dat)
[1] 6 34183
> dat[1:4,1:4]
(+)E1A_r60_1 (+)E1A_r60_3 (+)E1A_r60_a104 (+)E1A_r60_a107
GSM2068840 -0.10223055 -0.06468773 -0.10492134 -0.32794570
GSM2068841 -0.02216482 0.12406254 -0.09622002 -0.03319979
GSM2068842 -0.04937220 -0.41914177 -0.27500343 -0.02773619
GSM2068843 0.02216482 0.01200771 0.15058613 0.06135511
> dat=cbind(dat,group_list)##将group_list加入
# pca的统一操作走起
library(FactoMineR)#画主成分分析图需要加载这两个包
library(factoextra)
> dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)###PCA函数：对补充个体、补充数量变量和补充类别变量进行主成分分析(PCA)。
> fviz_pca_ind(dat.pca,
+ geom.ind = "point", # show points only (nbut not "text") ###加入点
+ col.ind = dat$group_list, # color by groups ###用group_list进行颜色划分
+ palette = c("#00AFBB", "#E7B800"),
+ addEllipses = TRUE, # Concentration ellipses
+ legend.title = "Groups"###名字
+ )
> ggsave('all_samples_PCA.png')
PCA
分组之间能分得开，则该数据可以使用；如果重合严重，则不可用。
热图
> cg=names(tail(sort(apply(exp,1,sd)),1000))###选取方差最大的1000个基因（方差大，则差异显著）
> library(pheatmap)
> pheatmap(exp[cg,],
+ show_colnames = F,
+ show_rownames = F)
image.png
归一化处理
> n=t(scale(t(exp[cg,]))) # 'scale'可以对log-ratio数值进行归一化
归一化方法有两种形式，一种是把数变为（0，1）之间的小数，一种是把有量纲表达式变为无量纲表达式。主要是为了数据处理方便提出来的，把数据映射到0～1范围之内处理，更加便捷快速，应该归到数字信号处理范畴之内。
例子
把数变为（0，1）之间的小数
例1：{2.5 3.5 0.5 1.5}归一化后变成了{0.3125 0.4375 0.0625 0.1875}
image
解：2.5+3.5+0.5+1.5=8，
2.5/8=0.3125，
3.5/8=0.4375，
0.5/8=0.0625，
1.5/8=0.1875.
这个归一化就是将括号里面的总和变成1.然后写出每个数的比例。
设置阈值
> thr=2
> n[n>thr]=thr
> n[n< -thr]= -thr
> n[1:4,1:4]
GSM2068840 GSM2068841 GSM2068842 GSM2068843
A_23_P212800 -0.9920418 -0.8542348 -0.8871113 0.9854177
A_33_P3333777 0.8698776 0.8780263 0.9860411 -0.8442771
A_32_P84728 -0.9628674 -0.8859174 -0.8866393 0.9828763
A_23_P369654 0.7580830 1.0977170 0.8491186 -0.7066642
> pheatmap(n,show_colnames =F,show_rownames = F)
image.png
加注释分组
> ac=data.frame(group=group_list)
> rownames(ac)=colnames(n)
> pheatmap(n,show_colnames =F,show_rownames = F,
+ annotation_col=ac)
image.png
经过归一化处理的图，差异更加明显；不经归一化处理，容易将一些稍小些的波动覆盖掉。其实很多情况下，我们是希望通过比较基因的变化情况，找到样本中的关键基因，而不是表达量的绝对值。归一化处理后可以将数值范围锁定在一个较小的范围内，可以更加直观的比较同一基因在不同样本的差异情况。
差异分析
> rm(list = ls()) ## 魔幻操作，一键清空~
> options(stringsAsFactors = F)
> #options('scipen'=100, 'digits'=4)
> options(digits = 4)
> load(file = "step1output.Rdata")
> load(file = "step2output.Rdata")
> exp[1:4,1:4]
GSM2068840 GSM2068841 GSM2068842 GSM2068843
(+)E1A_r60_1 -0.10223 -0.02216 -0.04937 0.02216
(+)E1A_r60_3 -0.06469 0.12406 -0.41914 0.01201
(+)E1A_r60_a104 -0.10492 -0.09622 -0.27500 0.15059
(+)E1A_r60_a107 -0.32795 -0.03320 -0.02774 0.06136
> table(group_list)
group_list
carcinoma health
3 3
单个基因差异分析
> boxplot(exp[1,]~group_list)###对第一个基因画箱图，并分组
箱图图解
> bp=function(g){ #定义一个函数g，函数为{}里的内容
+ library(ggpubr)
+ df=data.frame(gene=g,group=group_list)
+ p <- ggboxplot(df, x = "group", y = "gene",
+ color = "group", palette = "jco",
+ add = "jitter")
+ # Add p-value
+ p + stat_compare_means(label.y = 8)
+ }
> bp(exp[1,]) ## 调用上面定义好的函数，避免同样的绘图代码重复多次敲。
> bp(exp[2,])
> bp(exp[3,])
> bp(exp[4,])
image.png
批量差异分析，用limma包来做
需要表达矩阵和group_list，其他都不要动
> library(limma)
> design=model.matrix(~factor(group_list))###design是limma需要的输入值，无需过多理解
> View(design)###?lmFit###拟合线性模型为每个基因给定一系列阵列
?eBayes###给定一个微阵列线性模型拟合，计算慢化t统计量，慢化f统计量，并通过经验贝叶斯调节标准误差到一个共同的值的微分表达式的log-odds。
?topTable###从线性模型拟合中提取顶级基因表。
> fit=lmFit(exp,design)
> fit=eBayes(fit)
差异基因排名
> deg=topTable(fit,coef=2,number = Inf)
> deg[,1] <- -(deg[,1])
> head(deg)
logFC AveExpr t P.Value adj.P.Val B
A_23_P383009 -11.83 -0.023508 145.1 1.179e-11 1.092e-07 15.31
A_23_P83098 10.64 0.002649 -142.0 1.339e-11 1.092e-07 15.27
A_33_P3362008 -8.63 -0.008725 133.8 1.895e-11 1.092e-07 15.16
A_23_P111132 -10.81 -0.027914 124.3 2.915e-11 1.092e-07 15.00
A_23_P419714 -10.07 -0.037299 122.7 3.154e-11 1.092e-07 14.97
A_23_P214011 -10.39 -0.020587 121.8 3.281e-11 1.092e-07 14.96
> bp(exp[rownames(deg)[1],])
> bp(exp[rownames(deg)[2],])
logFC=treat/control的log2值，负值代表基因表达上调，正值代表基因表达下调；
p.value值越小，表示差异表达越大。该数据是按差异从大到小排名的。
image.png
image.png
为deg数据框添加几列
> #1.加probe_id列，把行名变成一列
> library(dplyr)
> deg <- mutate(deg,probe_id=rownames(deg))
> #tibble::rownames_to_column()
> head(deg)
logFC AveExpr t P.Value adj.P.Val B probe_id
1 -11.83 -0.023508 145.1 1.179e-11 1.092e-07 15.31 1
2 10.64 0.002649 -142.0 1.339e-11 1.092e-07 15.27 2
3 -8.63 -0.008725 133.8 1.895e-11 1.092e-07 15.16 3
4 -10.81 -0.027914 124.3 2.915e-11 1.092e-07 15.00 4
5 -10.07 -0.037299 122.7 3.154e-11 1.092e-07 14.97 5
6 -10.39 -0.020587 121.8 3.281e-11 1.092e-07 14.96 6
加symbol列，火山图要用
> #id转换，查找芯片平台对应的包
> #http://www.bio-info-trainee.com/1399.html
> eSet[[1]]@annotation###获取平台信息
[1] "GPL13497"###注释包为HsAgilentDesign026652
> ids <- toTable(hugene10sttranscriptclusterSYMBOL)
> head(ids)
probe_id symbol
1 7896759 LINC01128
2 7896761 SAMD11
3 7896779 KLHL17
4 7896798 PLEKHN1
5 7896817 ISG15
6 7896822 AGRN
网友评论