數(shù)據(jù)庫背景介紹
TCGA洁桌，全稱：The cancer genome altas
官網(wǎng)：https://cancergenome.nih.gov/
是由National Cancer Institute ( NCI, 美國國家癌癥研究所) 和 National Human Genome Research Institute (NHGRI, 國家人類基因組研究所) 合作建立的癌癥研究項(xiàng)目元莫，通過收集整理癌癥相關(guān)的各種組學(xué)數(shù)據(jù)妈拌。The Cancer Genome Atlas (TCGA) has quantified gene expression levels in >12000 samples from >33 cancer types.
TCGA一般測兩種數(shù)據(jù) 如：胃癌 一般測癌組織和癌旁組織
GTEx，全稱： The Genotype-Tissue Expression (GTEx) project
首次被提出來是2013年胶惰，上百位科學(xué)家聯(lián)名在Nature Genetics雜志發(fā)表的文章首次介紹了“基因型-組織表達(dá)工程”鹉戚，并成立了“基因型-組織表達(dá)研究聯(lián)盟”（Genotype-Tissue Expression Consortium，GTEx）以下簡稱“GTEx”）师幕。The GTEx has catalogued gene expression in >9,000 samples across 53 tissues from 544 healthy individuals.
但是通常在UCSC上面下載癌癥組織的測序數(shù)據(jù)
兩個(gè)小問題
問題一：既然癌癥組織的測序數(shù)據(jù)在TCGA粟按，為什么不在TCGA官網(wǎng)下載數(shù)據(jù)，而是在UCSC上面下載呢霹粥？
1灭将、 UCSC整合了多個(gè)癌癥公共數(shù)據(jù)庫的資源，所以數(shù)據(jù)整理起來很方便（如果在TCGA官網(wǎng)下載就很麻煩后控，首先是下載的時(shí)候麻煩庙曙，第二是下載之后整合數(shù)據(jù)沒法。
2浩淘、TCGA是在美國的官方網(wǎng)站矾利，正常情況下我們連Google都打不開，所以……………….
3馋袜、 UCSC在生信分析方面搭建了很多平臺男旗，包括以后看你會接觸到的UCSC基因組瀏覽器等等，所以UCSC下載數(shù)據(jù)是國際認(rèn)可的欣鳖，發(fā)文章的時(shí)候不會有reviewer質(zhì)疑察皇。
問題二：既然TCGA有癌旁的數(shù)據(jù)，為什么還要下載GTEx的正常組織數(shù)據(jù)呢泽台？
1什荣、 不是所有的癌腫都有癌旁數(shù)據(jù)的，比如骨肉瘤這些就沒有癌旁
2怀酷、最重要的原因是TCGA中的癌旁數(shù)據(jù)太少了稻爬，比如胃癌只有100多癌旁，癌組織樣本卻有400多個(gè)蜕依。
3桅锄、有時(shí)候也是為了作圖漂亮
數(shù)據(jù)下載
從UCSC Xena官網(wǎng)下載TCGA數(shù)據(jù)庫中胃癌RNA-Seq data和GTEX數(shù)據(jù)庫中的正常組織RNA-Seq data琉雳。
進(jìn)入U(xiǎn)CSC官網(wǎng)網(wǎng)址：https://xenabrowser.net/
點(diǎn)擊DATA SETS：https://xenabrowser.net/datapages/
選擇GDC TCGA Stomach Cancer (STAD) (15 datasets)：https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Stomach%20Cancer%20(STAD)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443
下載TCGA數(shù)據(jù)庫中胃癌RNA-Seq data
進(jìn)入gene expression RNAseq中的HTSeq - FPKM (n=407) GDC Hub：
怎么分辨癌組織和癌旁組織呢？這時(shí)候需要下載臨床信息
進(jìn)入phenotype中的Phenotype (n=544) GDC Hub：
進(jìn)入phenotype中的survival data (n=502) GDC Hub
下載GTEX數(shù)據(jù)庫中的正常組織RNA-Seq data
數(shù)據(jù)的整合
####讓報(bào)錯(cuò)變成英文 方便google
Sys.setenv(LANGUAGE = "en")
#禁止chr轉(zhuǎn)成factor
options(stringsAsFactors = FALSE)
###清空環(huán)境
rm(list=ls())
##加載包 這個(gè)包是用來加快數(shù)據(jù)讀取的
library(data.table)
library(dplyr)
library(tidyverse)
###讀取胃癌的臨床信息
stad.phe=fread("TCGA-STAD.GDC_phenotype.tsv",header = T, sep = '\t',data.table = F)
###查看一下數(shù)據(jù)類型
class(stad.phe)#"data.frame"
###讀取表達(dá)譜信息 rnaseq的數(shù)據(jù)
stad.fkpm=fread("TCGA-STAD.htseq_fpkm.tsv",header = T, sep = '\t',data.table = F)
###查看列名 方便merger的合并列設(shè)置
colnames(stad.fkpm)
####讀取探針信息 友瘤，目的是為了將ensemble名字轉(zhuǎn)為基因名
stad.pro=fread("gencode.v22.annotation.gene.probeMap",header = T, sep = '\t',data.table = F)
##查看一下數(shù)據(jù)
colnames(stad.pro)
###我們只要前兩列進(jìn)行轉(zhuǎn)換
stad.pro=stad.pro[,c(1,2)]
colnames(stad.pro)
###用merge函數(shù)將探針轉(zhuǎn)化的信息和表達(dá)譜信息進(jìn)行合并
stad.fkpm.pro=merge(stad.pro,stad.fkpm,by.y ="Ensembl_ID",by.x = "id" )
dim(stad.fkpm.pro) # 60483 409
rownames(stad.fkpm.pro)=stad.fkpm.pro$gene #把基因名轉(zhuǎn)換為行名
stad.fkpm.pro=distinct(stad.fkpm.pro,gene,.keep_all = T)#去重復(fù)
dim(stad.fkpm.pro) #去重復(fù)以后的數(shù)目 58387 409
stad.fkpm.pro <- column_to_rownames(stad.fkpm.pro,"gene")#相當(dāng)于把gene列移動 轉(zhuǎn)換為行名
##此時(shí)翠肘，已構(gòu)建好表達(dá)矩陣
###因?yàn)閠cga數(shù)據(jù)中有癌和癌旁，所以我們先根據(jù)臨床信息把癌和癌旁區(qū)分一下
View(stad.phe) #臨床信息中行名與表達(dá)矩陣的列名相同
x=stad.phe$submitter_id.samples
rownames(stad.phe)=stad.phe$submitter_id.samples
##通過查看臨床信息辫秧，我們發(fā)現(xiàn)在sample_type.samples列中束倍，Primary Tumor為癌，Solid Tissue Normal為癌旁
x2=stad.phe$sample_type.samples
table(stad.phe$sample_type.samples)#癌組織443個(gè) 癌旁組織101個(gè)
##將癌組織和癌旁組織提取出來
stad.phe.t=filter(stad.phe,sample_type.samples=="Primary Tumor")
stad.phe.n=filter(stad.phe,sample_type.samples=="Solid Tissue Normal")
#樣本信息中行數(shù)為544 而表達(dá)矩陣列數(shù)為408 說明有100多個(gè)樣本沒有表達(dá)矩陣
#看一下臨床信息與實(shí)際表達(dá)矩陣的交集
z1=intersect(rownames(stad.phe.t),colnames(stad.fkpm.pro))#腫瘤的臨床信息與表達(dá)矩陣的交集 375個(gè)
z2=intersect(rownames(stad.phe.n),colnames(stad.fkpm.pro))#癌旁的臨床信息與表達(dá)矩陣的交集 32個(gè) 說明很多組織測序不成功
stad.t=stad.fkpm.pro[,z1]##所有癌組織的表達(dá)矩陣
stad.n=stad.fkpm.pro[,z2]##所有癌旁組織的表達(dá)矩陣
##改變一下stad.t和stad.n的列名 方便查看
colnames(stad.n)=paste0("N",1:32)#癌旁的列名重新命名
colnames(stad.t)=paste0("T",1:375)#腫瘤的列名重新命名
stad.exp=merge(stad.n,stad.t,by.x = 0,by.y = 0)##合并
colnames(stad.exp)
stad.exp <- column_to_rownames(stad.exp,"Row.names")#58387 407
#表達(dá)矩陣的數(shù)據(jù)完成
save(stad.exp,file='stas.exp.Rdata')
library(data.table)
###讀取gtex的表達(dá)矩陣盟戏，注意解壓和不解壓都是可以讀取的
##電腦內(nèi)存小的話绪妹，會出現(xiàn)error。柿究。邮旷。。
memory.limit(size=100000)##我也不知道科不科學(xué)笛求。廊移。#60498 7863
gtex.exp=fread("gtex_RSEM_gene_fpkm",header = T, sep = '\t',data.table = F)
save(gtex.exp,file='gtex.exp.Rdata')
###
dim(gtex.exp) #60498 7863
gtex.exp[1:5,1:5]
###讀取gtex的臨床樣本注釋信息
gtex.phe=fread("GTEX_phenotype",header = T, sep = '\t',data.table = F)
##查看一下
View(gtex.phe)
###讀取gtex的基因注釋信息 也就是探針信息
gtex.pro=fread("probeMap_gencode.v23.annotation.gene.probemap",header = T, sep = '\t',data.table = F)
dim(gtex.pro) #60498 6
####我們比較一下v23和v22的差異 一個(gè)是60498 一個(gè)是60483
###現(xiàn)在我們先合并gtex的基因信息
###合并之前先看一下列名 找到共同的合并列
colnames(gtex.pro)
colnames(gtex.exp)
gtex.pro=gtex.pro[,c(1,2)]
###我們發(fā)現(xiàn)sample和id是共同的列
head(gtex.pro)
head(gtex.exp)
gtex.exp[1:4,1:4]
gtex.fkpm.pro=merge(gtex.pro,gtex.exp,by.y ="sample",by.x = "id" )#60498 7864
###取一波交集 目的是為了決定之后的gtex和stad合并是按照symbol還是enseble來合并
length(intersect(gtex.pro$id,stad.pro$id)) #42566
length(intersect(rownames(stad.exp),gtex.fkpm.pro$gene)) #57993
###我們可以看到 如果按照基因合并會有57793個(gè)交集 如果按照Ensembl卻只有42566個(gè)，所以最后還是按照gene來合并
###現(xiàn)在要提取正常的胃組織的表達(dá)矩陣探入，我們要根據(jù)gtex的臨床信息來匹配胃組織的sample
colnames(gtex.phe)
rownames(gtex.phe)=gtex.phe$Sample
table(gtex.phe$_primary_site)
colnames(gtex.phe)=c("Sample","body_site_detail (SMTSD)","primary_site","gender","patient","cohort")
colnames(gtex.phe)
table(gtex.phe$primary_site)#stomach 209
gtex.phe.s=filter(gtex.phe,primary_site=="Stomach")
x1=intersect(rownames(gtex.phe.s),colnames(gtex.fkpm.pro))
gtex.s=gtex.fkpm.pro[,c("gene",x1)]
rownames(gtex.s) <- gtex.s$gene
gtex.s1 <- distinct(gtex.s,gene,.keep_all = T)
gtex.s2 <- column_to_rownames(gtex.s1,"gene")
###我們發(fā)現(xiàn)一共有209個(gè)胃組織
###我們從官網(wǎng)可以看到gtex是按照log2(fpkm+0.001)處理的 stad是按照log2(fpkm+1)狡孔，所以我們在合并之前 先要把他們的處理方式變成一樣。
gtex.s3=2^gtex.s2
log2(0.001) #-9.965784
gtex.s5=log2(gtex.s3-0.001+1)
colnames(gtex.s5)=paste0("G",1:174)
###現(xiàn)在數(shù)據(jù)的處理方式都相同了 就有可比性了 將gtex的胃組織的數(shù)據(jù)與tcga的數(shù)據(jù)進(jìn)行合并
all.data=merge(gtex.s5,stad.exp,by= 0) #57793 582
all.data <- column_to_rownames(all.data,"Row.names")
save(all.data,file = 'all.data.Rdata')
我也暈了蜂嗽。苗膝。。
library(limma)##去除批次效應(yīng)
nromalized.data=normalizeBetweenArrays(all.data)
as.data.frame()
?normalizeBetweenArrays
##http://www.gsea-msigdb.org/gsea/index.jsp GSEA網(wǎng)址
###############另外一種gmt
BiocManager::install("GSEABase")
library(GSEABase)
library(clusterProfiler)
library("devtools")
install_github("GSEA-MSigDB/GSEA_R")
library(GSEA)
library(dplyr)
kegggmt2 <- read.gmt("c2.cp.kegg.v7.4.symbols.gmt")#12797 2
kegg_list = split(kegggmt2$gene, kegggmt2$term)
library(GSVA)
?gsva
#method:gsva;zscore;ssgsea(計(jì)算免疫浸潤)
expr=as.matrix(expr)
kegg2 <- gsva(nromalized.data, kegg_list, kcdf="Gaussian",method = "gsva",parallel.sz=0)
#mx.diff False:正態(tài)分布（差異度不明顯）（ True：雙峰分布（差異度更大植旧，更明顯辱揭。）
#########################自定義的基因集
gene.set=read.table("5.4.GENE.SET.txt",
header =F,sep = '\t',quote = '')
kegg.123=read.table("5.4.GENE.NAME.txt",
header =F,sep = '\t',quote = '')
gene.set1=as.matrix(gene.set)
gene.set2=t(gene.set1) #轉(zhuǎn)置以后每一列代表一個(gè)通路
gmt=list() #GSCA輸入需要list
for (i in 1:19) {
y=as.character(gene.set2[,i])
b=which(y=="")
gmt[[i]]=y[-b]
}
names(gmt)=kegg.123$V1
gmt=gmt[-19]
View(gmt)
getwd()
library(GSVA)
es.dif.nromalized <- gsva(nromalized.data, gmt, mx.diff=TRUE,kcdf="Poisson",parallel.sz=8)
es.max.nromalized <- gsva(nromalized.data, gmt, mx.diff=FALSE)
#把每個(gè)樣本對應(yīng)的通路的基因集都計(jì)算出來了
? data.frame
annotation_col = data.frame(
Tissuetype =c(rep("Stomach",174),rep("Solid Tissue Normal",32),rep("Tumor",375)),
Database =c(rep("gtex",174), rep("TCGA",407))
)
rownames(annotation_col)=colnames(es.max.nromalized)
pheatmap::pheatmap(es.max.nromalized, #熱圖的數(shù)據(jù)
cluster_rows = F,#行聚類
cluster_cols =F,#列聚類，可以看出樣本之間的區(qū)分度
annotation_col = annotation_col,
show_colnames=F,
scale = "row", #以行來標(biāo)準(zhǔn)化病附，這個(gè)功能很不錯(cuò)
color =colorRampPalette(c("green", "black","red"))(100))
###############特定基因分析
exprSet.all.r=all.data[c("RORC", "RORB", "RORA"),]
exprSet.all.r=t(exprSet.all.r)
exprSet.all.r=as.data.frame(exprSet.all.r)
x=c(rep("GTEX",174),rep("N",32),rep("T",375))
exprSet.all.r$Type=x
exprSet.rorc=exprSet.all.r[,c(1,4)]
exprSet.rorc$Gene=rep("RORC")
colnames(exprSet.rorc)[1]="Relative Expression"
exprSet.rorb=exprSet.all.r[,c(2,4)]
exprSet.rorb$Gene=rep("RORB")
colnames(exprSet.rorb)[1]="Relative Expression"
exprSet.rora=exprSet.all.r[,c(3,4)]
exprSet.rora$Gene=rep("RORA")
colnames(exprSet.rora)[1]="Relative Expression"
x.all=rbind(exprSet.rorc,exprSet.rorb,exprSet.rora)
colnames(x.all)
library(ggsignif)
library(ggpubr)
library(ggplot2)
p <- ggboxplot(x.all, x = "Gene", y = "Relative Expression",
color = "Type", palette = "Type",
add = "Type")
p + stat_compare_means(aes(group = Type))
table(x.all$Gene)
my_comparisons <- list(c("RORA","RORB"), c("RORA","RORC"),c("RORB", "RORC"))
p +geom_signif(comparisons = my_comparisons,
step_increase = 0.2,map_signif_level = F,
test = t.test,size=0.8,textsize =4)
?geom_signif
x.c.b=cbind(exprSet.rorc,exprSet.rorb)
GSEA
GSVA
GO/KEGG
colnames(x.c.b)=c("RORC","Type","Gene", "RORB","Type","Gene" )
x.c.b=x.c.b[,c(1,4)]
library(ggplot2)
library(ggpubr)
## Loading required package: magrittr
p1 <- ggplot(data = x.c.b, mapping = aes(x = RORC, y = RORB)) +
geom_point(colour = "red", size = 2) +
geom_smooth(method = lm, colour='blue', fill='gray') #添加擬合曲線
p1
p1 + stat_cor(method = "pearson", label.x = -0.4, label.y = 0.2) #添加pearson相關(guān)系數(shù)
