生信文章复现学习 PMID：36898287(一)

平淡的牛马

已于 2023-11-09 22:20:56 修改

阅读量694

点赞数 3

分类专栏：文献复现文章标签：学习 r语言

于 2023-11-09 22:17:04 首次发布

本文链接：https://blog.csdn.net/weixin_40081460/article/details/134300289

版权

文献复现专栏收录该内容

3 篇文章 5 订阅

订阅专栏

文章学习，并尽量复现文章Identification of potential feature genes in non-alcoholic fatty liver disease using bioinformatics analysis and machine learning strategieshttps://pubmed.ncbi.nlm.nih.gov/36898287/

获取数据

本文直接通过GEO网站下载原始数据，根据原文数据结构下载多个GSE数据集，首先下载GSE89632_series_matrix以及平台注释文件GPL14951。本文没有用GEOquery包

数据预处理

rm(list = ls())##清空数据
library(tidyr)##加载各类R包
library(dplyr)
library(tibble)
library(data.table)
##一 数据预处理
##数据的处理转换贯穿整个分析流程
##（1）读取数据芯片数据和注释文件
setwd("D:\\BioR\\new\\NAFLD_Bio\\01preprocess") 
exprSet <- read.table("GSE89632_series_matrix.txt",
                      comment.char="!",
                      stringsAsFactors=F,
                      header=T)
annona <- data.table::fread("GPL14951.txt",data.table = F)
annona <- annona[,c(1,12)]##注释文件只需要探针名和gene symbol
colnames(annona) <- c('ID_REF','symbol')##方便合并，重命名
exprSet <- inner_join(annona,exprSet,by='ID_REF')##注释文件和表达谱合并
exprSet <- exprSet[,-1]
library(limma)
##本次采用平均值法处理重复基因
exprSet <- as.data.frame(avereps(exprSet[,-1],ID=exprSet$symbol))
##（2）读取分组信息，来自GEO
group <- fread("group.txt",data.table = F)
library(stringr)
##观察分组信息，发现‘_’为分隔，采用split函数分开字符串，并取第二个
group$group2=trimws(str_split(group$group,'_',simplify = T)[,2])
table(group$group2)
##此处提取HC和SS的分组信息，grep和grepl都可以，grep是提取包含某一或某些字符的行
##提取group2列里包含HC或者SS的行
groupf <-group[grep('HC|SS',group$group2),] 

ex_ma <- exprSet[,c(groupf$symbol)]##因为symbol变成了行名，所以不需要提取第一列

分组信息明确+gene symbol明确下一步可以做差异基因表达分析

差异基因分析（limma包）

##（2）数据正态化
library(limma) 
boxplot(ex_ma,outline=FALSE, notch=T, las=2)
### 该函数默认使用quntile 矫正差异 
ex_ma=normalizeBetweenArrays(ex_ma)
boxplot(ex_ma,outline=FALSE, notch=T, las=2)
## 这步把矩阵转换为数据框很重要
ex_ma <- as.data.frame(ex_ma)

###（4）差异分析
library(limma)
### 1.创建分组
### 这一步根据样本来就行，原则就是: 跟样本匹配，取决于样本的排序
group_df <- groupf$group2
### 分组变成向量，并且限定leves的顺序
### levels里面，把对照组放在前面
group_df <- factor(group_df,levels = c("HC","SS"))

### 主成分分析PCA：提前预测结果
### 行是样本列是基因
res.pca <- prcomp(t(ex_ma), scale = TRUE)
library(factoextra)
fviz_pca_ind(res.pca,col.ind = group_df)

### 构建比较矩阵
design <- model.matrix(~group_df)
### 比较矩阵命名
colnames(design) <- levels(group_df)

### 2.线性模型拟合
fit <- lmFit(ex_ma,design)
### 3.贝叶斯检验
fit2 <- eBayes(fit)
### 4.输出差异分析结果,其中coef的数目不能操过design的列数
### 此处的2代表的是design中第二列和第一列的比较
allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf) 
###5.选出符合条件的差异基因，以下两条代码意义相同
diffgene <- subset(allDiff,abs(logFC) >1.0 & adj.P.Val < 0.05)
#test <- allDiff[allDiff$adj.P.Val < 0.05 & abs(allDiff$logFC)>1,]  
save(allDiff,file='D:\\BioR\\new\\NAFLD_Bio\\01preprocess\\alldiff.Rdata')
save(ex_ma,file='D:\\BioR\\new\\NAFLD_Bio\\01preprocess\\ex_ma.Rdata')

存储差异基因及全部基因的结果，还有标准化后的表达数据，以备下一步分析。

本处limma包共得到409个DEGs和原文数量410很接近，继续加油，尝试WGCNA

WGCNA分析

WGCNA的目的是什么？

一言以蔽之，对划分不同基因，将相似表达的基因放在一起，探讨基因群和表型的关系

注意事项

