当前位置: 首页 > news >正文

策划公司网站建设国际新闻最新消息美国

策划公司网站建设,国际新闻最新消息美国,dz做分类网站,用wordpress建站难吗检索到目标数据集后,开始数据挖掘,本文以阿尔兹海默症数据集GSE1297为例 目录 离群值处理 删除 低表达基因 函数归一化,矫正差异 数据标准化—log2处理 完整代码 上节围绕着探针ID和基因名称做了一些清洗工作,还做了重复值检查…

检索到目标数据集后,开始数据挖掘,本文以阿尔兹海默症数据集GSE1297为例

目录

离群值处理

删除 低表达基因

函数归一化,矫正差异

数据标准化—log2处理

 完整代码


上节围绕着探针ID和基因名称做了一些清洗工作,还做了重复值检查,空值删除操作。

#查看重复值
table(duplicated(matrix$Gene.Symbol))#去掉缺失值
matrix_na = na.omit(matrix)  #基因名称为空删除
matrix_final = matrix_na[matrix_na$Gene.Symbol != "",]

离群值处理

用于处理异常值,将超出一定范围的值替换为中位数,以减少异常值对后续分析的影响。

#数据离群处理
#处理极端值
#定义向量极端值处理函数
#用于处理异常值,将超出一定范围的值替换为中位数,以减少异常值对后续分析的影响。
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_finalboxplot(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)

删除 低表达基因

方案1 :仅去除在所有样本里表达量都为零的基因(或者平均值小于1)

方案2:仅保留在一半(50%,75%...自己选择)以上样本里表达的基因

#删除 低表达基因#仅去除在所有样本里表达量都为零的基因(平均值小于1)# 计算基因表达矩阵的平均值
gene_avg <- apply(matrix_final, 1, mean)
# 根据阈值筛选低表达基因
filtered_genes_1 <- matrix_final[gene_avg >= 1, ] # 表达量平均值小于1的过滤
dim(filtered_genes_1)#+================================
#常用过滤标准2(推荐):
#仅保留在一半以上样本里表达的基因# 计算基因表达矩阵每个基因的平均值
gene_means <- rowMeans(matrix_final)# 计算基因平均值的排序百分位数
gene_percentiles <- rank(gene_means) / length(gene_means)# 获取阈值
threshold <- 0.25  # 删除后25%的阈值
#threshold <- 0.5  # 删除后50%的阈值
# 根据阈值筛选低表达基因
filtered_genes_2 <- matrix_final[gene_percentiles > threshold, ]# 打印筛选后的基因表达矩阵
dim(filtered_genes_2)boxplot(filtered_genes_2,outline=FALSE, notch=T, las=2)  ##出箱线图
#dim(filtered_genes)#+删除 低表达基因   得到filtered_genes_1 删除平均表达小于1的   得到filtered_genes_2 删除后25%的

函数归一化,矫正差异


#1.归一化不是绝对必要的,但是推荐进行归一化。
#有重复的样本中,应该不具备生物学意义的外部因素会影响单个样品的表达,
#例如中第一批制备的样品会总体上表达高于第二批制备的样品,假设所有样品表达值的范围和分布都应当相似,
#需要进行归一化来确保整个实验中每个样本的表达分布都相似。
#2.归一化要在log2标准化之前做
 

library(limma) exprSet=normalizeBetweenArrays(filtered_genes_2)boxplot(exprSet,outline=FALSE, notch=T, las=2)  ##出箱线图## 这步把矩阵转换为数据框很重要
class(exprSet)   ##注释:此时数据的格式是矩阵(Matrix)
exprSet <- as.data.frame(exprSet)

数据标准化—log2处理

如果表达量的数值在50以内,通常是经过log2转化后的。如果数字在几百几千,则是未经转化的。

GSE数据集的注释部分会有说明,常见的标准化处理方法有3种:RMA算法、GC-RMA算法、MAS5算法

