重复一篇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,用来存放所有的压缩文件。

dir.create("data_in_one")

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

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格式的，很复杂，需要专门的函数读取。

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 文件 复现这篇文章，发邮件找我申请

## [e-mail] 即可

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的数据也是分的很开，并不用强求细节，掌握处理数据和画图是关键所在。