样本数 >= 15
基因过滤方法 采用均值、方差、中位数、绝对中位差MAD等方法，过滤低表达或样本间变化小的基因，但不建议用差异分析的方法进行过滤
输入数据形式 如果有批次效应,需要先进行去除；处理RNAseq数据，需要采用DESeq2的varianceStabilizingTransformation方法，或将基因标准化后的数据（如FPKM、CPM等）进行log2(x+1)转化
经验软阈值power 当无向网络在power小于15或有向网络power小于30内，计算出的power无法达到要求时（即没有一个power值可以使无标度网络图谱结构R^2达到0.8且平均连接度降到100以下），可采用经验power值：
代码主要参考医学笔记，感激不尽。

代码开始

数据预处理

过滤基因数据

首先是过滤低质量的基因，官方建议先自己过滤一下，通过均值、方差、绝对中位差都可以，不建议通过差异分析筛选基因。因为原文纳入了所有基因（>20000个基因）构建WGCNA网络，所以此处我也纳入全部基因，同样，也可以采用官网方法来先过滤基因。

library(WGCNA)
#ex_ma是经过limma包标准化后的芯片数据
load('D:\\BioR\\new\\NAFLD_Bio\\02WGCNA\\ex_ma.Rdata')
datExpr0 <- as.data.frame(t(ex_ma))##转置数据
###########start handling the data
##1.看看有没有为空的值，需要剔除。这里下载的都是原作处理好的。
gsg <- goodSamplesGenes(datExpr0,verbose = 3)
gsg$allOK ##返回为true证明数据可以，如果false需要运行一段代码，网上可得

返回FALSE就用以下代码删除不符合条件的基因。

# 不OK就需要筛选
if (!gsg$allOK){
  # Optionally, print the gene and sample names that were removed:
  if (sum(!gsg$goodGenes)>0) 
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
  if (sum(!gsg$goodSamples)>0) 
    printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
  # Remove the offending genes and samples from the data:
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}

gsg <- goodSamplesGenes(datExpr0)
gsg$allOK

得到true就可以

过滤样本数据

datExpr1 <- datExpr0 #3防止出错，我们先建一个复制一个数据
##2.判断是否有离群样本
#（1）通过聚类分析，能看出是否有个别sample跟多数sample的距离较远，决定是否需要去除离群样本。
sampleTree = hclust(dist(datExpr1), method = "average")
par(cex = 0.7);
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
##本处有两个离群样本，需要剔除，根据Y轴提示选定hegiht为150可以剔除离群样本
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2) +
  #想用哪里切，就把“h = 150”和“cutHeight = 150”中的“150”换成你的cutoff
  abline(h = 150, col = "red") 
##（2）运行下段代码，剔除离群样本2个，得到新的表达量datExpr
clust = cutreeStatic(sampleTree, cutHeight = 150, minSize = 10)
keepSamples = (clust==1)
datExpr = datExpr1[keepSamples, ]
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
dim(datExpr)
##（3）再次画图验证是否已成功剔除
sampleTree = hclust(dist(datExpr), method = "average")
par(cex = 0.7);
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,cex.axis = 1.5, cex.main = 2)

这个地方第一次分类，发现有两个离群样本，因为远离其他样本，所以删掉，但是删掉后和原文作图不一致，不清楚原文作者是否删除，以及那两个样本是否属于离群样本。

datExpr就是准备好的转置后的表达矩阵了，接下来就用它进行后续的分析。

处理分组或表型信息

此处预先处理分组或表型等信息

#1.首先得到分组信息，分组信息直接拷贝的GEO信息，因为规律明显，容易使用R处理
group <- fread("group.txt",data.table = F)
library(stringr)
##观察分组信息，发现‘_’为分隔，采用split函数分开字符串，并取第二个
group$group2=trimws(str_split(group$group,'_',simplify = T)[,2])
table(group$group2)
##此处提取HC和SS的分组信息，grep和grepl都可以，grep是提取包含某一或某些字符的行
##提取group2列里包含HC或者SS的行
group <-group[grep('HC|SS',group$group2),] 
group <- group[,-2]
group <- subset(group,group$symbol !='GSM2385767') 
##此处删除两个离群样本，不知道为啥报错，只能用subset一个一个删
rownames(group) <- group[,1]##第一列生成列名，然后删除原第一列
samples <- dplyr::select(group,'group2')
samples$group2 <- ifelse(samples$group2=='SS',1,0)
samples$group3 <- ifelse(samples$group2=='0',1,0)
colnames(samples) <- c('SS','HC')

构建gene module

确定软阈值

软阈值一般要求R^2在0.9以上，最小也要在0.8以上。下面是挑选的代码，基本不用改。

