《GEO数据挖掘课程》配套练习题

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是《GEO数据挖掘课程》的配套练习题

我把3年前的收费视频课程：3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月，免费到B站:

• 这个课程超级棒，B站免费学习咯：https://m.bilibili.com/video/BV1dy4y1C7jz
• 配套代码在GitHub哈：https://github.com/jmzeng1314/GSE76275-TNBC
• TCGA数据库挖掘，代码在：https://github.com/jmzeng1314/TCGA_BRCA
• GTEx数据库挖掘，代码在：https://github.com/jmzeng1314/gtex_BRCA
• METABRIC数据库挖掘，代码在：https://github.com/jmzeng1314/METABRIC

然后马上就有了3千多学习量，而且有学员给出来了图文并茂版本万字笔记，让我非常感动！

扫描下面二维码马上就可以学习起来啦，笔记需要至少半个小时来阅读哦！

本次配套练习题，需要复现的文章是 lncRNAs PVT1 and HAR1A are prognosis biomarkers and indicate therapy outcome for diffuse glioma patients，请自行仔细研读，拿到数据处理流程。

问题1：表达芯片里面的lncRNA

如何从Affymetrix HG-U133 Plus 2.0 arrays挑出1950 probe sets (corresponding to 1303 lncRNA genes),不包括 pseudogenes 。（数量差不多就好了，不要纠结于精准的1950个探针。）

library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
length(unique(ids$symbol))
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
plot(table(sort(table(ids$symbol))))
# 可以看到很多基因都是有着多个对应的探针

参考：http://www.biotrainee.com/thread-626-1-1.html  找到每个基因的分类属性。去gencode数据库下载gtf文件后，使用R包读入进行分类。或者其它bioconductor解决方案。

LncRNA的分类：

分类依据是lncRNA在基因组上面的位置和与附件一些编码基因的位置关系来分类的）

• Intergenic lncRNA（这种lncRNA完全位于基因的间距的和旁侧的编码基因是没有重叠的）
• bidirectional lncRNA（这种lncRNA和编码基因mRNA的位置小于1kp,转录的方向是相反的）
• Intronic lncRNA（这种lncRNA位于编码基因内含子的区域，不与该内含子的外显子存在重叠）
• Antisense lncRNA（这种lncRNA与蛋白质编码基因的外显子有重叠区，但是转录方向相反）
• sense overlapping lncRNA（这种lncRNA与蛋白质编码基因的外显子有重叠区，而且转录方向相同）

lncRNA的数据库：

NCBI 、UCSC 、Ensembl 、Gencode 、Lncpedia 、lncrnadb 、NONCODE 、the lncRNAand disease  datebase

还没有公认的权威的数据库，需要用到不同的数据库,各种数据库各有优劣。

问题2：GEO数据集下载

下载GEO dataset (GSE4290) 数据集，使用GEOquery包或者其它。然后筛选  77个 glioblastoma samples and 23 non-tumor controls

问题3：差异分析并且绘制热图

在下载的GSE4290表达矩阵里面提取1303 lncRNA基因的热图。（数量级是1000即可，不要求精确数量）

问题4：绘制火山图

仅仅是对1303个 lncRNA基因的表达矩阵进行差异分析，并且根据 (|logFC| >1; P <0.05),  阈值来画火山图

问题5：R包和GPL的soft信息差异

比较hgu133plus2.db里面的基因的注释信息和https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570的数据的区别。

主要是更新时间问题，看：hgu133plus2_dbInfo()

1558290_a_at	Homo sapiens	BG200951	PVT1

问题6: 指定基因（这里是PVT1）画boxplot，多个分组

下面仅仅是示例代码而已：

n_expr=raw_exprSet[,tmp1==" non-tumor"]
t_expr=raw_exprSet[,tmp1==" glioblastoma, grade 4"]
a=raw_exprSet[, grep('astrocytoma',tmp1)]
od=raw_exprSet[, grep('oligodendroglioma',tmp1)]

raw_exprSet=cbind(n_expr,t_expr,a,od)
group_list=c(rep('NC',ncol(n_expr)),
             rep('GBM',ncol(t_expr)),
             rep('a',ncol(a)),
             rep('od',ncol(od)))
raw_exprSet=log(raw_exprSet)
library(ggstatsplot)
dat=as.data.frame(t(raw_exprSet))
dat$group=group_list
# 1558290_a_at	Homo sapiens	BG200951	PVT1
ggbetweenstats(data = dat, x = group,  y = '1558290_a_at')

问题7：摸索TCGA的GBM的临床信息，找到classical, neural, proneural and mesenchymal分类方式。

建议使用ucsc的xena浏览器探索，下载数据。如果你网速太差，也可以看我这边备份的TCGA数据，就是来源于xena，ucsc的，都在，https://share.weiyun.com/5zLnKmO

需求最大的是tcga数据库的生存分析和表达量差异，看看这两个视频：

• https://www.bilibili.com/video/av25643438?p=9

• https://www.bilibili.com/video/av49363776?p=6

问题8：摸索TCGA的GBM的mRNA-seq表达矩阵（完成PVT1的boxplot）

下载mRNA-seq表达矩阵，counts矩阵，然后提取PVT1基因的表达量信息，根据临床信息，分类汇总后绘图。

写到最后

如果你也想开启自己的生物信息学数据处理生涯，加入生信技能树小圈子，还等什么呢，赶快行动起来吧！ 发邮件（[e-mail]）给生信技能树创始人jimmy就有惊喜哦！当然了，不能是辣鸡或者骚扰邮件啦，带上自己的简历和想学习交流的诚心吧！