问题还没解决的学习记录

一期数据挖掘-1

文章是lncRNAs PVT1 and HAR1A are prognosis biomarkers and indicate
therapy outcome for diffuse glioma patients

OD:Oligodendroglioma少突胶质瘤

OA:oligoastrocytoma 少星形细胞瘤

A:astrocytoma星形细胞瘤

BGM:glioblastoma恶性胶质瘤

摘要

1.通过基因表达量找出在肿瘤发生和不同肿瘤发展中的差异基因

2.通过TCGA和GEO数据库发现，PVT1和CYTOR上调时，HAR1A和MIAT肿瘤下调，是通过已有数据获得

3.PVT1和CYTOR高表达，HAR1A和MIAT低表达与Ki-67水平和TP53突变有关

4.通过生存分析和cox回归，PVT1高表达和HAR1A低表达的患者生存差

5.降低PVT1和升高HAR1A可以提供放化疗后患者的生存时间

思路

1.用GEO和TCGA现有数据得出：PVT1和CYTOR在GBM中比低级别LGG（lower grade gliomas）的表达量高，同时HAR1A和MIAT在GBM中的表达量又比低级别(lower grade gliomas)的表达量低，因此提出，这四种因子或许在肿瘤发生发展中具有重要作用。

image-20190808141305898 image-20190808142638911

2.从TCGA获得的数据，作者通过基因表达量（PVT1、CYTOR、HAR1A、MIAT）将GBM分成亚型，本文作者因此用TCGA中的数据，想看看这四种因子是否能也将gliblastoma subtyps进行亚型分型，并且得到结果为：CYTOR, HAR1A 和MIAT 的基因表达在这四种亚型中差异显著）

image-20190808152217766

3.找到了差异基因，那么哪些差异基因对病人生存分析有影响呢？

上一步作者得到了CYTOR, HAR1A andMIAT同样用TCGA和GEO中的数据集，用kaplan-meier生存分析+log-rank来研究这四种基因对病人生存分析的影响。

最后得到结果：PVT1和HAR1A对病人生存分析影响显著，同时在GSE43378中也有同样结果。综上，高表达PVT1和低表达HAR1A的病人生存差

image-20190808153918314

4.通过前面分析的结果用自己的数据集进行PCR验证

验证结果与前面数据相同

image-20190808163519354

5.用TCGA中的数据验证放化疗病人的差异基因的表达

低表达PVT1和高表达HAR1A的胶质瘤患者在放化疗中能够获得更好的疗效，同时也表明，那些即将放化疗的患者，查看体内PVT1和HAR1A的表达量可以在放化疗前评估疗效，避免不必要的放化疗，做到精准医疗。

image-20190808163935357

目标是复现文章中的火山图，热图

文中gse4290

先把表达矩阵做好

rm(list = ls()) 
options()$repos 
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$repos 
options()$BioC_mirror

rm(list = ls())#清空环境变量
options(stringsAsFactors = F)##字符不作为因子读入
#####数据下载
library(GEOquery)
f<-'GSE4290.Rdata'
####getGPL获得平台的注释信息，但下载速度会慢很多
####而且注释文件格式大多不如bioconductor包好用
if(!file.exists(f)){
  gset<-getGEO('GSE4290',destdir='.',
               AnnotGPL=F,
               getGPL=F)
  save(gset,file=f)
}
load('GSE4290.Rdata')
#数据提取
ex<- exprs(gset[[1]])
pd <- pData(gset[[1]])

save(ex,pd,file='ex_pd.Rdata')
group_list<-pd[,35]

############代码查看矩阵是否需要log（来自GEO2R）
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
LogC
ex <- log2(ex+1)

##下次希望不要再卡在这了，现在的group_list是pData中的一列，也就是说和表达矩阵的列名是对应的，这点很重要！
##只要知道符合条件的group_list的坐标，这里用grep函数提取出坐标，然后表达矩阵提取这些坐标（列）就得到了新的
##表达矩阵，然后把想要的表达矩阵用cbind合起来就可以变成一个是按分组信息来的也就是说，表达矩阵的列名（样本名）##和grouplist的顺序是一致的，因为表达矩阵就是按group_list提取的。我觉得这个思维太好了，因为完全不用再考虑##表达矩阵的样本名和分组信息的样本是不是对应的了，希望自己以后也能顺手的这么用起来。
exprSet<-ex
  n_expr = exprSet[ , grep( "non-tumor",         group_list )]
  g_expr = exprSet[ , grep( "glioblastoma",      group_list )]
  a_expr = exprSet[ , grep( "astrocytoma",       group_list )]
  o_expr = exprSet[ , grep( "oligodendroglioma", group_list )]
  

## 样本分组，新的表达矩阵只有normal和gbm样本
{
  exprSet = cbind( n_expr, g_expr )
  group_list = c(rep( 'normal', ncol( n_expr ) ), 
                 rep( 'gbm',    ncol( g_expr ) ) )
}

dim( exprSet )
exprSet[ 1:5, 1:5 ]
table( group_list )

save( exprSet, group_list, file = 'exprSet_by_group.Rdata')

接下来PCA，来看一下STHDA上面主成分分析的介绍：然后开始漫长的折腾过程，但是还是没有做出来PCA

我用之前的能做出来的晓娟老师的GEO模版里的PCA的代码跑