#WGCNA分析的关键是找gene module。先选择合适的阈值，通过构建网络找gene module，
#找出来的gene module可信度如何？要做Preservation，去除not preserved module。
#这样找出的共表达的gene module就可以用于下一步分析了。
#（1）选择构建网络的合适阈值
#通过这步计算，找出scale free topology modle fit接近0.9的最小power（soft threshold），用于下一步构建网络。
powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
##做两个图看趋势，得出power=16（原文是16）
pdf("1Threshold.pdf",width = 10, height = 5)
par(mfrow = c(1,2))
cex1 = 0.9
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")+
  abline(h=0.90,col="red")

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()

结果就是看这两个图，选择达到0.9以上的软阈值，右边的图是连接度，选择逐渐平稳的阈值。

代码也会自动给出一个合适的软阈值：注意代码推荐的和你最终选择的不一定一致，根据原文，我们选择16

sft$powerEstimate

构建module

###（2）构建网络，找出gene module
net = blockwiseModules(datExpr, power = 16,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "MyTOM",
                       verbose = 3)##本代码运行完毕可找的划分的基因模块
table(net$colors)##具体看各个模块的基因数目
names(net)##可以查看net的项目

allowWGCNAThreads()#RStudio不支持enableWGCNAThreads()
##一步法构建网络和识别模块。计算时用到了BLAS，所以这一步是可以提速的，R自带的blas非常慢

module可视化

####（3）gene module的可视化
mergedColors = labels2colors(net$colors)

pdf("2module.pdf",width = 10, height = 5)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]], "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()

和原文差别较大，希望指正。

导出模块基因

##（4）把gene module输出到文件
#输出每个module内的基因，用于后续分析作图，随便你DIY
moduleColors = labels2colors(net$colors)

color<-unique(moduleColors)
for (i  in 1:length(color)) {
  y=t(assign(paste(color[i],"expr",sep = "."),datExpr[moduleColors==color[i]]))
  write.csv(y,paste(color[i],"csv",sep = "."),quote = F)
}

导出每一个模块的基因及其表达量，注意此处是删除了离群样本的表达量

Module深入分析

此处深入分析主要是指和表型的联系

模块与表型相关性热图的构建

###通过计算表型与module的相关性，找出相关性高的gene module，推测可能是因为它们造成了表型差异。
moduleLabelsAutomatic = net$colors
moduleColorsAutomatic = labels2colors(moduleLabelsAutomatic)
moduleColorsWW = moduleColorsAutomatic
MEs0 = moduleEigengenes(datExpr, moduleColorsWW)$eigengenes##表达矩阵不要错
MEsWW = orderMEs(MEs0)
modTraitCor = cor(MEsWW, samples, use = "p")##分组samples不要错
colnames(MEsWW)

modlues=MEsWW

modTraitP = corPvalueStudent(modTraitCor, nSamples)
textMatrix = paste(signif(modTraitCor, 2), "\n(", signif(modTraitP, 1), ")", sep = "")
dim(textMatrix) = dim(modTraitCor)
##做相关性热图
pdf("4Module-trait.pdf",width = 6, height = 6)
labeledHeatmap(Matrix = modTraitCor, xLabels = colnames(samples), yLabels = names(MEsWW), cex.lab = 0.5,  yColorWidth=0.01, 
               xColorWidth = 0.03,
               ySymbols = colnames(modlues), colorLabels = FALSE, colors = blueWhiteRed(50), 
               textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.5, zlim = c(-1,1)
               , main = paste("Module-trait relationships"))
dev.off()

计算GS及MM

（这一部分原文没有）

基因和性状的相关性：gene significance GS
基因和模块的相关性：module membership MM

选择我们感兴趣的性状（临床信息），这里是SS，也就是是否有NASH。

MEs = orderMEs(MEs0)
# Define variable weight containing the weight column of datTrait
cluster = as.data.frame(samples$SS) # 我们选SS
names(cluster) = "cluster"

# names (colors) of the modules
modNames = substring(names(MEs), 3)

# 基因和模块的相关性及P值
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));

names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");

# 基因和性状的相关性，这里是和样本亚型的相关性
geneTraitSignificance = as.data.frame(cor(datExpr, cluster, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));

names(geneTraitSignificance) = paste("GS.", names(cluster), sep="");
names(GSPvalue) = paste("p.GS.", names(cluster), sep="");

选择GS和MM都很高的分子

#画个图展示下每个基因的GS和MM的分布情况：
module = "yellow"
column = match(module, modNames);
moduleGenes = moduleColors==module;

sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for body weight",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
#接下来就是提取出GS和MM都很高的基因：
tmp1 <- rownames(geneModuleMembership)[abs(geneModuleMembership[moduleGenes, column])>0.7]
tmp2 <- rownames(geneTraitSignificance)[abs(geneTraitSignificance[moduleGenes, 1])>0.7] 
##此处0.7是可以升高或降低的

genes <- unique(intersect(tmp1,tmp2))
length(genes)

按照原文是采用相关系数最大的模块，最终采用yellow模块，共1500左右基因，和原文有出入。