内容如下:
1.外泌体和肝癌TCGA数据下载
2.数据格式整理
3.差异表达基因筛选
4.预后相关外泌体基因确定
5.拷贝数变异及突变图谱
6.外泌体基因功能注释
7.LASSO回归筛选外泌体预后模型
8.预后模型验证
9.预后模型鲁棒性分析
10.独立预后因素分析及与临床的相关性分析
11.列线图,ROC曲线,校准曲线,DCA曲线
12.外部数据集验证
13.外泌体模型与免疫的关系
14.外泌体模型与单细胞测序
############################## 03.差异表达基因筛选 #############################
上次讲了差异表达基因,通过正则表达式将肿瘤组和正常组提取出来,然后肿瘤组放前面,正常组放后面进行差异表达分析,下面讲一个万能的办法,不用正则表达式,就是自己生成一个表达,对样本ID进行分组。
如我们肝癌的样本,我们复制到一个新的表格里面,然后对每个样本进行分组标记:
生成这个表格以后,然后存成csv格式。
下面读取代码:
setwd("E:\\blog外泌体相关预测模型\\Figure 1")
dir()
dir()
data <- read.csv("LIHC_TPM.csv",header = T,sep = ",")
data[1:5,1:5]
annotation <- read.csv("gencode.v22.annotation.gene.probeMap",header = T,sep = "\t")
head(annotation)
match <- match(data$X,annotation$id)
head(match)
annotation <- annotation[match,]
head(annotation)
data[1:5,1:5]
identical(data$X,annotation$id)
data$X <- annotation$gene
data[1:5,1:5]# > data[1:5,1:5]
# X TCGA.DD.A4NG.01A TCGA.G3.AAV4.01A TCGA.2Y.A9H1.01A TCGA.BC.A10Y.01A
#1 RP11-368I23.2 0.000000 0.00000000 0.0000000 0.000000
#2 RP11-742D12.2 0.000000 0.03631085 0.0000000 0.000000
#3 RAB4B 1.549102 2.59344450 3.0610666 1.880890
#4 AC104183.2 0.000000 0.00000000 0.0000000 0.000000
#5 C12orf5 2.058289 1.82597669 0.9323438 1.505011
然后读取刚才生成的样本数据,匹配肿瘤和正常样本,进行差异表达分析:
rownames <- as.data.frame(data$X)
head(rownames)
names(rownames) <- "Symbol"sample <- read.csv("Sample_LIHC.csv",header = T,sep = ",")
head(sample)
tumorsample <- sample[sample$Type %in% "Tumor",] ### 提取肿瘤样本名
head(tumorsample)
normalsample <- sample[sample$Type %in% "Normal",] ### 提取正常样本名
head(normalsample)
tumor <- data[,colnames(data) %in% tumorsample$SampleID]
dim(tumor)[2]## 记住肿瘤样本数量
normal <- data[,colnames(data) %in% normalsample$SampleID]
dim(normal)[2]## 记住正常样本数量
## 合并数据
data <- cbind(rownames,tumor,normal)
data[1:5,1:5]### 然后就可以进行后续的差异表达分析colnames(data)
library(limma)rt=as.matrix(data)
rownames(rt) <- data[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data33=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
## View(data33)library(limma)
dim(data33)
library(limma)
dim(data33)
data33=avereps(data33)
dim(data33)
data33=data33[rowMeans(data33)>0,]
data33[1:5,1:5]data <- normalizeBetweenArrays(data33)data <- as.data.frame(data)
data[1:5,1:5]
colnames(data)tumorNum <- dim(tumor)[2]
NormalNum <- dim(normal)[2]modType=c(rep("tumor",tumorNum),rep("normal",NormalNum))### 样本量需要改
design <- model.matrix(~0+factor(modType))
design
colnames(design) <- c("con","treat")
fit <- lmFit(data,design)
cont.matrix<-makeContrasts(treat-con,levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)allDiff=topTable(fit2,adjust='fdr',number=200000)write.csv(allDiff,file="TCGA_LIHC_diff_expression2.csv")
下一节讲预后相关外泌体基因确定。