专栏名称: 生信技能树

生物信息学学习资料分析，常见数据格式及公共数据库资料分享。常见分析软件及流程，基因检测及癌症相关动态。

数据挖掘课程能带给你什么收获

生信技能树 · 公众号 · · 2024-02-26 09:32

正文

在数据挖掘课程研发的早期（2020开春，疫情期间），我们就有预感，它会成为医学生/临床医师首选技能提高课。诚意满满的3周马拉松式授课，第一期就收到了非常棒的评价，见：数据挖掘第一期学习反馈。转眼间4年过去了，总计40期的数据挖掘课程也顺利落幕，非常多的学员乐于跟大家分享自己参加学习班的笔记和心得体会，我们没办法一一精选和展示。如果你看完后觉得有必要，可以考虑一下我们生信技能树举办的最适合医学生/临床医师的数据挖掘课程。

真情实感流露最是动人，前面我们精选了很多学员的心得体会：

这样的分享形式已经成为了咱们马拉松授课后亮丽的风景线，本月继续！

以下是正文

1.R与Rstudio

主要是学习到了会创建project啊，之前不会如此高效整理自己的项目....都是直接复制粘贴代码进去，所以各种报错，唉。

然后自从开了别人的一个文件project之后，然后回不到最开始的界面了，也不知道之前在哪个文件夹...setwd（）都不知道set去哪里

嗯嗯嗯？竟然光标放在那一行的任意位置都行... 我蠢啊 ，虽然没什么区别，但是就像是知道了了不起的小技巧一样。快捷键是ctrl+enter

2.数据类型

确实，学习的过程中不管学习到什么，肯定的说出自己的结论，错了也还有改正的机会呀！

Alt 加 - 可以直接得到 <-，但是一般情况下其实还是可以用=代替。

额我突然想起来在某一次代码出问题时，改过这个数据框的名字问题，明天在茫茫数据中找找，可能还比较显眼，应该是红色的错误。

match的知识点在7-7-3数据结构介绍的17min

7-9-2可以再看几遍。

3.数据结构

两个不同的向量可以用cbind组合为矩阵，但是矩阵的长度相同，数据类型相同，所以下面x，y是不同的数据类型，是数据框咯？

x=c(1,2,3,4,5,6)
y=c("Y","S","D","C","B","N")
cbind(x,y)

（PS: 学员的推理是错的，虽然最开始，x和y确实数据类型不一样，理论上不可能存在一个矩阵里面，但是呢，用cbind起来后，x和y就都成为了字符，所以可以存在一个矩阵了。）

x=c(1,2,3,4,5,6)
y=c("Y","S","D","C","B","N")
data.frame(x,y)

如果打函数时，遇到比较长需要自动补齐的函数，在出现选项时，可以直接上下键来挑选，然后点Tab或者enter键自动补齐

老师推荐了一个everything的软件，找东西比电脑内带的快很多。

4.函数和R包

安装R包可能会出现”is not available (for R version 3.5.2)“

主要有下面三个原因：

1.包名写错
2.安装命令使用错误
本机的R语言版本与包所要求的版本不符（极少）

browseVignettes("ggplot2") 如果不知道某个包的使用方法，可以直接给出这样的代码来查看使用方法和示例，当然也是可以用?ggplot2和example(ggplot2)

5.进阶知识

如果有什么包装不上，什么空间原因，可能是包的版本更新了，直接去原始的包存放的地方删掉这个包，重新装就好了。

GEO挖掘流程里，在之前我自己通过东抄一点西抄一点的组成的完整的一篇代码，总觉得自己的分组信息可能出了点问题，但是经过自己多次验证和GEO2R验证以及这次跟小洁老师的代码做对比，发现结果还是没有问题的，让我对自己的数据分析更有信心了。

8.下游分析常见图表介绍

以前我以为我了解的很清楚了的地方现在一听原理，也许不是那么回事， 但是更清晰明了了 。

P.value经常性等于这种格式3.495e-12，后面的12竟然表示的是10的多少次方，厉害鸭。我居然连科学计数法都不知道，果然每个人都有自己的知识盲区。

然后讲了一个新的思路，就是可以取两个不同的疾病联合在一起取差异基因，

热图聚类样本和正常没有聚类到一起，可以画图的矩阵标准化一下？什么意思7-18-1 6.45

rm(list = ls())
options(stringsAsFactors = F)
library(tinyarray)
library(stringr)
gse = "GSE42872"
geo = geo_download(gse)
group_list = ifelse(str_detect(geo$pd$title,"Control"),"control","treat")
group_list = factor(group_list,levels = c("control","treat"))

