####让报错变成英文 方便google
Sys.setenv(LANGUAGE = "en")
#禁止chr转成factor
options(stringsAsFactors = FALSE)
###清空环境
rm(list=ls())
##加载包 这个包是用来加快数据读取的
library(data.table)
library(dplyr)
library(tidyverse)
###读取胃癌的临床信息
stad.phe=fread("TCGA-STAD.GDC_phenotype.tsv",header = T, sep = '\t',data.table = F)
###查看一下数据类型
class(stad.phe)#"data.frame"
###读取表达谱信息 rnaseq的数据
stad.fkpm=fread("TCGA-STAD.htseq_fpkm.tsv",header = T, sep = '\t',data.table = F)
###查看列名 方便merger的合并列设置
colnames(stad.fkpm)
####读取探针信息 ，目的是为了将ensemble名字转为基因名
stad.pro=fread("gencode.v22.annotation.gene.probeMap",header = T, sep = '\t',data.table = F)
##查看一下数据
colnames(stad.pro)
###我们只要前两列进行转换
stad.pro=stad.pro[,c(1,2)]
colnames(stad.pro)
###用merge函数将探针转化的信息和表达谱信息进行合并
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 #把基因名转换为行名
stad.fkpm.pro=distinct(stad.fkpm.pro,gene,.keep_all = T)#去重复
dim(stad.fkpm.pro) #去重复以后的数目 58387 409
stad.fkpm.pro <- column_to_rownames(stad.fkpm.pro,"gene")#相当于把gene列移动 转换为行名
##此时，已构建好表达矩阵
###因为tcga数据中有癌和癌旁，所以我们先根据临床信息把癌和癌旁区分一下
View(stad.phe) #临床信息中行名与表达矩阵的列名相同
x=stad.phe$submitter_id.samples
rownames(stad.phe)=stad.phe$submitter_id.samples
##通过查看临床信息，我们发现在sample_type.samples列中，Primary Tumor为癌，Solid Tissue Normal为癌旁
x2=stad.phe$sample_type.samples
table(stad.phe$sample_type.samples)#癌组织443个 癌旁组织101个
##将癌组织和癌旁组织提取出来
stad.phe.t=filter(stad.phe,sample_type.samples=="Primary Tumor")
stad.phe.n=filter(stad.phe,sample_type.samples=="Solid Tissue Normal")
#样本信息中行数为544 而表达矩阵列数为408 说明有100多个样本没有表达矩阵
#看一下临床信息与实际表达矩阵的交集
z1=intersect(rownames(stad.phe.t),colnames(stad.fkpm.pro))#肿瘤的临床信息与表达矩阵的交集 375个
z2=intersect(rownames(stad.phe.n),colnames(stad.fkpm.pro))#癌旁的临床信息与表达矩阵的交集 32个 说明很多组织测序不成功
stad.t=stad.fkpm.pro[,z1]##所有癌组织的表达矩阵
stad.n=stad.fkpm.pro[,z2]##所有癌旁组织的表达矩阵
##改变一下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
#表达矩阵的数据完成
save(stad.exp,file='stas.exp.Rdata')
library(data.table)
###读取gtex的表达矩阵，注意解压和不解压都是可以读取的
##电脑内存小的话，会出现error。。。。
memory.limit(size=100000)##我也不知道科不科学。。#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的差异 一个是60498 一个是60483
###现在我们先合并gtex的基因信息
###合并之前先看一下列名 找到共同的合并列
colnames(gtex.pro)
colnames(gtex.exp)
gtex.pro=gtex.pro[,c(1,2)]
###我们发现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个交集 如果按照Ensembl却只有42566个，所以最后还是按照gene来合并
###现在要提取正常的胃组织的表达矩阵，我们要根据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")
###我们发现一共有209个胃组织
###我们从官网可以看到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)
###现在数据的处理方式都相同了 就有可比性了 将gtex的胃组织的数据与tcga的数据进行合并
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)##去除批次效应
nromalized.data=normalizeBetweenArrays(all.data)
as.data.frame()
?normalizeBetweenArrays
##http://www.gsea-msigdb.org/gsea/index.jsp GSEA网址
###############另外一种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(计算免疫浸润)
expr=as.matrix(expr)
kegg2 <- gsva(nromalized.data, kegg_list, kcdf="Gaussian",method = "gsva",parallel.sz=0)
#mx.diff False:正态分布（差异度不明显）（ 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) #转置以后每一列代表一个通路
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)
#把每个样本对应的通路的基因集都计算出来了
? 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, #热图的数据
cluster_rows = F,#行聚类
cluster_cols =F,#列聚类，可以看出样本之间的区分度
annotation_col = annotation_col,
show_colnames=F,
scale = "row", #以行来标准化，这个功能很不错
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相关系数
网友评论