🤠 WGCNA | 不止一个组的WGCNA怎么分析嘞！？~（一）

2023-05-06 本文已影响0人生信漫卷

写在前面

最近又是忙碌的一米，做不完的手术，收不完的病人，前途堪忧，收入更是不堪入目。🥲

把之前的WGCNA教程再补一补吧，之前介绍的是雌性鼠的表型数据分析，只有一组，相对简单。😅

今天开始我们介绍一下更为复杂的2组甚至多组分析。😘

用到的包

rm(list = ls())
library(tidyverse)      
library(WGCNA)

示例数据

用到的数据是雌性鼠和雄性鼠的肝脏表达谱，每个数据集包含大约130+个样本。🤨

femData <-  read.csv("./LiverFemale3600.csv")
dim(femData)

maleData <-  read.csv("./LiverMale3600.csv")
dim(maleData)

数据整合

我们先把两个数据集整合在一起，方便后续操作。🥳

nSets <-  2
setLabels <-  c("Female liver", "Male liver")
shortLabels <-  c("Female", "Male")

注意一下前面8列都是一些注释信息，第9列开始是表达谱。🤒

multiExpr <-  vector(mode = "list", length = nSets)

multiExpr[[1]] <-  list(data = as.data.frame(t(femData[-c(1:8)])))
names(multiExpr[[1]]$data) <-  femData$substanceBXH
rownames(multiExpr[[1]]$data) <-  names(femData)[-c(1:8)]

multiExpr[[2]] <-  list(data = as.data.frame(t(maleData[-c(1:8)])))
names(multiExpr[[2]]$data) <-  maleData$substanceBXH
rownames(multiExpr[[2]]$data) <-  names(maleData)[-c(1:8)]

检查一下数据集的格式是否正确吧。🧐

exprSize <-  checkSets(multiExpr)
exprSize

基本数据清理和异常值去除

5.1 去除不好的基因和样本

我们首先找到缺失样本数量过多的基因和样本。🥲

gsg <-  goodSamplesGenesMS(multiExpr, verbose = 3);
gsg$allOK

以下代码用于去除有问题的基因和样本，自取吧。🤒

if (!gsg$allOK)
{
if (sum(!gsg$goodGenes) > 0)
printFlush(paste("Removing genes:", paste(names(multiExpr[[1]]$data)[!gsg$goodGenes],
collapse = ", ")))
for (set in 1:exprSize$nSets)
{
if (sum(!gsg$goodSamples[[set]]))
printFlush(paste("In set", setLabels[set], "removing samples",
paste(rownames(multiExpr[[set]]$data)[!gsg$goodSamples[[set]]], collapse = ", ")))
multiExpr[[set]]$data = multiExpr[[set]]$data[gsg$goodSamples[[set]], gsg$goodGenes];
}
exprSize = checkSets(multiExpr)
}

5.2 样本聚类

接着我们对每个数据集分别进行cluster处理。🥰

sampleTrees = list()
for (set in 1:nSets)
{
sampleTrees[[set]] = hclust(dist(multiExpr[[set]]$data), method = "average")
}

可视化一下吧🤩

par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
xlab="", sub="", cex = 0.7)

5.3 去除异常值

可以看出来雌性鼠数据集中有一个异常值，雄性鼠数据集中没有，我们先可视化一下。✂️

baseHeight <-  16

cutHeights <-  c(16, 16*exprSize$nSamples[2]/exprSize$nSamples[1])

par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
xlab="", sub="", cex = 0.7)
abline(h=cutHeights[set], col = "red")
}

OK。我们现在真正的从数据中去除掉它！~🐵

for (set in 1:nSets)
{
labels = cutreeStatic(sampleTrees[[set]], cutHeight = cutHeights[set])
keep = (labels==1)
multiExpr[[set]]$data = multiExpr[[set]]$data[keep, ]
}
collectGarbage()

确定一下数据集的大小。🫶

exprSize <-  checkSets(multiExpr)
exprSize

载入特征数据

traitData <-  read.csv("ClinicalTraits.csv")
dim(traitData)
names(traitData)

去除一下我们不需要的数据。😭

allTraits <-  traitData[, -c(31, 16)]
allTraits <-  allTraits[, c(2, 11:36) ]
dim(allTraits)
names(allTraits)

把traits成list的格式来方便后面的应用。😉

Traits <-  vector(mode="list", length = nSets)
for (set in 1:nSets)
{
setSamples = rownames(multiExpr[[set]]$data)
traitRows = match(setSamples, allTraits$Mice)
Traits[[set]] = list(data = allTraits[traitRows, -1])
rownames(Traits[[set]]$data) = allTraits[traitRows, 1]
}
collectGarbage()

把genes数和samples数存起来，后面会用到。😜

nGenes <-  exprSize$nGenes
nSamples <-  exprSize$nSamples

save一下

初步处理完了数据我们就保存一下吧，毕竟辛辛苦苦清理好了数据。🎩

save(multiExpr, Traits, nGenes, nSamples, setLabels, shortLabels, exprSize,
file = "./Consensus-dataInput.RData")

<center>最后祝大家早日不卷!~</center>

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

本文由mdnice多平台发布