写在前面
最近又是忙碌的一米，做不完的手术，收不完的病人，前途堪忧，收入更是不堪入目。🥲
把之前的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>
点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
<center> <b>📍 往期精彩 <b> </center>
📍 <font size=1>🤩 WGCNA | 值得你深入学习的生信分析方法！~</font>
📍 <font size=1>🤩 ComplexHeatmap | 颜狗写的高颜值热图代码！</font>
📍 <font size=1>🤥 ComplexHeatmap | 你的热图注释还挤在一起看不清吗！？</font>
📍 <font size=1>🤨 Google | 谷歌翻译崩了我们怎么办！？（附完美解决方案）</font>
📍 <font size=1>🤩 scRNA-seq | 吐血整理的单细胞入门教程</font>
📍 <font size=1>🤣 NetworkD3 | 让我们一起画个动态的桑基图吧~</font>
📍 <font size=1>🤩 RColorBrewer | 再多的配色也能轻松搞定！~</font>
📍 <font size=1>🧐 rms | 批量完成你的线性回归</font>
📍 <font size=1>🤩 CMplot | 完美复刻Nature上的曼哈顿图</font>
📍 <font size=1>🤠 Network | 高颜值动态网络可视化工具</font>
📍 <font size=1>🤗 boxjitter | 完美复刻Nature上的高颜值统计图</font>
📍 <font size=1>🤫 linkET | 完美解决ggcor安装失败方案（附教程）</font>
📍 <font size=1>......</font>
本文由mdnice多平台发布
