TCGA数据库和GTEx数据库下载和表达矩阵整理

下载数据：

UCSC Xena (xenabrowser.net)依次点击：

• ​​​​​​GDC TCGA Breast Cancer (BRCA) (20 datasets)
• HTSeq - FPKM (n=1,217) GDC Hub
• downloadhttps://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-BRCA.htseq_fpkm.tsv.gz; Full metadata
• ID/Gene Mappinghttps://gdc-hub.s3.us-east-1.amazonaws.com/download/gencode.v22.annotation.gene.probeMap; Full metadata
• 返回上一级
• Phenotype (n=1,284) GDC Hub
• downloadhttps://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-BRCA.GDC_phenotype.tsv.gz; Full metadata
• 返回初始页，依次点击
• GTEX (11 datasets)
• TOIL RSEM fpkm (n=7,862) UCSC Toil RNA-seq Recompute
• downloadhttps://toil-xena-hub.s3.us-east-1.amazonaws.com/download/gtex_RSEM_gene_fpkm.gz; Full metadata
• ID/Gene Mappinghttps://toil-xena-hub.s3.us-east-1.amazonaws.com/download/probeMap%2Fgencode.v23.annotation.gene.probemap; Full metadata
• 返回上一级
• GTEX phenotype (n=9,783) UCSC Toil RNA-seq Recompute
• downloadhttps://toil-xena-hub.s3.us-east-1.amazonaws.com/download/GTEX_phenotype.gz; Full metadata

library(data.table)
library(dplyr)
library(tidyverse)
#读取乳腺癌临床信息
brca.phe=fread("TCGA-BRCA.GDC_phenotype.tsv.gz",header = T, sep = '\t',data.table = F)
#读取基因表达数据
brca.fkpm=fread("TCGA-BRCA.htseq_fpkm.tsv.gz",header = T, sep = '\t',data.table = F)
#读取探针信息，修改ensembl
brca.pro=fread("gencode.v22.annotation.gene.probeMap",header = T, sep = '\t',data.table = F)
#取前两列进行转换
brca.pro=brca.pro[,c(1,2)]
#将探针信息和表达数据进行合并
brca.fkpm.pro=merge(brca.pro,brca.fkpm,by.y ="Ensembl_ID",by.x = "id" )
#distinct函数去除重复基因。
#还有一种方法是用aggregate函数，将相同基因归为一组进行处理
brca.fkpm.pro=distinct(brca.fkpm.pro,gene,.keep_all = T)
#将基因名设置为行名
brca.fkpm.pro <- column_to_rownames(brca.fkpm.pro,"gene")
#根据临床信息把癌和癌旁区进行区分
View(brca.phe)
##通过查看临床信息，我们发现在sample_type.samples列中，Primary Tumor为癌，Solid Tissue Normal为癌旁
#把行名换成样本名
rownames(brca.phe)=brca.phe$submitter_id.samples
#系统查看组织数量
table(brca.phe$sample_type.samples)
brca.phe.t=filter(brca.phe,sample_type.samples=="Primary Tumor")#筛选肿瘤组织，1114个
brca.phe.n=filter(brca.phe,sample_type.samples=="Solid Tissue Normal")#筛选正常组织，162个
#临床数据和表达矩阵取交集
z1=intersect(rownames(brca.phe.t),colnames(brca.fkpm.pro))#1097个
z2=intersect(rownames(brca.phe.n),colnames(brca.fkpm.pro))#113个
#将交集的表达矩阵取出
brca.t=brca.fkpm.pro[,z1]
brca.n=brca.fkpm.pro[,z2]
colnames(brca.n)=paste0("N",1:113)#将正常组织表达矩阵样本名修改N
colnames(brca.t)=paste0("T",1:1097)#将肿瘤组织表达矩阵样本名修改T
#合并两个表达矩阵
brca.exp=merge(brca.n,brca.t,by.x = 0,by.y = 0)
#得到TCGA数据库的表达矩阵，前面为正常组织113个，后面为肿瘤组织1097个，行名为基因名，列名为样本名
brca.exp=column_to_rownames(brca.exp,"Row.names")

 
#读取gtex的表达矩阵
gtex.exp=fread("gtex_RSEM_gene_fpkm.gz",header = T, sep = '\t',data.table = F)
#读取gtex的临床样本注释信息
gtex.phe=fread("GTEX_phenotype.gz",header = T, sep = '\t',data.table = F)
#读取gtex探针信息
gtex.pro=fread("probeMap_gencode.v23.annotation.gene.probemap",header = T, sep = '\t',data.table = F)
#发现sample和id是共同的列
gtex.pro=gtex.pro[,c(1,2)]
gtex.exp=merge(gtex.exp,gtex.pro,by.x ="sample",by.y = "id" )
#把下划线都去掉，名字都不可以下划线开头
colnames(gtex.phe)=c("Sample","body_site_detail (SMTSD)","primary_site","gender","patient","cohort")
#查看组织各组织样本数
table(gtex.phe$primary_site)
#把乳腺癌的样本筛选出来
gtex.phe.s=filter(gtex.phe,primary_site=="Breast")#221个
rownames(gtex.phe.s)=gtex.phe.s$Sample
#表达矩阵和临床信息取交集
x1=intersect(rownames(gtex.phe.s),colnames(gtex.exp))#179个
gtex.s=gtex.exp[,c("sample",x1)]
#修改ensemble
gtex.exp=merge(gtex.pro,gtex.s,by.x ="id",by.y ="sample")
gtex.s1=distinct(gtex.exp,gene,.keep_all = T)
gtex.s2 <- column_to_rownames(gtex.s1,"gene")
#得到GTEx表达矩阵，一共发现179个乳腺组织
gtex.s2 =gtex.s2[,-1]
#从官网可以看到GTEx是按照log2(fpkm+0.001)处理的，TCGA是按照log2(fpkm+1)，在合并之前，先要把他们的处理方式变成一样。
gtex.s3=2^gtex.s2
gtex.s5=log2(gtex.s3-0.001+1)
colnames(gtex.s5)=paste0("G",1:179)