#标准化 表达矩阵自动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)## 开始判断
if (LogC) { exprSet [which(exprSet  <= 0)] <- NaN## 取log2exprSet_clean <- log2(exprSet)print("log2 transform finished")
}else{print("log2 transform not needed")
}

筛选后的基因表达矩阵的箱线图

经过数据清洗后的箱线图

 完整代码


# 安装并加载GEOquery包
library(GEOquery)# 指定GEO数据集的ID
gse_id <- "GSE1297"# 使用getGEO函数获取数据集的基础信息
gse_info <- getGEO(gse_id, destdir = ".", AnnotGPL = F ,getGPL = F) # Failed to download ./GPL96.soft.gz!# 提取基因表达矩阵
expression_data <- exprs(gse_info[[1]])# 提取注释信息
#annotation <- featureData(gse_info[[1]])#查看平台文件列名
colnames(annotation)#打印项目文件列表
dir() # 读取芯片平台文件txt
platform_file <- read.delim("GPL96-57554.txt", header = TRUE, sep = "\t", comment.char = "#")#查看平台文件列名
colnames(platform_file)# 假设芯片平台文件中有两列,一列是探针ID,一列是基因名
#probe_names <- platform_file$ID
#gene_symbols <- platform_file$Gene.Symbol
platform_file_set=platform_file[,c(1,11)]#一个探针对应多个基因名,保留第一个基因名
ids = platform_file_set
library(tidyverse)
test_function <- apply(ids,1,function(x){paste(x[1],str_split(x[2],'///',simplify=T),sep = "...")})
x = tibble(unlist(test_function))colnames(x) <- "ttt" 
ids <- separate(x,ttt,c("ID","Gene.Symbol"),sep = "\\...")
dim(ids)#将Matrix格式表达矩阵转换为data.frame格式
exprSet <- data.frame(expression_data)
dim(exprSet)#给表达矩阵新增加一列ID
exprSet$ID <- rownames(exprSet) # 得到表达矩阵,行名为ID,需要转换,新增一列
dim(exprSet)
#矩阵表达文件和平台文件有相同列‘ID’,使用merge函数合并
express <- merge(x = exprSet, y = ids, by.x = "ID")#删除探针ID列
express$ID =NULLdim(express) matrix = express
dim(matrix)
#查看多少个基因重复了
table(duplicated(matrix$Gene.Symbol))#把重复的Symbol取平均值
matrix <- aggregate(.~Gene.Symbol, matrix, mean)  ##把重复的Symbol取平均值
row.names(matrix) <- matrix$Gene.Symbol  #把行名命名为SYMBOLdim(matrix)matrix_na = na.omit(matrix)   #去掉缺失值
dim(matrix_na)matrix_final = matrix_na[matrix_na$Gene.Symbol != "",]
dim(matrix_final)matrix_final <- subset(matrix_final, select = -1)  #删除Symbol列(一般是第一列)
dim(matrix_final)
#+  经过注释、探针名基因名处理、删除基因名为空值、删除缺失值 得到最终 matrix_final
#+==================================================================================
#+========================================================================================#+==================================================================================
#+========================================================================================
#+增加#(2)离群值处理#数据离群处理
#处理极端值
#定义向量极端值处理函数
#用于处理异常值,将超出一定范围的值替换为中位数,以减少异常值对后续分析的影响。
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_finalboxplot(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)#+增加#(2)离群值处理
#+==================================================================================
#+========================================================================================#+==================================================================================
#+========================================================================================
#删除 低表达基因#仅去除在所有样本里表达量都为零的基因(平均值小于1)# 计算基因表达矩阵的平均值
gene_avg <- apply(matrix_final, 1, mean)
# 根据阈值筛选低表达基因
filtered_genes_1 <- matrix_final[gene_avg >= 1, ] # 表达量平均值小于1的过滤
dim(filtered_genes_1)#+================================
#常用过滤标准2(推荐):
#仅保留在一半以上样本里表达的基因# 计算基因表达矩阵每个基因的平均值
gene_means <- rowMeans(matrix_final)# 计算基因平均值的排序百分位数
gene_percentiles <- rank(gene_means) / length(gene_means)# 获取阈值
threshold <- 0.25  # 删除后25%的阈值
#threshold <- 0.5  # 删除后50%的阈值
# 根据阈值筛选低表达基因
filtered_genes_2 <- matrix_final[gene_percentiles > threshold, ]# 打印筛选后的基因表达矩阵
dim(filtered_genes_2)boxplot(filtered_genes_2,outline=FALSE, notch=T, las=2)  ##出箱线图
#dim(filtered_genes)#+删除 低表达基因   得到filtered_genes_1 删除平均表达小于1的   得到filtered_genes_2 删除后25%的
#+==================================================================================
#+========================================================================================#+==================================================================================
#+========================================================================================
### limma的函数归一化,矫正差异  ,表达矩阵自动log2化#1.归一化不是绝对必要的,但是推荐进行归一化。
#有重复的样本中,应该不具备生物学意义的外部因素会影响单个样品的表达,
#例如中第一批制备的样品会总体上表达高于第二批制备的样品,假设所有样品表达值的范围和分布都应当相似,
#需要进行归一化来确保整个实验中每个样本的表达分布都相似。
#2.归一化要在log2标准化之前做library(limma) exprSet=normalizeBetweenArrays(filtered_genes_2)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)## 开始判断
if (LogC) { exprSet [which(exprSet  <= 0)] <- NaN## 取log2exprSet_clean <- log2(exprSet)print("log2 transform finished")
}else{print("log2 transform not needed")
}boxplot(exprSet_clean,outline=FALSE, notch=T, las=2)  ##出箱线图
### limma的函数归一化,矫正差异  ,表达矩阵自动log2化 得到exprSet_clean
#+==================================================================================
#+========================================================================================# saving all data to the path
save(exprSet_clean, file ="exprSet_clean_75percent_filter.RData")

