重复一篇Cell文献的PCA图

这天，接到了生信技能树创始人jimmy老师的一个任务，要重复一篇CELL文章中的一个图示：

完成任务历程有点坎坷……

看到图第一步去找到了原文，《Tumor Evolution and Drug Response in Patient-Derived Organoid Models of Bladder Cancer》，找到了图示的地方，在补充材料部分，有一些基本信息，介绍了数据的存储，GEO数据库中的GSE103990, 还有用到了TCGA数据库中的bladder cancer数据。

这是一张PCA图，之前没有接触过，所以去查了一些资料，我这里就不多介绍了，网上资料一大堆，不过看过一些资料后，了解了个大概，涉及到很多知识点，还得去好好研究一下……

• 这两好玩的算法（PCA,EFA）
• 一文看懂主成分分析

1

TCGA-数据下载

由于GEO数据之前挖掘过，所以这次从TCGA开始。最好的教程在《生信技能树》，这话一点不假，跟着做就对了，下载TCGA数据有好多种方法，本次我尝试了最原始的方法，直接从网站下载。

网址为【https://cancergenome.nih.gov/】。

打开是这个样子的。

网页中下部（改版了，和教程不太一样）。

点开是这个样子。

再点箭头所示，进去是TCGA“超市”，和淘宝一样要加购物车。

在这里选择要下载的数据选项，我要下载bladder cancer数据，就在“CASES”里找到bladder cancer，然后在下面选择合适的选项，再在“Files”中选好合适选项，最后选好后就生成了你需要的数据了。下面进入“购物车”下载数据。

先点箭头所示部位下载电脑系统相应的下载工具。

再回到“购物车”下载箭头对应文件。

这时候文件夹中有三个文件，然后解压压缩文件。

然后在此文件夹中直接按“shift“ 右键，会出现下图，点箭头部分会出现对话框。

在对话框中写入图中红线所示文字，等一会就会开始下载文件。

下载好后在文件夹中就会看到很多的文件夹

把这些下载的文件先复制在一个rawdata文件中，这些文件都是一个个独立的文件夹，还不能直接用，需要合成到一个文件中，后期操作需要在R中实现。

2

TCGA-数据处理

首先新建一个文件夹data_in_one,用来存放所有的压缩文件。

代码语言：javascript复制

dir.create("data_in_one")

用for循环来把这些文件放到一起。

代码语言：javascript复制

for (dirname in dir("rawdata/")){  
  ## 使用list.files函数找到rawdata里面单个文件夹下面的压缩文件
  file <- list.files(paste0(getwd(),"/rawdata/",dirname),pattern = "*.counts")  #找到对应文件夹中的内容，pattern可以是正则表达式
  ## 使用file.copy函数复制粘贴压缩文件到data_in_one
  file.copy(paste0(getwd(),"/rawdata/",dirname,"/",file),"data_in_one")  #复制内容到新的文件夹
}

所有的文件被复制到了新的文件夹。

接下来把数据读入R语言中，找出文件名对应的TCGA id。 这个对应关系在上次下载的metadata文件中，这个文件是json格式的，很复杂，需要专门的函数读取。

代码语言：javascript复制

metadata <- jsonlite::fromJSON("metadata.cart.2019-03-03.json")

我们再用for循环提取对应的两者对应关系。

代码语言：javascript复制

naid_df <- data.frame()
for (i in 1:nrow(metadata)){
  naid_df[i,1] <- metadata$file_name[i]
  naid_df[i,2] <- metadata$associated_entities[i][[1]]$entity_submitter_id
}
colnames(naid_df) <- c("filename","TCGA_id")

批量读取数据。

代码语言：javascript复制

test <- data.table::fread(paste0("data_in_one/",naid_df$filename[1]))
expr_df <- data.frame(matrix(NA,nrow(test),nrow(naid_df)))

for (i in 1:nrow(naid_df)) {
  print(i)
  expr_df[,i]= data.table::fread(paste0("data_in_one/",naid_df$filename[i]))[,2]
}

