> rm(list = ls())
> load(file = "step2output.Rdata")
> #查看右侧变量
>
> #差异分析，用limma包来做
> #需要表达矩阵和Group，不需要改
> library(limma)
Warning message:
程辑包‘limma’是用R版本4.0.3 来建造的
> design=model.matrix(~Group)
> fit=lmFit(exp,design)
> fit=eBayes(fit)#贝叶斯检验
> deg=topTable(fit,coef=2,number = Inf)
> dim(deg)
[1] 54675 6
> head(deg)
logFC AveExpr t P.Value adj.P.Val B
224563_at -3.419608 10.462003 -16.58978 4.947714e-14 2.358899e-09 21.39728
208629_s_at -1.633518 11.204696 -16.14895 8.628803e-14 2.358899e-09 20.92048
212417_at -1.810946 8.760968 -15.16339 3.139295e-13 4.882577e-09 19.79803
201861_s_at -1.713884 11.965832 -15.06767 3.572073e-13 4.882577e-09 19.68466
208325_s_at -1.840898 9.321510 -13.60230 2.819698e-12 3.083340e-08 17.84515
203479_s_at -2.274611 7.651265 -13.43343 3.618635e-12 3.297481e-08 17.61994
> identical(rownames(deg),rownames(exp))#已按照p值列排序 不同了
[1] FALSE
> head(rownames(deg))
[1] "224563_at" "208629_s_at" "212417_at" "201861_s_at" "208325_s_at" "203479_s_at"
> head(rownames(exp))
[1] "1007_s_at" "1053_at" "117_at" "121_at" "1255_g_at" "1294_at"
>
> #为deg数据框添加几列
> #1.加probe_id列，把行名变成一列
> library(dplyr)
载入程辑包：‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Warning message:
程辑包‘dplyr’是用R版本4.0.3 来建造的
> deg <- mutate(deg,probe_id=rownames(deg))
> head(deg)
logFC AveExpr t P.Value adj.P.Val B probe_id
1 -3.419608 10.462003 -16.58978 4.947714e-14 2.358899e-09 21.39728 224563_at
2 -1.633518 11.204696 -16.14895 8.628803e-14 2.358899e-09 20.92048 208629_s_at
3 -1.810946 8.760968 -15.16339 3.139295e-13 4.882577e-09 19.79803 212417_at
4 -1.713884 11.965832 -15.06767 3.572073e-13 4.882577e-09 19.68466 201861_s_at
5 -1.840898 9.321510 -13.60230 2.819698e-12 3.083340e-08 17.84515 208325_s_at
6 -2.274611 7.651265 -13.43343 3.618635e-12 3.297481e-08 17.61994 203479_s_at
> #2.加上探针注释
> #查看
> dim(ids)
[1] 41905 2
> head(ids)
probe_id symbol
1 1053_at RFC2
2 117_at HSPA6
3 121_at PAX8
4 1255_g_at GUCA1A
5 1316_at THRA
6 1320_at PTPN21
> table(!duplicated(ids$probe_id))#无重复的 第一列probe_id
TRUE
41905
> table(!duplicated(ids$symbol))#无重复的 第二列 symbol
FALSE TRUE
21744 20161
去重
按symbol列去重，常见标准有3个：最大值/平均值/随机去重
方法1.随机去重
> ids = ids[!duplicated(ids$symbol),]#谁排第一个就取谁
>
> deg <- inner_join(deg,ids,by="probe_id")#去交集
> head(deg)
logFC AveExpr t P.Value adj.P.Val B probe_id symbol
1 -1.633518 11.204696 -16.14895 8.628803e-14 2.358899e-09 20.92048 208629_s_at HADHA
2 -1.713884 11.965832 -15.06767 3.572073e-13 4.882577e-09 19.68466 201861_s_at LRRFIP1
3 -2.274611 7.651265 -13.43343 3.618635e-12 3.297481e-08 17.61994 203479_s_at OTUD4
4 -1.935657 8.628444 -13.04331 6.501205e-12 3.554534e-08 17.08858 203391_at FKBP2
5 -1.893808 8.896532 -12.38173 1.812689e-11 6.607251e-08 16.15073 205457_at ILRUN
6 -1.691681 10.302132 -12.00890 3.290947e-11 1.124578e-07 15.60094 201392_s_at IGF2R
> nrow(deg)
[1] 20161
方法2.求平均值
> #rm(list = ls())
> #load("step2output.Rdata")
> exp3 = exp[ids$probe_id,]
> exp3[1:4,1:4]
GSM1366348 GSM1366349 GSM1366350 GSM1366351
1053_at 8.932805 8.679543 8.625015 8.637085
117_at 9.383421 8.605809 9.462774 9.898573
121_at 7.916751 8.500635 8.258467 8.553656
1255_g_at 5.085221 2.414033 1.718570 4.311794
> # 有重复的基因,求平均值
> s = unique(ids$symbol[duplicated(ids$symbol)]);length(s)
[1] 10753
> exp4 = sapply(s, function(i){
+ colMeans(exp3[ids$symbol==i,])
+ })
> dim(exp4)
[1] 22 10753
> exp4 = t(exp4)
> dim(exp4)
[1] 10753 22