#取交集，看用gene合并还是ensembl合并,根据结果采用gene合并
length(intersect(gtex.exp$id,brca.fkpm.pro$id))#40597
length(intersect(gtex.exp$gene,rownames(brca.fkpm.pro)))#57793
#按照rownames取交集，有些负数是指数转换没有完全回去，但值都特别小，可以忽略
all.data=merge(gtex.s5,brca.exp,by.x = 0,by.y = 0)
all.data <- column_to_rownames(all.data,"Row.names")
fwrite(all.data,file = "stomach.cancer.gtex.tcga.txt",sep = "\t",row.names = T)
all.data1<- as.matrix(all.data)

下载KEGG通路：

GSEA (gsea-msigdb.org)依次点击：

• Downloads
• Individual Human Gene Set GMTs--> Human Collections page.
• KEGG subset of CP
  (browse 186 gene sets)-->Gene Symbols

GSVA算法：

基因集变异分析，是一种非参数的无监督分析方法，主要用来评估芯片核转录组的基因集富集结果。主要是通过将基因在不同样本间的表达矩阵转换为基因集在样本间的表达矩阵，从而来评估不同代谢通路在不同样本间是否富集。其实就是研究感兴趣的基因在不同样本间的差异，作为一种分析方法，主要是为了从生物信息学的角度去解释导致表型差异的原因。
library(GSVA)
library(GSEABase)
library(clusterProfiler)
library(GSEA)
kegggmt2 <- read.gmt("c2.cp.kegg.v2022.1.Hs.symbols.gmt")
kegg_list = split(kegggmt2$gene, kegggmt2$term)
kegg2 <- gsva(nromalized.data, kegg_list, kcdf="Gaussian",method = "gsva")
rownames(annotation_col)=colnames(nromalized.data)
pheatmap::pheatmap(kegg2[1:20,], #热图的数据
                   cluster_rows =T,#行聚类
                   cluster_cols =F,#列聚类，可以看出样本之间的区分度
                   annotation_col = annotation_col,
                   show_colnames=F,
                   scale = "row", #以行来标准化，这个功能很不错
                   color =colorRampPalette(c("#FF7744", "white","#AAAAAA","#0044BB"))(100))