GEO生信数据挖掘（六）实践案例——四分类结核病基因数据预处理分析

前面五节，我们使用阿尔兹海默症数据做了一个数据预处理案例，包括如下内容：

GEO生信数据挖掘（一）数据集下载和初步观察

GEO生信数据挖掘（二）下载基因芯片平台文件及注释

GEO生信数据挖掘（三）芯片探针ID与基因名映射处理

GEO生信数据挖掘（四）数据清洗（离群值处理、低表达基因、归一化、log2处理）

GEO生信数据挖掘（五）提取临床信息构建分组，分组数据可视化（绘制层次聚类图，绘制PCA图）

本节目录

结核病基因表达数据（GSE107994）观察

临床形状数据预处理

基因表达数据预处理

绘图观察数据

结核病基因表达数据（GSE107994）观察

由于，在数据分析过程，你拿的数据样式可能会有不同，本节我们以结核病基因表达数据（GSE107994）为例，做一个实践案例。该数据集的临床形状数据和基因表达数据是单独分开的，读取，和处理都需自己改动代码。

先来看看基因表达数据，这个探针注释工作已经完成了，不需要处理。

再看看临床形状数据，需要手工删除前面的注释，把后半部分规整的数据保留下来。

临床形状数据预处理

# 手工删除前面的注释，读文件，转置
pdata <- t(read.delim("GSE107994_series_matrix_clean.txt", header = TRUE, sep = "\t"))

# 手工删除前面的注释，读文件，转置
pdata <- t(read.delim("GSE107994_series_matrix_clean.txt", header = TRUE, sep = "\t"))
pdata  <-pdata[-1,]
pdata_info = pdata[,c(1,7)]
colnames(pdata_info) = c('geo_accession','type')#观察样本类型的取值都有哪些（结核，潜隐进展，对照和潜隐）
unique(pdata_info[,2])  
#"Leicester_Active_TB"  "Longitudnal_Leicester_LTBI_Progressor" "Leicester_Control"    #"Leicester_LTBI" group_data = as.data.frame(pdata_info)

处理前

处理后

增加不同的类型标签，根据需要，选取实验组和对照组


# 使用grepl函数判断字符串是否包含'Control'，并进行相应的修改
group_data$group_easy <- ifelse(grepl("Control", group_data$type ), "Control", "TB")# 使用grepl函数判断字符串是否包含特定内容，然后进行相应的修改
group_data$group_easy <- ifelse(grepl("Control", group_data$type), "Control",ifelse(grepl("LTBI", group_data$type), "LTBI","TB"))
# 使用grepl函数判断字符串是否包含特定内容，然后进行相应的修改
group_data$group_more <- ifelse(grepl("Control", group_data$type), "Control",ifelse(grepl("LTBI_Progressor", group_data$type), "LTBI_Progressor",ifelse(grepl("LTBI", group_data$type), "LTBI","TB")))#尝试把进展组排除出去save(group_data,file = "group_data.Rdata")

例如 我们可以进行 TB（结核） 和LTBI(潜隐结核)实验对照分析。

基因表达数据预处理


读取数据集

install.packages("openxlsx")
library(openxlsx)# 读基因表达矩阵，第一列为基因名ID
gse_info<- as.data.frame(read.xlsx("GSE107994_Raw.xlsx", sheet = 1))
colnames(gse_info)

后续运行代码过程中，发现基因名称中有全数字的情况，这里做删除操作。

library(dplyr)
dim(gse_info)
#基因里面有数字
gse_info <- gse_info[!grepl("^\\d+$", gse_info$ID), ]  #有效#基因名全为空
gse_info = gse_info[gse_info$ID != "",]  #无剔除
dim(gse_info) #[1] 58023   176#负值处理
gse_info[gse_info <= 0] <- 0.0001#重复值检查
table(duplicated(gse_info$ID))

分组数据条件筛选，TB（结核） 和LTBI(潜隐结核)

