往前一步是优秀，退后一步是懵懂
我们的生信入门班和数据挖掘线上直播课程已经有了三年多的历史，培养了一波又一波优秀的生信人才。前面提到R语言授课时的超纲练习题，已经分享过两位优秀学员的答案。
下面继续来看优秀学员Dr.luka的分享：
R语言超纲练习题
(生信技能树优秀学员Dr.luka)
- 数据挖掘（GEO,TCGA,单细胞）2022年6月场，快速了解一些生物信息学应用图表
- 生信入门课-2022年6月场，你的生物信息学第一课
❝最好的输入就是输出。这篇笔记大量参考了同学们的分享，并结合自己的思路进行修改和补充。 ❞
❝十分感谢 佳男同学 的分享笔记系统学习＋主动探索，是最舒适的入门学习方式！ 徐谦同学 的分享笔记超纲练习题不超纲 ❞
1.读取文件
exp <- read.csv('exp.csv') #第一次读入不要用row.names=1，防止报错
soft <- read.table("soft.txt",sep = "\t",header = T)
exp内容展示
soft内容展示
2.探针过滤
由于实际情况可能存在注释文件的探针(probe_id)与表达矩阵探针不能完全对应的情况，因此在进行基因名转换之前，需要把探针进行过滤，留下有效的探针
table(exp$X %in% soft$ID) #如果有FALSE则说明有不对应情况
dim(exp) #过滤前探针数
exp <- exp[exp$X %in% soft$ID,]
dim(exp) #过滤后探针数
3.删除重复的基因名，整理表达矩阵
方法1.直接删除重复基因，保留下标最小的
#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.重复基因的表达量取平均值
> 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.重复基因的表达量取最大值
> 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. 基本语法
aggregate(x = any_data, by = group_list, FUN = any_function)
# x: 进行运算的数据
# by: 进行运算的分组(以list形式)
# FUN: 执行运算的函数
2. 基本用法
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
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
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值
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
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
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
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2022-05-21，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录
相关产品与服务
云直播
云直播（Cloud Streaming Services，CSS）为您提供极速、稳定、专业的云端直播处理服务，根据业务的不同直播场景需求，云直播提供了标准直播、快直播、云导播台三种服务，分别针对大规模实时观看、超低延时直播、便捷云端导播的场景，配合腾讯云视立方·直播 SDK，为您提供一站式的音视频直播解决方案。