學(xué)習(xí)內(nèi)容（WCGNA）
之前在GEO數(shù)據(jù)挖掘的時候出現(xiàn)了問題，所以就去跑了一遍，又熟悉了一下流程，剛好現(xiàn)在來學(xué)習(xí)一下TCGA
基本概念
WGCNA其譯為加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析。該分析方法旨在尋找協(xié)同表達(dá)的基因模塊(module)书聚，并探索基因網(wǎng)絡(luò)與關(guān)注的表型之間的關(guān)聯(lián)關(guān)系，以及網(wǎng)絡(luò)中的核心基因藻雌。
適用于復(fù)雜的數(shù)據(jù)模式雌续，推薦5組(或者15個樣品)以上的數(shù)據(jù)。一般可應(yīng)用的研究方向有：不同器官或組織類型發(fā)育調(diào)控胯杭、同一組織不同發(fā)育調(diào)控驯杜、非生物脅迫不同時間點(diǎn)應(yīng)答、病原菌侵染后不同時間點(diǎn)應(yīng)答
基本原理
從方法上來講做个，WGCNA分為表達(dá)量聚類分析和表型關(guān)聯(lián)兩部分鸽心，主要包括基因之間相關(guān)系數(shù)計(jì)算、基因模塊的確定居暖、共表達(dá)網(wǎng)絡(luò)顽频、模塊與性狀關(guān)聯(lián)四個步驟。
第一步計(jì)算任意兩個基因之間的相關(guān)系數(shù)（Person Coefficient）太闺。為了衡量兩個基因是否具有相似表達(dá)模式糯景，一般需要設(shè)置閾值來篩選，高于閾值的則認(rèn)為是相似的省骂。但是這樣如果將閾值設(shè)為0.8蟀淮，那么很難說明0.8和0.79兩個是有顯著差別的。因此钞澳，WGCNA分析時采用相關(guān)系數(shù)加權(quán)值怠惶，即對基因相關(guān)系數(shù)取N次冪，使得網(wǎng)絡(luò)中的基因之間的連接服從無尺度網(wǎng)絡(luò)分布(scale-freenetworks)轧粟，這種算法更具生物學(xué)意義策治。
第二步通過基因之間的相關(guān)系數(shù)構(gòu)建分層聚類樹，聚類樹的不同分支代表不同的基因模塊逃延，不同顏色代表不同的模塊览妖≡簦基于基因的加權(quán)相關(guān)系數(shù)揽祥，將基因按照表達(dá)模式進(jìn)行分類，將模式相似的基因歸為一個模塊檩电。這樣就可以將幾萬個基因通過基因表達(dá)模式被分成了幾十個模塊拄丰，是一個提取歸納信息的過程府树。
直接上代碼,用了健明老師的數(shù)據(jù)先認(rèn)識一下流程及做什么
load('GSE48213-wgcna-input.RData')
if(T){
fpkm[1:4,1:4]
head(datTraits)
table(datTraits$subtype)
RNAseq_voom <- fpkm
## 因?yàn)閃GCNA針對的是基因進(jìn)行聚類，而一般我們的聚類是針對樣本用hclust即可料按，所以這個時候需要轉(zhuǎn)置
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)
#設(shè)置網(wǎng)絡(luò)構(gòu)建參數(shù)選擇范圍奄侠，計(jì)算無尺度分布拓?fù)渚仃? 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 自己按照流程期間有點(diǎn)問題，下午的時候跟另一位學(xué)徒再跑關(guān)于GEO數(shù)據(jù)挖掘载矿，把之前一點(diǎn)點(diǎn)小尾巴給解決掉垄潮。
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]]) #（提取表達(dá)矩陣）
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')
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
## 下面是畫PCA的必須操作弯洗，需要看說明書。
dat=t(dat)#畫PCA圖時要求是行名時樣本名逢勾，列名時探針名牡整，因此此時需要轉(zhuǎn)換
dat=as.data.frame(dat)#將matrix轉(zhuǎn)換為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)#現(xiàn)在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注釋轉(zhuǎn)換
rm(list = ls()) ## 魔幻操作逃贝，一鍵清空~
options(stringsAsFactors = F)
library(AnnoProbe)
load('GSE41804.Rdata')
# 選擇 pipe 獲取的是 冗余注釋，也就是說一個探針很有可能會對應(yīng)多個基因迫摔。
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這個矩陣按行操作句占，取每一行的中位數(shù)迫皱，將結(jié)果給到median這一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#對ids$symbol按照ids$median中位數(shù)從大到小排列的順序排序，將對應(yīng)的行賦值為一個新的ids
ids=ids[!duplicated(ids$symbol),]#將symbol這一列取取出重復(fù)項(xiàng)辖众，'!'為否卓起，即取出不重復(fù)的項(xiàng)峻堰，去除重復(fù)的gene 背亥，保留每個基因最大表達(dá)量結(jié)果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第一次出現(xiàn)的信息
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)
####構(gòu)建比較矩陣,設(shè)置用來對比的組別
contrast.matrix<-makeContrasts(contrasts=
c('HCC-Normal'),levels = design)
######limma三部曲,只需要?dú)w一化后的數(shù)據(jù)具钥、實(shí)驗(yàn)矩陣什燕、比較矩陣的輸入
fit <- lmFit(dat,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
####上面得到了p值等統(tǒng)計(jì)的結(jié)果载迄，topTable對p值校驗(yàn)岛啸，對基因排序
tT <- topTable(fit2, adjust="fdr", number=nrow(fit2))
dim(tT)
####提取出我們需要的指標(biāo)
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是指你的數(shù)據(jù)ID類型是屬于哪一類的
toType = c('ENTREZID'), #toType是指你要轉(zhuǎn)換成哪種ID類型徘钥，可以寫多種变骡，也可以只寫一種
OrgDb = org.Hs.eg.db)#Orgdb是指對應(yīng)的注釋包是哪個
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 因?yàn)橹白鲞^ 就直接一串代碼下來离赫，然后流程跑
image.png
image.png
image.png
image.png
