往前一步是优秀,退后一步是懵懂

2022-06-27 20:28:32 浏览数 (1)

我们的生信入门班和数据挖掘线上直播课程已经有了三年多的历史,培养了一波又一波优秀的生信人才。前面提到R语言授课时的超纲练习题,已经分享过两位优秀学员的答案。

  • 系统学习+主动探索,是最舒适的入门学习方式!
  • 超纲练习题不超纲

下面继续来看优秀学员Dr.luka的分享:

R语言超纲练习题

(生信技能树优秀学员Dr.luka)

  • 数据挖掘(GEO,TCGA,单细胞)2022年6月场,快速了解一些生物信息学应用图表
  • 生信入门课-2022年6月场,你的生物信息学第一课

❝最好的输入就是输出。这篇笔记大量参考了同学们的分享,并结合自己的思路进行修改和补充。 ❞

❝十分感谢 佳男同学 的分享笔记系统学习+主动探索,是最舒适的入门学习方式! 徐谦同学 的分享笔记超纲练习题不超纲 ❞

1.读取文件
代码语言:javascript复制
exp <- read.csv('exp.csv') #第一次读入不要用row.names=1,防止报错
soft <- read.table("soft.txt",sep = "t",header = T)

exp内容展示

soft内容展示

2.探针过滤

由于实际情况可能存在注释文件的探针(probe_id)与表达矩阵探针不能完全对应的情况,因此在进行基因名转换之前,需要把探针进行过滤,留下有效的探针

代码语言:javascript复制
table(exp$X %in% soft$ID) #如果有FALSE则说明有不对应情况
dim(exp) #过滤前探针数
exp <- exp[exp$X %in% soft$ID,]
dim(exp) #过滤后探针数
3.删除重复的基因名,整理表达矩阵
方法1.直接删除重复基因,保留下标最小的
代码语言:javascript复制
#1.合并探针信息
colnames(exp)[1] <- 'ID'
exp_new <- merge(exp,soft,by = 'ID') #由于soft$GeneName和exp$ID内容相同,可以直接合并
#dplyr包inner_join()和base包merge()用法相近

#2.删除重复基因
exp_new <- exp_new[!duplicated(exp_new$GeneName),] #直接删除重复基因
rownames(exp_new) <- exp_new$GeneName
exp_new <- exp_new[,paste0('S',1:6)] %>% as.matrix() #提取表达矩阵

> head(exp_new)
                       S1        S2        S3        S4        S5        S6
LOC641522        6.536468  7.739101  4.407368 10.569904  9.306856  7.845223
RPL31           12.711015  5.611023  8.932942 11.769962 11.554846 12.940682
ENST00000292530  9.624848 10.393935 11.535562  8.479768  6.753438  6.582785
HARS2            5.921307  8.075039  8.010578  8.548929  5.495655  7.350851
KLK3             8.710169  6.204116 11.232571 10.322088  5.984473 12.062789
ADD1            13.510201 10.956088  7.861020  9.365741  4.412002  7.038853
方法2.重复基因的表达量取平均值
代码语言:javascript复制
> library(tibble)
> library(dplyr)

# 1.合并探针信息并整理
> colnames(exp)[1] <- 'ID'
> exp_new <- merge(exp,soft,by = 'ID')
> exp_new <- exp_new[,c('GeneName',paste0('S',1:6))]

# 2.重复基因表达量取平均值
> dim(exp_new)
[1] 1000    7
> exp_mean <- aggregate(x = exp_new[,colnames(exp_new)!= 'GeneName'], # x是要进行分析的数据
                    by = list(exp_new$GeneName), # by是进行运算的分组(list形式出现)
                    FUN = mean) # FUN是运算函数
> dim(exp_mean)
[1] 946   7
> rownames(exp_mean) <- exp_mean$Group.1
> exp_mean$Group.1 <- NULL
> exp_mean[1:5,]
                S1       S2        S3       S4        S5        S6
15E1.2    5.826481 6.609382  6.807679 7.814539  8.950446  5.897408
AB016902  9.543493 9.385374  5.995681 8.729387 12.121148  9.337484
ABCA1     8.773911 8.026560  8.090828 7.892680  6.316645  6.701494
ABCC6     9.228167 9.337685  8.517037 6.703833 10.921183 10.143636
ACBD3    10.042010 8.634697 14.158641 8.164689 12.386997 11.999694

aggregate()函数详解见下方

方法3.重复基因的表达量取最大值
代码语言:javascript复制
> library(tibble)
> library(dplyr)
> colnames(exp)[1] <- 'ID'
> exp_new <- merge(exp,soft,by = 'ID')
> exp_new <- exp_new[,c('GeneName',paste0('S',1:6))]
> dim(exp_new)
[1] 1000    7

# 方法1:aggregate()
> exp_max <- aggregate(x = exp_new[,colnames(exp_new)!= 'GeneName'], 
                    by = list(exp_new$GeneName), 
                    FUN = max) # FUN=max即可
> dim(exp_max)
[1] 946   7
> rownames(exp_max) <- exp_max$Group.1
> exp_max$Group.1 <- NULL
> exp_max[1:5,]
                S1        S2        S3        S4        S5        S6