给读入的数据添加列名和基因名称，每一个文件读取时都对应了一个TCGA id。

代码语言：javascript复制

colnames(expr_df) <- naid_df$TCGA_id
gene_id <- data.table::fread(paste0("data_in_one/",naid_df$filename[1]))$V1
expr_df <- cbind(gene_id=gene_id,expr_df)

因为后5行是我们不需要的，所以去掉后5行。

代码语言：javascript复制

expr_df <- expr_df[1:(nrow(expr_df)-5),]

保持数据为Rdata格式。

代码语言：javascript复制

save(expr_df,file = "expr_df.Rdata")

去掉gene_id的点号。

代码语言：javascript复制

library(tidyr)
expr_df_nopoint <- expr_df %>% 
  tidyr::separate(gene_id,into = c("gene_id"),sep="\.") 

标准化和差异分析都是用Deseq2这个包来完成，文中也有介绍他们是用这个包来做的。 首先把样本名称变成数据框格式。

代码语言：javascript复制

metadata <- data.frame(TCGA_id =colnames(expr_df)[-1])

下面用table这个函数统计一下膀胱癌样本的分类。

代码语言：javascript复制

table(substring(metadata$TCGA_id,14,15))

01代表原发灶，11代表正常固体组织，教程里在这里是分组做的，现在就跟着往下做。

代码语言：javascript复制

sample <- ifelse(substring(metadata$TCGA_id,14,15)=="01","cancer","recur")
metadata$sample <- as.factor(sample)

这里必须强行分割

如果认真跟了我们TCGA教程就不会耗费如此长时间在下载膀胱癌RNA-seq表达矩阵上面，UCSC的XENA浏览器有现成的，见：送你一篇TCGA数据挖掘文章

继续看表演

构建dds对象。

代码语言：javascript复制

dds <-DESeqDataSetFromMatrix(countData=mycounts, 
                             colData=metadata, 
                             design=~sample,
                             tidy=TRUE)

然后是漫长的DEseq分析。

代码语言：javascript复制

dds <- DESeq(dds)

vst标准化,时间也很长。

代码语言：javascript复制

vsd <- vst(dds, blind = FALSE)

用Deseq2内置的主成分分析来看一下样本分布（这个和任务没有关系，只是看看结果如何）。

代码语言：javascript复制

plotPCA(vsd, "sample")

这图就其实很有问题了，normal和tumor几乎分不开，需要详细解读。

3

GEO数据

接下来是GEO数据库数据的下载分析了。 最开始还是按着技能树的视频及代码做了处理，但是在处理过程中就一直出错，这里就不赘述了。后来经jimmy老大指点了一下，因为这是RNAseq数据，所以不需要用之前处理芯片的方法，直接在网页下载counts数据就可以了。

但是下载下来还需要一些处理，在这里试了不少方法依然报错，所以最后还是请教了健明老师。 下面是健明老师提供的代码，“大神一出手就知有没有”这话一点不错，现在还在学习摸索中，希望早日能写出这样的代码。

代码语言：javascript复制

rm(list = ls())
options(stringsAsFactors = F)
library(DESeq2)
library(edgeR)
library(limma)
library(airway)
## 在GEO数据库现在这个文件
a=read.table('GSE103990_Normalized_counts.txt.gz',
             header = T,sep = 't',row.names = 1)
boxplot(a[,1])
exprSet=a
print(dim(exprSet))
exprSet=exprSet[apply(exprSet,1, 
                              function(x) sum(x>1) > floor(ncol(exprSet)/20)),]
print(dim(exprSet))
colnames(exprSet)

a=read.table('group.txt')
library(GEOquery)  
gset <- getGEO('GSE103990', destdir=".",
                 AnnotGPL = F,     
                 getGPL = F)  
gset[[1]]
pd=pData(gset[[1]])

