《GEO数据挖掘课程》配套练习题
本次配套练习题,需要复现的文章是 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_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 sapiensBG200951PVT1
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基因的表达量信息,根据临床信息,分类汇总后绘图。
写到最后
如果你也想开启自己的生物信息学数据处理生涯,加入生信技能树小圈子,还等什么呢,赶快行动起来吧! 发邮件(jmzeng1314@163.com)给生信技能树创始人jimmy就有惊喜哦!当然了,不能是辣鸡或者骚扰邮件啦,带上自己的简历和想学习交流的诚心吧!