df.pca <- PCA(t(ex), graph = FALSE) #graph=FALSE不出图
fviz_pca_ind(df.pca,
             geom.ind = "point",
             col.ind = group_list, 
             addEllipses = TRUE, #不加圈圈起来因为样本少
             legend.title = "Groups"
             ）

提示让我用imputePCA

但是出来下面的提示，说让我用imputaPCA，搜索后知道这个imputePCA是missMDA包中的一个函数，于是需要然后下载missMDA包，并把其中的PCA换成imputePCA。那么什么是missMDA?如下图所示，就是一种matrix中含有缺失值，但是可以用这个包是这些缺失值对PCA的结果影像不大，这个函数imputePCA就是能把含有缺失值的matrix变成一个这些缺失值可以变成一个可以估算的值，然后仍然可以用PCA中的FactoMineR来计算。也可以参考：http://blog.sciencenet.cn/home.php?mod=space&uid=267448&do=blog&id=1182057

image-20190809073317096

上图网站：http://factominer.free.fr/missMDA/PCA.html

而missMDA中含有的函数如下：https://www.rdocumentation.org/packages/missMDA/versions/1.14

MIPCA	Multiple Imputation with PCA
TitanicNA	Categorical data set with missing values: Survival of passengers on the Titanic
Overimpute	Overimputation diagnostic plot
imputeMFA	Impute dataset with variables structured into groups of variables (groups of continuous or categorical variables)
imputeMultilevel	Impute a multilevel mixed dataset
geno	Genotype-environment data set with missing values
estim_ncpFAMD	Estimate the number of dimensions for the Factorial Analysis of Mixed Data by cross-validation
imputePCA	Impute dataset with PCA
imputeCA	Impute contingency table
missMDA-package	Handling missing values with/in multivariate data analysis (principal component methods)
vnf	Questionnaire done by 1232 individuals who answered 14 questions
prelim	Converts a dataset imputed by MIMCA, MIPCA or MIFAMD into a mids object
snorena	Characterization of people who snore
estim_ncpMCA	Estimate the number of dimensions for the Multiple Correspondence Analysis by cross-validation
estim_ncpMultilevel	Estimate the number of dimensions for the Multilevel PCA, multlevel MCA or Multilevel FAMD by cross-validation

##这个是一个关于missMDA的示例
data(orange)
#extim_ncpPCA是missMDA中的一个函数,Estimate the number of dimensions for the Factorial Analysis of Mixed Data by cross-validation
nb <- estim_ncpPCA(orange,ncp.min=0,ncp.max=5,method.cv="Kfold",nbsim=20,pNA=0.05)
res.comp <- imputePCA(orange,ncp=nb$ncp)
require(FactoMineR)
res.pca <- PCA(res.comp$completeObs)
resMI <- MIPCA(orange,ncp=nb$ncp)
plot(resMI)

orange image-20190809072549660 image-20190809074518107

观察上面的orange的纵坐标是分类，因此由于我的目标依然是看normal组和tumor组是否分的开，就像上面这张图一样，我就把样本名换成对应的分组，集normal和tumor，expr如下图。

image-20190809081239842

nb <- estim_ncpPCA(expr,ncp.min=0,ncp.max=5,method.cv="Kfold",nbsim=20,pNA=0.05)
res.comp <- imputePCA(expr,ncp=nb$ncp)
require(FactoMineR)
res.pca <- PCA(res.comp$completeObs)
resMI <- MIPCA(expr,ncp=nb$ncp)
plot(resMI)

事实上是我too yong too sample，做出来的图是下面这样的没这个

image

那么当时为什么没有继续用露露的代码呢？就是接着第一个代码块往下？是因为有个报错，我没有搜索到答案，往下看。

data = as.data.frame(t(exprSet))#这里为什么要转换过来？因为此时的列名是样本名，而我们要做的PCA是需要行名是样本名
data$group = group_list #data新建一列是group_list
png( 'pca_plot.png', res=80 )
autoplot(prcomp(data[,1:(ncol(data)-1)]), data = data, colour = 'group', 
          label=T,frame=T) + theme_bw()

image-20190809082132706

然后我想可能是前面某一步和露露的代码不一样，那就是用GEOquery下载的时候，因为露露下载是把平台信息GPL下载了，是为了注释，把那些没有对应上探针的去掉了，我想可能是矩阵不同，所以导致了上面的那个misssing values in "x"的报错。我没有下载平台的信息，我还是晓娟老师讲课时用的GEO模版，于时我觉定重下载一下，就是下面这样。

library(GEOquery)
f<-'GSE4290.Rdata'
####getGPL获得平台的注释信息，但下载速度会慢很多
####而且注释文件格式大多不如bioconductor包好用
if(!file.exists(f)){
  gset<-getGEO('GSE4290',destdir='.',
               AnnotGPL=F,
               getGPL=T)
  save(gset,file=f)
}
load('GSE4290.Rdata')
GPL = gset@featureData@data

可是让我百思不得其解的是，什么都一样，竟然出现下面的报错，我怎么就没有featureData槽了？下面的图里明明有的，而且和露露的一样代码，不可能有错，找不到任何原因，搜索也没搜索到。

image-20190809082835194 image-20190809083045926

所有我才从这个露露的做PCA的autoplot转而去之前的GEO模版里去用GEO模版里的df.PCA那个PCA的代码跑。

现在有个重点是，用的这个imputePCA是没有错的，因为咩有报错，只是我没有找到用哪个函数跑，并且我的分组信息就是两个，和missMDA的包给的示例中的orange矩阵是不太一样的，他有7、8个分组，另外，我的矩阵是不是合乎要求，同时对应哪个函数，官方文档pdf文件把上面的那几个function都列出来了，那么多，我也不能一个个看不是，最重要的是，么有搜索到现场的 missMDA用什么函数的结果，所以现在仅停留在获得表达矩阵和分组信息，PCA都还没做得，老大。