pd=pd[match(colnames(exprSet),pd$description.1), ]

title=pd$title
colnames(exprSet)=title
library(stringr)
passages=as.numeric(str_split(pd[,12],':',simplify = T)[,2])
stage=str_split(pd[,13],':',simplify = T)[,2]
grade=str_split(pd[,14],':',simplify = T)[,2]
group_list=ifelse(grepl('tumor',title),'tumor','organoids')
table(group_list)

colD=data.frame(group=group_list,
                stage=stage,
                grade=grade,
                passages=passages )
rownames(colD)=colnames(exprSet)

save(exprSet,colD,file = 'input.Rdata')

if(F){
  colnames(exprSet)
  pheatmap::pheatmap(cor(exprSet)) 
  # 组内的样本的相似性应该是要高于组间的！
  pheatmap::pheatmap(cor(exprSet),
                     annotation_col = colD,
                     show_rownames = F,
                     filename = 'cor_all.png')
  dim(exprSet)
  exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
  dim(exprSet)

  exprSet=log(edgeR::cpm(exprSet) 1)
  dim(exprSet)
  exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
  dim(exprSet)
  M=cor(log2(exprSet 1)) 
  pheatmap::pheatmap(M,annotation_col = colD)
  pheatmap::pheatmap(M,
                     show_rownames = F,
                     annotation_col = colD,
                     filename = 'cor_top500.png')



  library(pheatmap)
  pheatmap(scale(cor(log2(exprSet 1))))

}

library(stringr)
exprSet[1:4,1:4]
boxplot(exprSet[,1])
rownames(exprSet)=str_split(rownames(exprSet),'_',simplify = T)[,1]

## 这个文件，就是UCSC的XENA数据库的表达矩阵被R处理了，非常简单。# 如果需要 TCGA-BLCA.counts.Rdata  文件 复现这篇文章，发邮件找我申请## jmzeng1314@1314.com 即可
load("file =TCGA-BLCA.counts.Rdata")
RNAseq_expr[1:4,1:4]
rownames(RNAseq_expr)=str_split(rownames(RNAseq_expr),'[.]',simplify = T)[,1]

gid=intersect(rownames(exprSet),rownames(RNAseq_expr))

dat=cbind(exprSet[gid,],RNAseq_expr[gid,])
group=c(colD$group,rep('TCGA',ncol(RNAseq_expr)))
tmp=data.frame(group=group)
rownames(tmp)=colnames(dat)

dat=log(edgeR::cpm(dat) 1)
dat_top1000=dat[names(sort(apply(dat, 1,mad),decreasing = T)[1:1000]),]
dim(dat_top1000)
M=cor(log2(dat_top1000 1)) 
pheatmap::pheatmap(M,show_rownames = F,show_colnames = F,
                   annotation_col=tmp)

library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")  
dat=t(dat)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
dat=as.data.frame(dat)#将matrix转换为data.frame
dat.pca <- PCA(dat,graph = FALSE)

fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = group, # color by groups
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)
ggsave('all_samples_PCA.png')

一张漂亮的图出现了，和原文中的图有点出入，因为大家挑选的基因不一样，但是展现出来的规律是一样的，TCGA的样本跟作者的数据区分的很好，而且organoids的数据也是分的很开，并不用强求细节，掌握处理数据和画图是关键所在。

4

写在最后

TCGA数据下载有好多种方法，其他方法在生信技能树公众号推文中都有介绍，以后都要尝试一遍，找到最好用的方法。经过实战演练才知道自己差距有多大，觉得看过视频看过教程就会了，实际应用中会遇到各种挑战，有时候你上百度、谷歌搜索都不一定会找到解决方法，我觉得只有扎扎实实听前辈的教导，好好看看5本以上的R语言书，才会知道大神们写的代码是什么意思，还要多实战，从简单到复杂，最重要的是多看生信技能树博客和微信推文，学习最新最好的技能。

数据库 sql 数据分析 r语言

0 人点赞