TCGA的乳腺癌RNAseq数据WGCNA分析示例Word文档格式.docx
- 文档编号:17965386
- 上传时间:2022-12-12
- 格式:DOCX
- 页数:5
- 大小:18.45KB
TCGA的乳腺癌RNAseq数据WGCNA分析示例Word文档格式.docx
《TCGA的乳腺癌RNAseq数据WGCNA分析示例Word文档格式.docx》由会员分享,可在线阅读,更多相关《TCGA的乳腺癌RNAseq数据WGCNA分析示例Word文档格式.docx(5页珍藏版)》请在冰豆网上搜索。
,sep='
\t'
row.names=1)
dim(samples)
#[1]1003
expro=read.csv('
Merge_matrix.txt.cv.txt'
dim(expro)
#[1]24991100数据读取完成,从上述结果可以看出100个样本,有24991个基因,这么多基因全部用来做WGCNA很显然没有必要,我们只要选择一些具有代表性的基因就够了,这里我们采取的方式是选择在100个样本中方差较大的那些基因(意味着在不同样本中变化较大)继续命令:
m.vars=apply(expro,1,var)
expro.upper=expro[which(m.vars>
quantile(m.vars,probs=seq(0,1,0.25))[4]),]##选择方差最大的前25%个基因作为后续WGCNA的输入数据集通过上述步骤拿到了6248个基因的表达谱作为WGCNA的输入数据集,进一步的我们需要看看样本之间的差异情况datExpr=as.data.frame(t(expro.upper));
gsg=goodSamplesGenes(datExpr,verbose=3);
gsg$allOK
sampleTree=hclust(dist(datExpr),method='
average'
plot(sampleTree,main='
Sampleclusteringtodetectoutliers'
sub='
'
xlab='
)从图中可看出大部分样本表现比较相近,而有两个离群样本,对后续的分析可能造成影响,我们需要将其去掉,共得到98个样本clust=cutreeStatic(sampleTree,cutHeight=80000,minSize=10)
table(clust)
#clust
#01
#298
keepSamples=(clust==1)
datExpr=datExpr[keepSamples,]
nGenes=ncol(datExpr)
nSamples=nrow(datExpr)
save(datExpr,file='
FPKM-01-dataInput.RData'
)得到最终的数据矩阵之后,我们需要确定软阈值,从代码中可以看出pickSoftThreshold很简单,就两个参数,其他默认即可powers=c(c(1:
10),seq(from=12,to=20,by=2))
sft=pickSoftThreshold(datExpr,powerVector=powers,verbose=5)
##画图##
par(mfrow=c(1,2));
cex1=0.9;
plot(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab='
SoftThreshold(power)'
ylab='
ScaleFreeTopologyModelFit,signedR^2'
type='
n'
main=paste('
Scaleindependence'
));
text(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col='
red'
);
abline(h=0.90,col='
plot(sft$fitIndices[,1],sft$fitIndices[,5],
MeanConnectivity'
type='
Meanconnectivity'
))
text(sft$fitIndices[,1],sft$fitIndices[,5],labels=powers,cex=cex1,col='
)从图中可以看出这个软阈值选择7比较合适,选择软阈值7进行共表达模块挖掘pow=7
net=blockwiseModules(datExpr,power=pow,maxBlockSize=7000,
TOMType='
unsigned'
minModuleSize=30,
reassignThreshold=0,mergeCutHeight=0.25,
numericLabels=TRUE,pamRespectsDendro=FALSE,
saveTOMs=TRUE,
saveTOMFileBase='
FPKM-TOM'
verbose=3)
table(net$colors)
#openagraphicswindow
#sizeGrWindow(12,9)
#Convertlabelstocolorsforplotting
mergedColors=labels2colors(net$colors)
#Plotthedendrogramandthemodulecolorsunderneath
plotDendroAndColors(net$dendrograms[[1]],mergedColors[net$blockGenes[[1]]],
groupLabels=c('
Modulecolors'
'
GS.weight'
),
dendroLabels=FALSE,hang=0.03,
addGuide=TRUE,guideHang=0.05)从图中可以看出大部分基因在灰色区域,灰色部分一般认为是没有模块接受的,从这里也可以看出其实咱们选择的这些基因并不是特别好那么做到这一步了基本上共表达模块做完了,每个颜色代表一个共表达模块,统计看看各个模块下的基因个数:
那么得到模块之后下一步该做啥呢,或许很多人到这就不知道如何继续分析了这里就需要咱们利用这些模块搞事情了,举个例子如果你是整合的数据(整合lnc与gene),那么同时在某个模块中的基因和lncRNA咱们可以认为是共表达的,这便是lnc-gene共表达关系的获得途径之一了,进一步你可以根据该模块的基因-lnc-基因之间的关系绘制出共表达网络今天咱们这里不讲这个,而是跟表型关联,咱们已经拿到了这98个样本的ER、PR、HER2阳性阴性信息,那么进一步的咱们可以看看哪些共表达模块跟ER、PR、HER2阴性最相关,代码如下:
moduleLabelsAutomatic=net$colors
moduleColorsAutomatic=labels2colors(moduleLabelsAutomatic)
moduleColorsFemale=moduleColorsAutomatic
MEs0=moduleEigengenes(datExpr,moduleColorsFemale)$eigengenes
MEsFemale=orderMEs(MEs0)
samples=samples[match(row.names(datExpr),paste0(gsub('
-'
'
.'
row.names(samples)),'
.01'
)),]#匹配98个样本数据
trainDt=as.matrix(cbind(ifelse(samples[,1]=='
Positive'
0,1),#将阴性的样本标记为1
ifelse(samples[,2]=='
ifelse(samples[,3]=='
ifelse(samples[,1]=='
Negative'
&
samples[,2]=='
samples[,3]=='
1,0))#将三阴性的样本标记为1
)#得到一个表型的0-1矩阵
modTraitCor=cor(MEsFemale,trainDt,use='
p'
colnames(MEsFemale)
modTraitP=corPvalueStudent(modTraitCor,nSamples)
textMatrix=paste(signif(modTraitCor,2),'
\n('
signif(modTraitP,1),'
)'
sep='
dim(textMatrix)=dim(modTraitCor)
labeledHeatmap(Matrix=modTraitCor,xLabels=colnames(trainDt),yLabels=names(MEsFemale),
ySymbols=colnames(modlues),colorLabels=FALSE,colors=greenWhiteRed(50),
textMatrix=textMatrix,setStdMargins=FALSE,cex.text=0.5,zlim=c(-1,1)
main=paste('
Module-traitrelationships'
))最终找到几个共表达网络与三阴性表型最相关的模块。
到了这一步其实还没完,咱们已经拿到了最相关的模块比如上图中的yellow模块,那么yellow模块中的基因,哪些基因最重要,我们该如何获取首先我们需要计算每个基因与模块的相关关系modTraitCor=cor(MEsFemale,datExpr,use='
corYellow=modTraitCor[which(row.names(modTraitCor)=='
MEyellow'
),]
head(corYellow[order(-corYellow)])
#RAD51AP1HDAC2FOXM1NCAPD2TPI1NOP2
#0.92493540.90808720.89917330.88726070.87170500.8708449从上诉结果中可以看出RAD51AP1,HDAC2两个基因与yellow相关性最好,也就是说这两个基因是yellow模块的hub-gene进一步的咱们导出yellow模块的基因共表达关系使用cytoscape进行可视化,代码如下:
TOM=TOMsimilarityFromExpr(datExpr,power=pow);
probes=names(datExpr)
mc='
yellow'
mcInds=which(match(moduleColorsAutomatic,gsub('
^ME'
mc))==1)
modProbes=probes[mcInds]
modTOM=TOM[mcInds,mcInds];
dimnames(modTOM)=list(modProbes,modProbes)
cyt=exportNetworkToCytoscape(modTOM,
edgeFile=paste('
edges-'
mc,'
.txt'
sep='
nodeFile=paste('
nodes-'
weighted=TRUE,
threshold=median(modTOM),
nodeNames=modProbes,
#altNodeNames=modGenes,
nodeAttr=moduleColorsAutomatic[mcInds]);
到此基本结束了,敲这么多字。
。
好累。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- TCGA 乳腺癌 RNAseq 数据 WGCNA 分析 示例