【原】善用公共数据库，从TCGA的乳腺癌多组学数据说起

健明 2021-07-14

展开全文

前天是全网第一个单细胞转录组视频课程的收官阶段，第四单元我们复现了单细胞转录组数据的结论是如何使用TCGA等公共数据库来验证的，这里值得重点分享一下。当然，使用TCGA数据库的前提是了解它！
如果你不想错过我们的精彩教程，请置顶我们：没看到通知？是不是五行缺星？
如果你不想漏掉我们往期教程，请学会搜索：历史宝藏这样

首先在UCSC Xena数据库下载TCGA计划的所有BRCA相关数据分析结果来进行下游挖掘。

UCSC Xena网址：https:///datapages/

我这里选择的是TCGA Breast Cancer (BRCA) (30 datasets) 而不是 GDC TCGA Breast Cancer (BRCA) (18 datasets) 一定要搞清楚哦！！！

下载的数据包括：

AgilentG4502A_07_3 (n=597) TCGA hub
IlluminaHiSeq (n=1,218) TCGA hub
gistic2 thresholded (n=1,080) TCGA hub
Phenotypes (n=1,247) TCGA hub

首先对芯片表达矩阵分析

这里仅仅是跑了PAM50分类，结果如下：

如果你感兴趣其它分析，可以看我安排给2018年学徒的数据挖掘任务，比如下载乳腺癌的芯片表达数据进行差异分析 https://mp.weixin.qq.com/s/CJb27qhbjdZadJDnK2vNLw

然后是针对RNA-seq表达矩阵的

因为GitHub容量限制，我仅仅是挑选了TNBC的病人，代码如下：

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('TCGA-BRCA.survival.tsv.gz',header = T,sep = '\t')
a=read.table('TCGA-BRCA.GDC_phenotype.tsv.gz',header = T,sep = '\t',quote = '')
(tmp=as.data.frame(colnames(a)))
tmp=a[,grepl('her2',colnames(a))]
table(a$breast_carcinoma_estrogen_receptor_status)
table(a$breast_carcinoma_progesterone_receptor_status)
table(a$lab_proc_her2_neu_immunohistochemistry_receptor_status)
eph=a[,grepl('receptor_status',colnames(a))]
eph=eph[,1:3]
## 挑选全部是阴性的
tnbc_s=a[apply(eph,1, function(x) sum(x=='Negative'))==3,1]
tnbc_s
save(tnbc_s,file = 'tnbc_s.Rdata')

然后在TNBC病人里面挑选那些既有normal又有tumor的样本，这样就只有9个TNBC病人了，他们的表达矩阵的主成分分析如下：

可以看到，normal和tumor在RNA-seq的表达水平上泾渭分明，就可以做差异分析流程啦，代码见：TCGA数据库中三阴性乳腺癌在亚洲人群中的差异表达， https://mp.weixin.qq.com/s/IOGfzzpcWkzyQPzMADKY4g

当然了，针对这么大的数据量，你可以任意开启自己的课题，比如我安排给2018年学徒的：

有PIK3CA基因突变的肿瘤病人的转录水平变化 https://mp.weixin.qq.com/s/MJLEZPWqzJe4LaKRDtiZQQ 从986个样本中挑出了327个有PIK3CA突变的样本，但是呢，同一个病人既有正常样本，又有肿瘤样本的表达信息 ，符合要求的样本就只有35个了
TP53突变型和TP53野生型BRCA病人的差异分析结果 https://mp.weixin.qq.com/s/Phu-MxA0d079HdtBWTHbWg 而且还有参考文章 https://www./content/10.1101/354779v1 可以比较容易把自己的代码跟已经发表的研究进行比较，发现 immune gene-sets had significantly higher enrichment levels in TP53-mutated BCs compared to TP53-wildtype BCs

somatic 突变分析

