学习内容（WCGNA）
之前在GEO数据挖掘的时候出现了问题，所以就去跑了一遍，又熟悉了一下流程，刚好现在来学习一下TCGA
基本概念
WGCNA其译为加权基因共表达网络分析。该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。
适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答
基本原理
从方法上来讲，WGCNA分为表达量聚类分析和表型关联两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。
第一步计算任意两个基因之间的相关系数（Person Coefficient）。为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。但是这样如果将阈值设为0.8，那么很难说明0.8和0.79两个是有显著差别的。因此，WGCNA分析时采用相关系数加权值，即对基因相关系数取N次幂，使得网络中的基因之间的连接服从无尺度网络分布(scale-freenetworks)，这种算法更具生物学意义。
第二步通过基因之间的相关系数构建分层聚类树，聚类树的不同分支代表不同的基因模块，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块。这样就可以将几万个基因通过基因表达模式被分成了几十个模块，是一个提取归纳信息的过程。
直接上代码,用了健明老师的数据先认识一下流程及做什么
load('GSE48213-wgcna-input.RData')
if(T){
fpkm[1:4,1:4]
head(datTraits)
table(datTraits$subtype)
RNAseq_voom <- fpkm
## 因为WGCNA针对的是基因进行聚类，而一般我们的聚类是针对样本用hclust即可，所以这个时候需要转置
WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:5000],])
datExpr0 <- WGCNA_matrix ## top 5000 mad genes
datExpr <- datExpr0
## 下面主要是为了防止临床表型与样本名字对不上
sampleNames = rownames(datExpr);
traitRows = match(sampleNames, datTraits$gsm)
rownames(datTraits) = datTraits[traitRows, 1]
}
## step 2
datExpr[1:4,1:4]
if(T){
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
#设置网络构建参数选择范围，计算无尺度分布拓扑矩阵
png("step2-beta-value.png",width = 800,height = 600)
# Plot the results:
##sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
dev.off()
}
感觉这个WCGNA 自己按照流程期间有点问题，下午的时候跟另一位学徒再跑关于GEO数据挖掘，把之前一点点小尾巴给解决掉。
rm(list = ls())
.libPaths(c("C:/Users/22349/R3.6.1", .libPaths()))
.libPaths()
install.packages("devtools")
install.packages("rlang")
install_github("jmzeng1314/idmap")
library(idmap)
library(AnnoProbe)
library(ggpubr)
getwd()
install.packages("missMDA")
require(missMDA)
suppressPackageStartupMessages(library(GEOquery))
gset=AnnoProbe::geoChina('GSE41804')
eset=gset[[1]]
dat <- exprs(gset[[1]]) #（提取表达矩阵）
dat[1:4,1:4]
pd <- pData(gset[[1]]) # (提取临床信息)
pd[1:4,1:4]
boxplot(dat,las=2)
############代码查看矩阵是否需要log（来自GEO2R）
qx <- as.numeric(quantile(probes_expr, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0) ||
(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
LogC
#dat=log2(dat+1)
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat,las=2)
#分组
pd1=pd[,apply(pd, 2, function(x){
length(unique(x))>1})] #缩小范围，从所有临床信息中筛选出含有大于2个元素的信息列
dim(pd1)
apply(pd1,2,table)
table(pd1$`tissue:ch1`)
HCC=rownames(pd1[grepl('resected liver tumor',as.character(pd$`tissue:ch1`)),])
Normal=rownames(pd1[grepl('resected non-tumor liver tissue',as.character(pd$`tissue:ch1`)),])
dat=dat[,c(HCC,Normal)]
dim(dat)
group_list=c(rep('HCC',length(HCC)) ,
rep('Normal',length(Normal))) #分组信息
table(group_list)
eset@annotation
GPL=eset@annotation
save(dat,eset,pd,gset,GPL,group_list,file='GSE41804.Rdata')
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'GSE41804.Rdata')
# 每次都要检测数据
dat[1:4,1:4]
## 下面是画PCA的必须操作，需要看说明书。
dat=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
dat=as.data.frame(dat)#将matrix转换为data.frame
dat=cbind(dat,group_list) #cbind横向追加，即将分组信息追加到最后一列
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
# The variable group_list (index = 54676) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = dat$group_list, # color by groups
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
#id注释转换
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
library(AnnoProbe)
load('GSE41804.Rdata')
# 选择 pipe 获取的是 冗余注释，也就是说一个探针很有可能会对应多个基因。
ids<- idmap(GPL,type = 'soft')
head(ids)
colnames(ids)=c('probe_id','symbol')
ids=ids[ids$probe_id %in% rownames(dat),]
ids$probe_id=as.character(ids$probe_id)
dat[1:4,1:4]
dat=dat[ids$probe_id,]
ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4] #保留每个基因ID第一次出现的信息
dim(dat)
library(clusterProfiler)
#DEG
group_list <- factor(group_list)
class(group_list)
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(dat)
library(limma)
####构建比较矩阵,设置用来对比的组别
contrast.matrix<-makeContrasts(contrasts=
c('HCC-Normal'),levels = design)
######limma三部曲,只需要归一化后的数据、实验矩阵、比较矩阵的输入
fit <- lmFit(dat,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
####上面得到了p值等统计的结果，topTable对p值校验，对基因排序
tT <- topTable(fit2, adjust="fdr", number=nrow(fit2))
dim(tT)
####提取出我们需要的指标
tT <- subset(tT, select=c('adj.P.Val',"P.Value","logFC"))
save(tT,file='DEG.Rdata')
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
load('tT1DEG.Rdata')
b <- rownames(tT)
gene<- bitr(b, fromType = "SYMBOL", #fromType是指你的数据ID类型是属于哪一类的
toType = c('ENTREZID'), #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
OrgDb = org.Hs.eg.db)#Orgdb是指对应的注释包是哪个
tT1$SYMBOL <- tT$symbol
tT2 <- merge(gene,tT)
d <- tT2
## assume that 1st column is ID
## 2nd column is fold change
## feature 1: numeric vector
geneList1 <- d[,6]
## feature 2: named vector
names(geneList1) <- as.character(d[,2])
## feature 3: decreasing order
geneList1 <- sort(geneList1, decreasing = TRUE)
save(geneList1,file='geneList1.Rdata')
head(geneList1)
这上面就是做的gene list 因为之前做过 就直接一串代码下来，然后流程跑
image.png
image.png
image.png
image.png