上述处理得到了干净的基因表达矩阵,数据部分已经没有问题,但是在做数据挖掘(差异分析、富集分析等)之前还有一项准备工作,要将数据样本进行分组,即患病组、对照组。

http://www.hkea.cn/news/901213/

相关文章:

  • 网站设计详细设计有没有好用的网站推荐
  • 没有货源可以开网店吗网站更新seo
  • 淄博有做网站的吗百度搜索排名怎么收费
  • wordpress页面添加自定义字段木卢seo教程
  • 长寿网站制作保定seo排名外包
  • 域名和网站一样吗电商运营推广怎么做
  • css个人简介网站怎么做b2b网站免费推广平台
  • 网站建设中企动力上海百度广告投诉电话客服24小时
  • 深圳靠谱的电商公司正版搜索引擎优化
  • 自己如何做团购网站腾讯云建站
  • 怀化招标网站磁力狗bt
  • 佛山网站建设服务公司培训机构查询网
  • 海尔集团电商网站建设考证培训机构
  • 动漫制作专业的高职实训室福州整站优化
  • 织梦商城网站模板免费下载怎么在网上做推广
  • asp做网站用什么写脚本温岭网络推广
  • 怎么建设外贸网站免费发seo外链平台
  • 郴州是几线城市武汉网站seo推广公司
  • 网站开发工程师求职信焊工培训内容
  • 铜陵公司做网站中国网站排名100
  • 我要建一个网站泰州百度公司代理商
  • php响应式网站模板vi设计公司
  • 随身wifi网站设置广告投放是做什么的
  • 中企动力做网站的优势网络销售平台有哪些软件
  • 网站建设的费用如何查看百度搜索指数
  • 自己做网站需要什么seo的基本步骤
  • 视频直播app开发网站南京最新消息今天
  • 溧阳手机网站哪里做万网域名注册官网查询
  • 网站维护收费推广产品吸引人的句子
  • 怎么用一个主机做多个网站许昌网络推广公司