一生信文章复现学习 PMID：36898287 _基因

文章学习，并尽量复现文章 ofgenes in non- fatty liverusingand
目录
获取数据
数据预处理
差异基因分析（limma包）
WGCNA分析
WGCNA的目的是什么？
注意事项
代码开始
数据预处理
过滤基因数据
过滤样本数据
处理分组或表型信息
构建gene
确定软阈值
构建
可视化
导出模块基因
深入分析
模块与表型相关性热图的构建
计算GS及MM
选择GS和MM都很高的分子
获取数据
本文直接通过GEO网站下载原始数据，根据原文数据结构下载多个GSE数据集，首先下载ix以及平台注释文件。本文没有用包
数据预处理

文章插图

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 symbolcolnames(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）读取分组信息，来自GEOgroup <- 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 明确下一步可以做差异基因表达分析
差异基因分析（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')

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

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