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

健明 2021-07-14

展开全文

本次配套练习题，需要复现的文章是 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./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./geo/query/acc.cgi?acc=GPL570的数据的区别。

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

1558290_a_atHomo sapiensBG200951PVT1

问题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_atHomo sapiensBG200951PVT1ggbetweenstats(data = dat, x = group, y = '1558290_a_at')