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

生信技能树

发布于 2020-09-24 11:13:02

1K0

发布于 2020-09-24 11:13:02

文章被收录于专栏：生信技能树生信技能树

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

下面是《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')