参考文章：
本文主要讲illumina beadchip 芯片原始数据处理，基于数据GSE73461
- 原始表达矩阵存在负值（无法进行下游分析），所以必须自己处理原始数据
> x <- read.ilmn(files="GSE73461_Raw_data.txt",
+ expr="Sample",
+ probeid="ID_REF",
+ other.columns="Detection Pval")
Reading file GSE73461_Raw_data.txt ... ...
Error in readGenericHeader(fname, columns = expr, sep = sep) :
Specified column headings not found in file
- 上一步现实原始数据，列名不统一，read.ilmn函数无法识别；所以就需要修改
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)#
# 读取原始数据
raw_data=read.table('GSE73461_series_raw_matrix.txt.gz',sep='\t',quote = "",fill = T,
comment.char = "!",header = T) # 提取表达矩阵
raw_data=raw_data[,-2]
a=raw_data
## 修改列明是为了统一名称，方便后续read.ilmn能够识别expr="Sample",other.columns="Detection Pval"
colnames(a)[seq(2,919,by=2)]=paste0("Sample",seq=".",1:459)
colnames(a)[seq(3,919,by=2)]=paste0("Detection Pval",seq=".",1:459)
#再次导出数据为TXT文本，方便read.ilmn读入
write.table(a,"limma.txt",sep='\t',quote = FALSE) #quote = FALSE取消列明加用引号
原始数据处理
#1.读取、背景校正和标准化
x <- read.ilmn(files="limma.txt",
expr="Sample",
probeid="ID_REF",
other.columns="Detection Pval")
## Reading file GSE16997_raw.txt ... ...
y <- neqc(x,detection.p="Detection Pval")
#2.探针过滤
x$other$Detection[1:4,1:4]
dim(y)
expressed <- rowSums(y$other$`Detection Pval` < 0.05) >= 115 ;table(expressed)
y <- y[expressed,]
dim(y)
#3.提取表达矩阵
exp = as.data.frame(y$E)
一些概念
1.关键基因筛选
FilterGenes= abs(GS1)> .2 & abs(datKME$MM.brown)>.8
- WGCNA如何从module中挖掘关键基因
- gene significance：对于基因定义了GS值，表征基因和表型之间的相关性。简称GS, 将该基因的表达量与对应的表型数值进行相关性分析，最终的相关系数的值就是GS, GS反映出基因表达量与表型数据的相关性，计算GS的代码如下
- module membership：1.Hub基因筛选需要确定模块成员(module membership，MM)，以确定基因与给定模块之间的相关性。2.简称MM, 将该基因的表达量与module的第一主成分，即module eigengene进行相关性分析就可以得到MM值，所以MM值本质上是一个相关系数，如果基因和某个module的MM值为0，说明二者根本不相关，该基因不属于这个module; 如果MM的绝对值接近1，说明基因与该module相关性很高。
2.样本聚类分析
- 使用dist()函数计算样品之间的欧式距离，并通过聚类查看是否存在异常样品。
sampleTree = hclust(dist(datExpr0), method = "average")
3.GO/KEGG
- 3分钟了解GO/KEGG功能富集分析
- 并用clusterProfile包对Hub基因在京都基因与基因组百科全书（Kyoto Encyclope⁃dia of Genes and Genomes，KEGG）中进行通路富集分析，在基因本体论（gene Ontology，GO）数据库中进行基因富集分析。
4.排序并取方差前25%的基因筛选
exp <- as.data.frame(dat)
exp$sd=apply(exp,1,sd)
exp=exp[order(exp$sd,decreasing = TRUE),]
dim(exp)
tj <- quantile(exp$sd,0.75);tj
exp <- exp[exp$sd >= tj,]
dim(exp)
