模仿GEO-TCGA数据库挖掘文章的全套分析流程

Posted on by

这是一个学徒作业，首发于生信技能树公众号，见：《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_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基因的表达量信息，根据临床信息，分类汇总后绘图。

    文末友情推荐

    要想真正入门生物信息学建议务必购买全套书籍，一点一滴攻克计算机基础知识，书单在：什么，生信入门全套书籍仅需160 。
    如果大家没有时间自行慢慢摸索着学习，可以考虑我们生信技能树官方举办的学习班：

  • 数据挖掘学习班第7期（线上直播3周，马拉松式陪伴，带你入门），原价4800的数据挖掘全套课程， 疫情期间半价即可抢购。
  • 生信爆款入门-第9期（线上直播4周，马拉松式陪伴，带你入门），原价9600的生信入门全套课程，疫情期间3.3折即可抢购。
    如果你课题涉及到转录组，欢迎添加一对一客服：详见：你还在花三五万做一个单细胞转录组吗？
    号外：生信技能树知识整理实习生招募，长期招募，也可以简单参与软件测评笔记撰写，开启你的分享人生！另外，：绝大部分生信技能树粉丝都没有机会加我微信，已经多次满了5000好友，所以我开通了一个微信好友，前100名添加我，仅需150元即可，3折优惠期机会不容错过哈。我的微信小号二维码在：0元，10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》
​

Comments are closed.