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

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

    目录

    离群值处理

    删除 低表达基因

    函数归一化，矫正差异

    数据标准化—log2处理

     完整代码

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

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

    离群值处理

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

    1. #数据离群处理
    2. #处理极端值
    3. #定义向量极端值处理函数
    4. #用于处理异常值，将超出一定范围的值替换为中位数，以减少异常值对后续分析的影响。
    5. dljdz=function(x) {
    6.   DOWNB=quantile(x,0.25)-1.5*(quantile(x,0.75)-quantile(x,0.25))
    7.   UPB=quantile(x,0.75)+1.5*(quantile(x,0.75)-quantile(x,0.25))
    8.   x[which(x<DOWNB)]=quantile(x,0.5)
    9.   x[which(x>UPB)]=quantile(x,0.5)
    10.   return(x)
    11. }
    12.  
    13. #第一列设置为行名
    14. matrix_leave=matrix_final
    15.  
    16. boxplot(matrix_leave,outline=FALSE, notch=T, las=2)  ##出箱线图
    17. dim(matrix_leave)
    18.  
    19. #处理离群值
    20. matrix_leave_res=apply(matrix_leave,2,dljdz)
    21.  
    22. boxplot(matrix_leave_res,outline=FALSE, notch=T, las=2)  ##出箱线图
    23. dim(matrix_leave_res)

    删除 低表达基因

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

    方案2：仅保留在一半（50%，75%...自己选择）以上样本里表达的基因

    1. #删除 低表达基因
    2.  
    3. #仅去除在所有样本里表达量都为零的基因（平均值小于1）
    4.  
    5. # 计算基因表达矩阵的平均值
    6. gene_avg <- apply(matrix_final, 1, mean)
    7. # 根据阈值筛选低表达基因
    8. filtered_genes_1 <- matrix_final[gene_avg >= 1, ] # 表达量平均值小于1的过滤
    9. dim(filtered_genes_1)
    10.  
    11.  
    12. #+================================
    13. #常用过滤标准2(推荐)：
    14. #仅保留在一半以上样本里表达的基因
    15.  
    16. # 计算基因表达矩阵每个基因的平均值
    17. gene_means <- rowMeans(matrix_final)
    18.  
    19. # 计算基因平均值的排序百分位数
    20. gene_percentiles <- rank(gene_means) / length(gene_means)
    21.  
    22. # 获取阈值
    23. threshold <- 0.25  # 删除后25%的阈值
    24. #threshold <- 0.5  # 删除后50%的阈值
    25. # 根据阈值筛选低表达基因
    26. filtered_genes_2 <- matrix_final[gene_percentiles > threshold, ]
    27.  
    28. # 打印筛选后的基因表达矩阵
    29. dim(filtered_genes_2)
    30.  
    31. boxplot(filtered_genes_2,outline=FALSE, notch=T, las=2)  ##出箱线图
    32. #dim(filtered_genes)
    33.  
    34. #+删除 低表达基因   得到filtered_genes_1 删除平均表达小于1的   得到filtered_genes_2 删除后25%的

    函数归一化，矫正差异


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

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

    数据标准化—log2处理

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

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

    1.  
    2.  
    3. #标准化 表达矩阵自动log2化
    4. qx <- as.numeric(quantile(exprSet, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
    5. LogC <- (qx[5] > 100) ||
    6.   (qx[6]-qx[1] > 50 && qx[2] > 0) ||
    7.   (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
    8.  
    9.  
    10. ## 开始判断
    11. if (LogC) { 
    12.   exprSet [which(exprSet  <= 0)] <- NaN
    13.   ## 取log2
    14.   exprSet_clean <- log2(exprSet)
    15.   print("log2 transform finished")
    16. }else{
    17.   print("log2 transform not needed")
    18. }

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

    经过数据清洗后的箱线图

     完整代码

    1.  
    2. # 安装并加载GEOquery包
    3. library(GEOquery)
    4.  
    5. # 指定GEO数据集的ID
    6. gse_id <- "GSE1297"
    7.  
    8. # 使用getGEO函数获取数据集的基础信息
    9. gse_info <- getGEO(gse_id, destdir = ".", AnnotGPL = F ,getGPL = F) # Failed to download ./GPL96.soft.gz!
    10.  
    11. # 提取基因表达矩阵
    12. expression_data <- exprs(gse_info[[1]])
    13.  
    14. # 提取注释信息
    15. #annotation <- featureData(gse_info[[1]])
    16.  
    17.  
    18. #查看平台文件列名
    19. colnames(annotation)
    20.  
    21. #打印项目文件列表
    22. dir() 
    23.  
    24. # 读取芯片平台文件txt
    25. platform_file <- read.delim("GPL96-57554.txt", header = TRUE, sep = "\t", comment.char = "#")
    26.  
    27. #查看平台文件列名
    28. colnames(platform_file)
    29.  
    30. # 假设芯片平台文件中有两列，一列是探针ID，一列是基因名
    31. #probe_names <- platform_file$ID
    32. #gene_symbols <- platform_file$Gene.Symbol
    33. platform_file_set=platform_file[,c(1,11)]
    34.  
    35.  
    36. #一个探针对应多个基因名，保留第一个基因名
    37. ids = platform_file_set
    38. library(tidyverse)
    39. test_function <- apply(ids,
    40.                        1,
    41.                        function(x){
    42.                          paste(x[1],
    43.                                str_split(x[2],'///',simplify=T),
    44.                                sep = "...")
    45.                        })
    46. x = tibble(unlist(test_function))
    47.  
    48. colnames(x) <- "ttt" 
    49. ids <- separate(x,ttt,c("ID","Gene.Symbol"),sep = "\\...")
    50. dim(ids)
    51.  
    52. #将Matrix格式表达矩阵转换为data.frame格式
    53. exprSet <- data.frame(expression_data)
    54. dim(exprSet)
    55.  
    56. #给表达矩阵新增加一列ID
    57. exprSet$ID <- rownames(exprSet) # 得到表达矩阵，行名为ID，需要转换，新增一列
    58. dim(exprSet)
    59. #矩阵表达文件和平台文件有相同列‘ID’，使用merge函数合并
    60. express <- merge(x = exprSet, y = ids, by.x = "ID")
    61.  
    62. #删除探针ID列
    63. express$ID =NULL
    64.  
    65.  
    66. dim(express) 
    67.  
    68. matrix = express
    69. dim(matrix)
    70. #查看多少个基因重复了
    71. table(duplicated(matrix$Gene.Symbol))
    72.  
    73.  
    74. #把重复的Symbol取平均值
    75. matrix <- aggregate(.~Gene.Symbol, matrix, mean)  ##把重复的Symbol取平均值
    76. row.names(matrix) <- matrix$Gene.Symbol  #把行名命名为SYMBOL
    77.  
    78. dim(matrix)
    79.  
    80. matrix_na = na.omit(matrix)   #去掉缺失值
    81. dim(matrix_na)
    82.  
    83. matrix_final = matrix_na[matrix_na$Gene.Symbol != "",]
    84. dim(matrix_final)
    85.  
    86. matrix_final <- subset(matrix_final, select = -1)  #删除Symbol列（一般是第一列）
    87. dim(matrix_final)
    88. #+  经过注释、探针名基因名处理、删除基因名为空值、删除缺失值 得到最终 matrix_final
    89. #+==================================================================================
    90. #+========================================================================================
    91.  
    92.  
    93.  
    94.  
    95. #+==================================================================================
    96. #+========================================================================================
    97. #+增加#(2)离群值处理
    98.  
    99. #数据离群处理
    100. #处理极端值
    101. #定义向量极端值处理函数
    102. #用于处理异常值，将超出一定范围的值替换为中位数，以减少异常值对后续分析的影响。
    103. dljdz=function(x) {
    104.   DOWNB=quantile(x,0.25)-1.5*(quantile(x,0.75)-quantile(x,0.25))
    105.   UPB=quantile(x,0.75)+1.5*(quantile(x,0.75)-quantile(x,0.25))
    106.   x[which(x<DOWNB)]=quantile(x,0.5)
    107.   x[which(x>UPB)]=quantile(x,0.5)
    108.   return(x)
    109. }
    110.  
    111. #第一列设置为行名
    112. matrix_leave=matrix_final
    113.  
    114. boxplot(matrix_leave,outline=FALSE, notch=T, las=2)  ##出箱线图
    115. dim(matrix_leave)
    116.  
    117. #处理离群值
    118. matrix_leave_res=apply(matrix_leave,2,dljdz)
    119.  
    120. boxplot(matrix_leave_res,outline=FALSE, notch=T, las=2)  ##出箱线图
    121. dim(matrix_leave_res)
    122.  
    123. #+增加#(2)离群值处理
    124. #+==================================================================================
    125. #+========================================================================================
    126.  
    127.  
    128.  
    129.  
    130.  
    131. #+==================================================================================
    132. #+========================================================================================
    133. #删除 低表达基因
    134.  
    135. #仅去除在所有样本里表达量都为零的基因（平均值小于1）
    136.  
    137. # 计算基因表达矩阵的平均值
    138. gene_avg <- apply(matrix_final, 1, mean)
    139. # 根据阈值筛选低表达基因
    140. filtered_genes_1 <- matrix_final[gene_avg >= 1, ] # 表达量平均值小于1的过滤
    141. dim(filtered_genes_1)
    142.  
    143.  
    144. #+================================
    145. #常用过滤标准2(推荐)：
    146. #仅保留在一半以上样本里表达的基因
    147.  
    148. # 计算基因表达矩阵每个基因的平均值
    149. gene_means <- rowMeans(matrix_final)
    150.  
    151. # 计算基因平均值的排序百分位数
    152. gene_percentiles <- rank(gene_means) / length(gene_means)
    153.  
    154. # 获取阈值
    155. threshold <- 0.25  # 删除后25%的阈值
    156. #threshold <- 0.5  # 删除后50%的阈值
    157. # 根据阈值筛选低表达基因
    158. filtered_genes_2 <- matrix_final[gene_percentiles > threshold, ]
    159.  
    160. # 打印筛选后的基因表达矩阵
    161. dim(filtered_genes_2)
    162.  
    163. boxplot(filtered_genes_2,outline=FALSE, notch=T, las=2)  ##出箱线图
    164. #dim(filtered_genes)
    165.  
    166. #+删除 低表达基因   得到filtered_genes_1 删除平均表达小于1的   得到filtered_genes_2 删除后25%的
    167. #+==================================================================================
    168. #+========================================================================================
    169.  
    170.  
    171.  
    172.  
    173.  
    174. #+==================================================================================
    175. #+========================================================================================
    176. ### limma的函数归一化，矫正差异  ，表达矩阵自动log2化
    177.  
    178. #1.归一化不是绝对必要的，但是推荐进行归一化。
    179. #有重复的样本中，应该不具备生物学意义的外部因素会影响单个样品的表达，
    180. #例如中第一批制备的样品会总体上表达高于第二批制备的样品，假设所有样品表达值的范围和分布都应当相似，
    181. #需要进行归一化来确保整个实验中每个样本的表达分布都相似。
    182. #2.归一化要在log2标准化之前做
    183.  
    184.  
    185. library(limma) 
    186.  
    187. exprSet=normalizeBetweenArrays(filtered_genes_2)
    188.  
    189. boxplot(exprSet,outline=FALSE, notch=T, las=2)  ##出箱线图
    190.  
    191. ## 这步把矩阵转换为数据框很重要
    192. class(exprSet)   ##注释：此时数据的格式是矩阵（Matrix）
    193. exprSet <- as.data.frame(exprSet)
    194.  
    195.  
    196. #标准化 表达矩阵自动log2化
    197. qx <- as.numeric(quantile(exprSet, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
    198. LogC <- (qx[5] > 100) ||
    199.   (qx[6]-qx[1] > 50 && qx[2] > 0) ||
    200.   (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
    201.  
    202.  
    203. ## 开始判断
    204. if (LogC) { 
    205.   exprSet [which(exprSet  <= 0)] <- NaN
    206.   ## 取log2
    207.   exprSet_clean <- log2(exprSet)
    208.   print("log2 transform finished")
    209. }else{
    210.   print("log2 transform not needed")
    211. }
    212.  
    213. boxplot(exprSet_clean,outline=FALSE, notch=T, las=2)  ##出箱线图
    214. ### limma的函数归一化，矫正差异  ，表达矩阵自动log2化 得到exprSet_clean
    215. #+==================================================================================
    216. #+========================================================================================
    217.  
    218.  
    219.  
    220. # saving all data to the path
    221. save(exprSet_clean, file ="exprSet_clean_75percent_filter.RData")

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

  • 相关阅读:
    一文看懂推荐系统：物品冷启02：简单的召回通道
    【python版CV】- 银行卡号识别项目
    Codeforcess834
    Chromium 开发指南2024 Mac篇-开始编译Chromium（五）
    虹科活动 | 探索全新AR应用时代，虹科AR VIP研讨会广州场回顾！
    测量工具----示波器
    关联容器(字典)map
    MYSQL CONCAT函数
    第3章 初识SqlSugarCore之ConfigureOptions注入实现
    DPDK代码目录结构
  • 原文地址：https://blog.csdn.net/zzh1464501547/article/details/133326080