重复一篇Cell文献的PCA图
这天，接到了生信技能树创始人jimmy老师的一个任务，要重复一篇CELL文章中的一个图示：
完成任务历程有点坎坷……
看到图第一步去找到了原文，《Tumor Evolution and Drug Response in Patient-Derived Organoid Models of Bladder Cancer》，找到了图示的地方，在补充材料部分，有一些基本信息，介绍了数据的存储，GEO数据库中的GSE103990, 还有用到了TCGA数据库中的bladder cancer数据。
这是一张PCA图，之前没有接触过，所以去查了一些资料，我这里就不多介绍了，网上资料一大堆，不过看过一些资料后，了解了个大概，涉及到很多知识点，还得去好好研究一下……
1
TCGA-数据下载
由于GEO数据之前挖掘过，所以这次从TCGA开始。最好的教程在《生信技能树》，这话一点不假，跟着做就对了，下载TCGA数据有好多种方法，本次我尝试了最原始的方法，直接从网站下载。
网址为【https://cancergenome.nih.gov/】。
打开是这个样子的。
网页中下部（改版了，和教程不太一样）。
点开是这个样子。
再点箭头所示，进去是TCGA“超市”，和淘宝一样要加购物车。
在这里选择要下载的数据选项，我要下载bladder cancer数据，就在“CASES”里找到bladder cancer，然后在下面选择合适的选项，再在“Files”中选好合适选项，最后选好后就生成了你需要的数据了。下面进入“购物车”下载数据。
先点箭头所示部位下载电脑系统相应的下载工具。
再回到“购物车”下载箭头对应文件。
这时候文件夹中有三个文件，然后解压压缩文件。
然后在此文件夹中直接按“shift“+右键，会出现下图，点箭头部分会出现对话框。
在对话框中写入图中红线所示文字，等一会就会开始下载文件。
下载好后在文件夹中就会看到很多的文件夹
把这些下载的文件先复制在一个rawdata文件中，这些文件都是一个个独立的文件夹，还不能直接用，需要合成到一个文件中，后期操作需要在R中实现。
2
TCGA-数据处理
首先新建一个文件夹data_in_one,用来存放所有的压缩文件。
dir.create("data_in_one")
用for循环来把这些文件放到一起。
for (dirname in dir("rawdata/")){
file <- list.files(paste0(getwd(),"/rawdata/",dirname),pattern = "*.counts")
file.copy(paste0(getwd(),"/rawdata/",dirname,"/",file),"data_in_one")
}
所有的文件被复制到了新的文件夹。
接下来把数据读入R语言中，找出文件名对应的TCGA id。
这个对应关系在上次下载的metadata文件中，这个文件是json格式的，很复杂，需要专门的函数读取。
metadata <- jsonlite::fromJSON("metadata.cart.2019-03-03.json")
我们再用for循环提取对应的两者对应关系。
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")
批量读取数据。
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。
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行。
expr_df <- expr_df[1:(nrow(expr_df)-5),]
保持数据为Rdata格式。
save(expr_df,file = "expr_df.Rdata")
去掉gene_id的点号。
library(tidyr)
expr_df_nopoint <- expr_df %>%
tidyr::separate(gene_id,into = c("gene_id"),sep="\\.")
标准化和差异分析都是用Deseq2这个包来完成，文中也有介绍他们是用这个包来做的。
首先把样本名称变成数据框格式。
metadata <- data.frame(TCGA_id =colnames(expr_df)[-1])
下面用table这个函数统计一下膀胱癌样本的分类。
table(substring(metadata$TCGA_id,14,15))
01代表原发灶，11代表正常固体组织，教程里在这里是分组做的，现在就跟着往下做。
sample <- ifelse(substring(metadata$TCGA_id,14,15)=="01","cancer","recur")
metadata$sample <- as.factor(sample)
这里必须强行分割
如果认真跟了我们TCGA教程就不会耗费如此长时间在下载膀胱癌RNA-seq表达矩阵上面，UCSC的XENA浏览器有现成的，见：送你一篇TCGA数据挖掘文章
继续看表演
构建dds对象。
dds <-DESeqDataSetFromMatrix(countData=mycounts,
colData=metadata,
design=~sample,
tidy=TRUE)
然后是漫长的DEseq分析。
dds <- DESeq(dds)
vst标准化,时间也很长。
vsd <- vst(dds, blind = FALSE)
用Deseq2内置的主成分分析来看一下样本分布（这个和任务没有关系，只是看看结果如何）。
plotPCA(vsd, "sample")
这图就其实很有问题了，normal和tumor几乎分不开，需要详细解读。
3
GEO数据
接下来是GEO数据库数据的下载分析了。
最开始还是按着技能树的视频及代码做了处理，但是在处理过程中就一直出错，这里就不赘述了。后来经jimmy老大指点了一下，因为这是RNAseq数据，所以不需要用之前处理芯片的方法，直接在网页下载counts数据就可以了。
但是下载下来还需要一些处理，在这里试了不少方法依然报错，所以最后还是请教了健明老师。
下面是健明老师提供的代码，“大神一出手就知有没有”这话一点不错，现在还在学习摸索中，希望早日能写出这样的代码。
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语言书，才会知道大神们写的代码是什么意思，还要多实战，从简单到复杂，最重要的是多看生信技能树博客和微信推文，学习最新最好的技能。
如果你完全没有看懂我们在讲什么，那你可能需要下面的课程：
生信技能树（爆款入门培训课）巡讲第一站-重庆 （已结束）
生物信息学全国巡讲之粤港澳大湾区专场 （已结束）
生信技能树（爆款入门培训课）巡讲第二站-济南 （已结束）
生信技能树（爆款入门培训课）巡讲- 广州和上海（正在报名）