#+=====================================================
#================================================
#+========type分组数据条件筛选step3===========
#+====================================#预处理之前，先筛选出TB组和LTBI组 的数据
unique(group_data[,"group_more"])  #"TB"  "LTBI_Progressor" "Control"    "LTBI" #"TB"  "LTBI" 对照，则剔除 "LTBI_Progressor" "Control" 
geo_accession_TB_LTBI <- group_data[group_data$group_more == "LTBI_Progressor" | group_data$group_more == "Control","geo_accession"]
gse_TB_FTBI = gse_info[,!(names(gse_info) %in% geo_accession_TB_LTBI)]gse_TB_FTBI

低表达过滤(平均值小于1)


#+=====================================================
#================================================
#+========删除 低表达(平均值小于1)基因 step4===========
#+====================================
#+==============================#新增一列计算平均
gene_avg_expression <- rowMeans(gse_TB_FTBI[, -1])  # 计算每个基因的平均表达量，排除第一列（基因名）
#仅去除在所有样本里表达量都为零的基因（平均值小于1）
gse_TB_FTBI_filtered_genes_1 <- gse_TB_FTBI[gene_avg_expression >= 1, ]

低表达过滤方案二（保留样本表达的排名前50%的基因）

#+=================================================================
#============================================================
#+========删除 低表达(排名前50%)基因 step5===========
#+==========================================
#+================================#仅保留在一半以上样本里表达的基因# 计算基因表达矩阵每个基因的平均值
gene_means <- rowMeans(gse_TB_FTBI_filtered_genes_1[,-1])# 计算基因平均值的排序百分位数
gene_percentiles <- rank(gene_means) / length(gene_means)# 获取阈值
threshold <- 0.25  # 删除后25%的阈值
#threshold <- 0.5  # 删除后50%的阈值
# 根据阈值筛选低表达基因
gse_TB_FTBI_filtered_genes_2 <- gse_TB_FTBI_filtered_genes_1[gene_percentiles > threshold, ]# 打印筛选后的基因表达矩阵
dim(gse_TB_FTBI_filtered_genes_2) #[1] 17049   176

删除重复基因，取平均

#+=================================================================
#============================================================
#+========重复基因，取平均值 step6===========
#+==========================================
#+================================dim(filtered_genes_2) 
table(duplicated(filtered_genes_2$ID))#把重复的Symbol取平均值
averaged_data <- aggregate(.~ID , filtered_genes_2, mean, na.action = na.pass)  ##把重复的Symbol取平均值#把行名命名为SYMBOL
row.names(averaged_data) <- averaged_data$ID  
dim(averaged_data)#去掉缺失值
matrix_na = na.omit(averaged_data)  #删除Symbol列（一般是第一列）
matrix_final <- subset(matrix_na, select = -1) 
dim(matrix_final) #[1] 22687   175

离群值处理


#+=================================================================
#============================================================
#+========离群值处理 step7==========================
#+==========================================
#+================================#数据离群处理
#处理极端值
#定义向量极端值处理函数
#用于处理异常值，将超出一定范围的值替换为中位数，以减少异常值对后续分析的影响。
dljdz=function(x) {DOWNB=quantile(x,0.25)-1.5*(quantile(x,0.75)-quantile(x,0.25))UPB=quantile(x,0.75)+1.5*(quantile(x,0.75)-quantile(x,0.25))x[which(x<DOWNB)]=quantile(x,0.5)x[which(x>UPB)]=quantile(x,0.5)return(x)
}#第一列设置为行名
matrix_leave=matrix_final_TB_LTBIboxplot(matrix_leave,outline=FALSE, notch=T, las=2)  ##出箱线图
dim(matrix_leave)#处理离群值
matrix_leave_res=apply(matrix_leave,2,dljdz)boxplot(matrix_leave_res,outline=FALSE, notch=T, las=2)  ##出箱线图
dim(matrix_leave_res)

log2 处理