这里主要是下载MAF文件记录的突变数据，值得一提的是TCGAmutations包整合了TCGA中全部样本的maf文件

 # TCGAmutations包整合了TCGA中全部样本的maf文件
  # devtools::install_github(repo = "PoisonAlien/TCGAmutations")
  library(TCGAmutations)
  tmp=as.data.frame(tcga_available())

有趣的是，这些信息是基于hg19参考基因组的，但是突变与否的全景图跟参考基因组没关系，如下：

CNV分析

可以用来分组，然后查看分好组的样本的其它组学的表现情况。

临床资料分析

主要是生存分析啦，不过我前天分析的是IHC与PAM50分型的探索，需要对乳腺癌有一定了解才能看懂。乳腺癌的IHC分类和PAM50分型的差异情况

重点是使用TCGA数据辅助自己的科研

比如前面我们介绍的一个单细胞转录组文章，单细胞转录组探索CAFs的功能和空间异质性

把CAFs分成了4组，如下：

然后两个重点组别代表性基因可以拿到，如下：

去TCGA数据库的BRCA癌症里面的表达矩阵里面探索它的表现情况，代码如下：

library(stringr)
vCAF='Esam, Gng11, Higd1b, Cox4i2, Cygb, Gja4, Eng'
vCAF=unlist(str_split(vCAF,', ') )
mCAF='Dcn, Col12a1, Mmp2, Lum, Mrc2, Bicc1, Lrrc15, Mfap5, Col3A1, Mmp14, Spon1, Pdgfrl, Serpinf1, Lrp1, Gfpt2, Ctsk, Cdh11, Itgbl1, Col6a2, Postn, Ccdc80, Lox, Vcan, Col1a1, Fbn1, Col1a2, Pdpn, Col6a1, Fstl1, Col5a2, Aebp1'
mCAF=unlist(str_split(mCAF,', ') )

if(F){
  library(data.table)
  # 文件BRCA.htseq_counts.tsv.gz从UCSC的XENA数据库下载，大于100M所以不提供在这里。
  # a=fread('TCGA_BRCA/UCSC_xena/TCGA-BRCA.htseq_counts.tsv.gz',data.table=F)
  a[1:4,1:4]
  # Ensembl_ID TCGA-3C-AAAU-01A TCGA-3C-AALI-01A TCGA-3C-AALJ-01A
  # 1 ENSG00000000003.13         9.348728         8.714246        10.356452
  # 2  ENSG00000000005.5         1.584963         1.584963         5.727920
  # 3 ENSG00000000419.11        10.874981        10.834471        10.329796
  # 4 ENSG00000000457.12        10.121534        11.512247         8.867279
  library(org.Hs.eg.db)
  library(stringr)
  esid=str_split(a$Ensembl_ID,
                 '[.]',simplify = T)[,1]
  columns(org.Hs.eg.db)
  head(esid)
  rownames(a)=esid
  e2s=select(org.Hs.eg.db,keys = esid,columns = c( "ENSEMBL" ,  "SYMBOL" ),keytype = 'ENSEMBL')
  vCAF=toupper(vCAF);vCAF=vCAF[vCAF %in% e2s$SYMBOL]
  mCAF=toupper(mCAF);mCAF=mCAF[mCAF %in% e2s$SYMBOL]
  ng=e2s[match(c(vCAF,mCAF),e2s$SYMBOL),1]
  mat=a[ng,]
  mat=mat[,-1]
  save(mat,file = 'TCGA-BRCA-vCAF-and-mCAF-genes-expression.Rdata')
}

load(file = 'TCGA-BRCA-vCAF-and-mCAF-genes-expression.Rdata')
mat[1:4,1:4] 
M=cor(t(mat))
colnames(M)=c(vCAF,mCAF)
rownames(M)=c(vCAF,mCAF)
pheatmap::pheatmap(M)
pheatmap::pheatmap(M,cluster_rows = F,cluster_cols = F)
n=t(scale(t( M )))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap::pheatmap(n,cluster_rows = F,cluster_cols = F)