find_anno(geo$gpl,install = T)
ids <- toTable(hugene10sttranscriptclusterSYMBOL)

try = get_deg_all(geo$exp,
                  group_list,
                  ids,
                  logFC_cutoff = 1,
                  entriz = F,
                  adjust = F,heat_id = 1
                  )
try$plots
en = quick_enrich(try$cgs$deg$diff$diffgenes)
en[[4]]
names(en)

然后小洁老师因为讲课时学员上课氛围良好，**而主动增加了授课时长。**讲解了ceRNA网路的数据挖掘，让我受益匪浅，我也想试着重复一下其它文章的图表。

自己随便找一篇文章复现图。

文章标题是： Comprehensive Analysis of Competing Endogenous RNA (ceRNA) Network Based on RNAs Differentially Expressed in Lung Adenocarcinoma Using The Cancer Genome Atlas (TCGA) Database

首先是下载数据

刚开始想用简便一点的办法，用TCGAbiolinks下载数据，飞快，但是卡在了expda这一步，重复了几次，dyplr也重新安装了，但是问题无法解决，一次运行超慢。遂用了复杂的办法，直接用gdc client下载，听课的时候没有跟着运行，自己跑起来一言难尽，也算是基本掌握了这个下载方式吧，其实也还行，就是下载数据比起上一个工具来，蜗牛一样慢。

Warning messages: 1: `select_()` is deprecated as of dplyr 0.7.0. Please use `select()` instead. This warning is displayed once every 8 hours. Call `lifecycle::last_warnings()` to see where this warning was generated.

我i还是不懂这个，明明没有选择miRNA数据，为什么作者说下载了594个样本，然后包含mRNA, lncRNA, and miRNA expression data？他用在线数据库做了CeRNA网络分析，这一点我又得去看下回放了，第一遍基本没起到什么作用，done.

下载完了数据，使用与作者一样的edgeR包取了差异基因，崩了，完全不一样，他“筛选了2551个差异表达的mRNA，包括2033个上调和518个下调，1359个差异表达的lncRNA，包括1185个上调和174个下调，99个差异表达的miRNA，包括78个上调和21个下调。”而我总共才3833个上调，772个下调（我猜测可能是因为上一步我取了594个样本数据的75％，而作者使用TCGA biolinks后并没有取这个值），不管了，继续往下，我还没有分类。

下面就是看看肿瘤样本和正常样本的差别怎么样，这好像有点差

差异基因后热图

热图和火山图拼到一起（画这个图卡死好几次，垃圾电脑口吐芬芳）

生存分析（年龄中竟然有好多缺失值？不过不是用年龄来分析，应该不要紧）

不，年龄中有缺失值还是影响了。那是1.删掉有缺失值的样本还是2.使缺失值在分组中也是缺失值呢？目前我两种方法好像都不会欸，找一找。

先是使用这个 meta3 = meta2[na.omit(meta2$age),]但是结果是少了30行确实数量是对的，但是并没有删除缺失值的行啊，删除了个什么？

table(is.na(meta2$age))
meta2 = meta2[complete.cases(meta2$age),]

ok，成功

生存分析出图

行吧，因为前面删除了30个数据，后面又出了问题，exprSet还是535列，但是这就要涉及到exprSet的列名的前12个和meta的行名要对应起来，并且寻找出我之前删掉了哪些确实的不匹配的.... 饶了我吧，我已经预感到很多东西都出了问题 ，但是如果不删除缺失值，那么年龄分组缺失值不就没有组？难道还单独分出一个组放缺失值吗？也不行。

但是这个匹配又删除怎么搞.....

exprSet1 = exprSet
colnames(exprSet1)
rownames(meta)


library(stringr)
tmp = str_split(colnames(exprSet1),'-01A',simplify = T)[,1]
str_count(tmp,"TCGA")
str_locate(tmp,"TCGA")[,1]
table(tmp)
length(tmp)
colnames(exprSet1) = tmp

exprSet2 = exprSet1
###这样为什么不行？删除了25个是什么
exprSet2 = exprSet2[,rownames(meta) %in% colnames(exprSet2)=="TRUE"]

exprSet2 = exprSet2[,colnames(exprSet2) %in% colnames(meta1)]

table(rownames(meta) %in% colnames(exprSet2))

？最后发现不管怎样删除都不行了，数量为什么会不对等呢。