15E1.2    5.826481  6.609382  6.807679  7.814539  8.950446  5.897408
AB016902  9.543493  9.385374  5.995681  8.729387 12.121148  9.337484
ABCA1    12.187228 14.106784 10.329787 11.635219  9.869883 12.745335
ABCC6     9.228167  9.337685  8.517037  6.703833 10.921183 10.143636
ACBD3    10.042010  8.634697 14.158641  8.164689 12.386997 11.999694

# 方法2:表达量排序后取最大值
> exp_max2 <- exp_new
> index <- order(rowMeans(exp_max2[,-ncol(exp_max2)]),decreasing = T) #取行平均值进行降序排序
> exp_max2 <- exp_max2[index,] #按行平均值对表达矩阵的进行排序
> exp_max2 <- exp_max2[!duplicated(exp_max2$GeneName),] #删除重复基因名,保留平均表达量最高的
> rownames(exp_max2) <- exp_max2$GeneName
> exp_max2$GeneName <- NULL
> exp_max2[1:5,]
                S1        S2       S3        S4        S5        S6
DCUN1D1  14.131177  6.638108 15.62183 15.459357 13.845734 12.301768
JPH3     13.767640 13.188044 14.54098 12.387238 13.578091  8.735123
SLC25A37  7.319557 11.190795 14.07802  9.833991 14.866505 13.296882
KIF5A    11.617515 17.033012 12.34930  6.605227 11.582761  9.196236
PCNXL2   12.426430 14.238976 11.85402 13.616833  3.654022 12.426465

# 方法3:使用dplyr包(思路与上述相近)
exp_max3 <- exp %>% 
  #合并探针的信息
  inner_join(anno,by="ID") %>% 
  #去掉多余信息,select支持按列名和列号同时选择
  select(c(GeneName,2:7)) %>%  
  #·增加一列,内容为每一行的平均数
  mutate(rowMean =rowMeans(.[,-1])) %>% 
  #把表达量的平均值按从大到小排序
  arrange(desc(rowMean)) %>% 
  # 去重,GeneName留下第一个
  distinct(GeneName,.keep_all = T) %>% 
  #GeneName转换为行名
  column_to_rownames(var="GeneName") %>% 
  #反向选择去掉平均值的那一列
  select(-rowMean)

# 管道符后面的.可以代表管道符前面传入的数据,如果调用tidyverse的函数应该都是可以省略的,
# 默认第一个参数,如果调用其他函数,用.代替就行。

❝方法3参考:超纲练习题不超纲 ❞

【补充】aggregate()函数
1. 基本语法
代码语言:javascript复制
aggregate(x = any_data, by = group_list, FUN = any_function)
# x: 进行运算的数据
# by: 进行运算的分组(以list形式)
# FUN: 执行运算的函数
2. 基本用法
代码语言:javascript复制
data <- data.frame(x1 = 1:5,                                  # Create example data
                   x2 = 2:6,
                   x3 = 1,
                   group = c("A", "A", "B", "C", "C"))
data                                                          # Print data
#   x1 x2 x3 group
# 1  1  2  1     A
# 2  2  3  1     A
# 3  3  4  1     B
# 4  4  5  1     C
# 5  5  6  1     C
代码语言:javascript复制
aggregate(x = data[ , colnames(data) != "group"],             
# Mean by group
          by = list(data$group),
          FUN = mean)
 
#   Group.1  x1  x2 x3        # by指定的列会被清除,生成新的列(Group.1)
# 1       A 1.5 2.5  1
# 2       B 3.0 4.0  1
# 3       C 4.5 5.5  1
代码语言:javascript复制
aggregate(x = data[ , colnames(data) != "group"],             # Sum by group
          by = list(data$group),
          FUN = sum)
 
#   Group.1 x1 x2 x3
# 1       A  3  5  2
# 2       B  3  4  1
# 3       C  9 11  2
3. 如果含有NA值
代码语言:javascript复制
data_NA <- data                                               
# Create data containing NAs
data_NA$x1[2] <- NA
data_NA$x2[4] <- NA
data_NA     
                                                  
# Print data
#   x1 x2 x3 group
# 1  1  2  1     A
# 2 NA  3  1     A
# 3  3  4  1     B
# 4  4 NA  1     C
# 5  5  6  1     C
代码语言:javascript复制
aggregate(x = data_NA[ , colnames(data_NA) != "group"],       
# aggregate without na.rm
          by = list(data_NA$group),
          FUN = mean)
 
#   Group.1  x1  x2 x3
# 1       A  NA 2.5  1
# 2       B 3.0 4.0  1
# 3       C 4.5  NA  1
代码语言:javascript复制
aggregate(x = data_NA[ , colnames(data_NA) != "group"],       
# Using na.rm option
          by = list(data_NA$group),
          FUN = mean,
          na.rm = TRUE)

#   Group.1  x1  x2 x3
# 1       A 1.0 2.5  1
# 2       B 3.0 4.0  1
# 3       C 4.5 6.0  1

0 人点赞