#+=================================================================
#============================================================
#+========log2 处理 step8==========================
#+==========================================
#+================================# limma的函数归一化，矫正差异  ，表达矩阵自动log2化#1.归一化不是绝对必要的，但是推荐进行归一化。
#有重复的样本中，应该不具备生物学意义的外部因素会影响单个样品的表达，
#例如中第一批制备的样品会总体上表达高于第二批制备的样品，假设所有样品表达值的范围和分布都应当相似，
#需要进行归一化来确保整个实验中每个样本的表达分布都相似。
#2.归一化要在log2标准化之前做library(limma) exprSet=normalizeBetweenArrays(matrix_leave_res)boxplot(exprSet,outline=FALSE, notch=T, las=2)  ##出箱线图## 这步把矩阵转换为数据框很重要
class(exprSet)   ##注释：此时数据的格式是矩阵（Matrix）
exprSet <- as.data.frame(exprSet)#标准化 表达矩阵自动log2化
qx <- as.numeric(quantile(exprSet, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||(qx[6]-qx[1] > 50 && qx[2] > 0) ||(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)#负值全部置为空
#exprSet[exprSet <= 0] <- 0.0001
#去掉缺失值
#exprSet = na.omit(exprSet)  #15654
#save (exprSet,file = "waitlog_data_TB_LTBI.Rdata")## 开始判断
if (LogC) { exprSet [which(exprSet  <= 0)] <- NaN## 取log2exprSet_clean <- log2(exprSet+1)  #@@@@是否加一 加一的话不产生负值@#@￥@#@#@%@%￥@@@@@@print("log2 transform finished")
}else{print("log2 transform not needed")
}boxplot(exprSet_clean,outline=FALSE, notch=T, las=2)  ##出箱线图dataset_TB_LTBI =exprSet_clean

绘图观察数据


#+=================================================================
#============================================================
#+========对照组不同颜色画箱线图 step9==========================
#+==========================================
#+================================# 使用grepl函数判断字符串是否包含'LTBI'，并进行颜色标记，为了画图
group_data_TB_LTBI$group_color <- ifelse(grepl("LTBI", group_data_TB_LTBI$group_more), "yellow", "blue")#画箱线图查看数据分布
group_list_color = group_data_TB_LTBI$group_color 
boxplot( data.frame(dataset_TB_LTBI),outline=FALSE,notch=T,col=group_list_color,las=2)dev.off()#+=================================================================
#============================================================
#+========绘制层次聚类图 step10==========================
#+==========================================
#+================================
#+#检查表达矩阵的样本名称，和分租信息的样本名称顺序,是否一致对应
colnames(dataset_TB_LTBI)
group_data_TB_LTBI$geo_accessionexprSet =dataset_TB_LTBI
#修改GSM的名字，改为分组信息
#colnames(exprSet)=paste(group_list,1:ncol(exprSet),sep = '')#定义nodePar
nodePar=list(lab.cex=0.6,pch=c(NA,19),cex=0.7,col='blue')
#聚类
hc=hclust(dist(t(exprSet))) #t()的意思是转置#绘图
plot(as.dendrogram(hc),nodePar = nodePar,horiz = TRUE)dev.off()#+=================================================================
#============================================================
#+========绘制PCA散点样本可视化图 step11===================
#+==========================================
#+================================##PCA图
#install.packages('ggfortify')
library(ggfortify)
df=as.data.frame(t(exprSet))  #转置后就变成了矩阵
dim(df)  #查看数据维度
dim(exprSet)df$group=group_data_TB_LTBI$group_more  #加入样本分组信息
autoplot(prcomp(df[,1:ncol(df)-1]),data=df,colour='group')  #PCA散点图dev.off()

至此，我们对两个数据集进行了预处理工作，下面我们可以对处理完毕的数据进行差异分析了。

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.rhkb.cn/news/154276.html

如若内容造成侵权/违法违规/事实不符，请联系长河编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

matlab相机标定实验

matlab相机标定实验

实验原理 1. 相机标定坐标系 相机的参数对目标的识别、定位精度有很大的影响&#xff0c;相机标定就是为了求出相机的内外参数。标定中有3个不同层次的坐标系&#xff1a;世界坐标系、相机坐标系和图像坐标系&#xff08;图像物理坐标系和图像像素坐标系&#xff09;。世界坐…
阅读更多...
清华智谱AI大模型ChatGLM-Pro申请开通教程

清华智谱AI大模型ChatGLM-Pro申请开通教程

清华智谱AI大模型ChatGLM-Pro申请开通教程 ChatGLM系列模型&#xff0c;包括ChatGLM-130B和ChatGLM-6B模型&#xff0c;支持相对复杂的自然语言指令&#xff0c;并且能够解决困难的推理类问题。其中&#xff0c;ChatGLM-6B模型吸引了全球超过 160 万人下载安装&#xff0c;该模…
阅读更多...
flutter 绘制原理探究

flutter 绘制原理探究

文章目录 Widget1、简介2、源码分析Element1、简介2、源码分析RenderObjectWidget 渲染过程总结思考Flutter 的核心设计思想便是“一切皆 Widget”，Widget 是 Flutter 功能的抽象描述，是视图的配置信息，同样也是数据的映射，是 Flutter 开发框架中最基本的概念。 在 Flutter…
阅读更多...
蓝桥杯基础---切面条

蓝桥杯基础---切面条

切面条 一根高筋拉面&#xff0c;中间切一刀&#xff0c;可以得到2根面条。 如果先对折1次&#xff0c;中间切一刀&#xff0c;可以得到3根面条。 如果连续对折2次&#xff0c;中间切一刀&#xff0c;可以得到5根面条。 那么&#xff0c;连续对折10次&#xff0c;中间切一刀…
阅读更多...
python curl2pyreqs 生成接口脚本

python curl2pyreqs 生成接口脚本

下载 curl2pyreqs 库 pip install curl2pyreqs -i https://pypi.tuna.tsinghua.edu.cn/simple 打开调试模式&#xff0c;在Network这里获取 接口的cURL 打开cmd窗口&#xff0c;输入curl2pyreqs&#xff0c;会自动生成接口代码 curl2pyreqs 执行接口脚本&#xff0c;返回响应…
阅读更多...
实现Promise所有核心功能和方法

实现Promise所有核心功能和方法

一直以来对Promise只是会用简单的方法&#xff0c;例如then&#xff0c;catch等&#xff0c;对于其余各种方法也只是简单了解&#xff0c;这次想要通过实现Promise来加深对Promise的使用 话不多说&#xff0c;直接开始&#xff0c;简单粗暴一步步来 一&#xff1a;了解Promise …
阅读更多...
Electron.js入门-构建第一个聊天应用程序

Electron.js入门-构建第一个聊天应用程序

什么是electron 电子是一个开源框架&#xff0c;用于使用web技术构建跨平台桌面应用程序&#xff1b;即&#xff1a; HTML、CSS和JavaScript&#xff1b;被集成为节点模块&#xff0c;我们可以为我们的应用程序使用节点的所有功能&#xff1b;组件&#xff0c;如数据库、Api休…
阅读更多...
FaceFusion：探索无限创意，创造独一无二的面孔融合艺术！

FaceFusion：探索无限创意，创造独一无二的面孔融合艺术！

FaceFusion&#xff1a;探索无限创意&#xff0c;创造独一无二的面孔融合艺术&#xff01; 它使用先进的图像处理技术&#xff0c;允许用户将不同的面部特征融合在一起&#xff0c;创造有趣和令人印象深刻的效果。这个项目的潜在应用包括娱乐、虚拟化妆和艺术创作&#xff0c;…
阅读更多...
(十五)VBA常用基础知识：正则表达式的使用

(十五)VBA常用基础知识：正则表达式的使用

vba正则表达式的说明 项目说明Pattern在这里写正则表达式&#xff0c;例&#xff1a;[\d]{2,4}IgnoreCase大小写区分&#xff0c;默认false&#xff1a;区分&#xff1b;true&#xff1a;不区分Globaltrue&#xff1a;全体检索&#xff1b;false&#xff1a;最小匹配Test类似p…
阅读更多...
自动求导，计算图示意图及pytorch实现

自动求导，计算图示意图及pytorch实现

pytorch实现 x1 torch.tensor(3.0, requires_gradTrue) y1 torch.tensor(2.0, requires_gradTrue) a x1 ** 2 b 3 * a c b * y1 c.backward() print(x1.grad) print(y1.grad) print(x1.grad 6 * x1 * y1) print(y1.grad 3 * (x1 ** 2))输出为&#xff1a; tensor(36.) …
阅读更多...
Go 复合数据类型之结构体与自定义类型

Go 复合数据类型之结构体与自定义类型

Go 复合数据类型之结构体与自定义类型 文章目录 Go 复合数据类型之结构体与自定义类型一、类型别名和自定义类型1.1 类型定义&#xff08;Type Definition&#xff09;简单示例 1.2 类型别名简单示例 1.3 类型定义和类型别名的区别 二、结构体2.1 结构体介绍2.2 结构体的定义2.…
阅读更多...
【Linux服务端搭建及使用】

【Linux服务端搭建及使用】

连接服务器的软件&#xff1a;mobaxterm 设置root 账号 sudo apt-get install passwd #安装passwd 设置方法 sudo passwd #设置root密码 su root #切换到root账户设置共享文件夹 一、强制删除原有环境 1.删除python rpm -qa|grep pytho…
阅读更多...
初学者如何选择：前端开发还是后端开发？

初学者如何选择：前端开发还是后端开发？

#开发做前端好还是后端好【话题征文】# 作为一名有多年开发经验的过来人&#xff0c;我认为前端开发和后端开发都有其独特的魅力和挑战。下面我将就我的个人经历和观点来分享一些关于前端开发和后端开发的看法。 首先&#xff0c;让我们将编程世界的大城市比作前端开发和后端开…
阅读更多...
vcf 文件如何修改染色体修改样本名称提取样本

vcf 文件如何修改染色体修改样本名称提取样本

大家好&#xff0c;我是邓飞。 对于vcf文件和plink文件是经常用的文件&#xff0c;对于基因型数据的处理&#xff0c;一般分为&#xff1a; 数据质控数据提取染色体修改名称样本修改名称 今天介绍一下vcf文件的三个处理方法&#xff1a; 1&#xff0c;染色体修改2&#xff…
阅读更多...
自助建站系统,一建建站系统api版,自动建站

自助建站系统,一建建站系统api版,自动建站

安装推荐php7.2或7.2以下都行 可使用虚拟主机或者服务器进行搭建。 分站进入网站后台 域名/admin 初始账号123456qq.com密码123456 找到后台的网站设置 将主站域名及你在主站的通信secretId和通信secretKey填进去。 即可正常使用 通信secretId和通信secretKey在主站的【账号…
阅读更多...
ruoyi 若依 前端vue npm install 运行vue前端

ruoyi 若依 前端vue npm install 运行vue前端

首次导入&#xff0c;需要先执行 npm install #进入到前端模块目录下 cd ruoyi-ui # 安装 npm install 启动后端项目 运行前端项目&#xff1a;运行成功后&#xff0c;会浏览器自动加载到前端首页&#xff08;或者 浏览器访问打印的两个地址&#xff09; # 运行 npm run dev 部…
阅读更多...
docker 安装 neo4j

docker 安装 neo4j

1. 安装所需的软件包 yum install -y yum-utils device-mapper-persistent-data lvm2 2. 设置阿里云仓库(国内仓库稳定) yum-config-manager --add-repo http://mirrors.aliyun.com/docker-ce/linux/centos/docker-ce.repo3. 查看docker容器版本 yum list docker-ce --showdupl…
阅读更多...
LeetCode-496-下一个更大元素

LeetCode-496-下一个更大元素

题目描述&#xff1a; 题目链接&#xff1a;LeetCode-496-下一个更大元素 解题思路&#xff1a; 方法一&#xff1a;暴力 方法二&#xff1a;单调栈 方法一代码实现&#xff1a; class Solution {public int[] nextGreaterElement(int[] nums1, int[] nums2) {// 最笨的方法&am…
阅读更多...
java基础 异常

java基础 异常

异常概述&#xff1a; try{ } catch{ }&#xff1a; package daysreplace;import com.sun.jdi.IntegerValue;import java.text.ParseException; import java.text.SimpleDateFormat; import java.util.Arrays; import java.util.Calendar; import java.util.Date; import java.…
阅读更多...
利用excel表格进行分包和组包

利用excel表格进行分包和组包

实际使用中&#xff0c;我们可能希望修改某几个数据之后&#xff0c;最终的数据包能够自动发动数据&#xff0c;类似于在给结构体变量修改数据&#xff0c;自动生成完整的结构体&#xff1b; excel语法 1&#xff1a;拆分数据 LEFT(A4,2) – 取A4单元格左边的两个数据 RIGHT(A4…
阅读更多...
最新文章
推荐文章

Copyright @ 2022~2023