TCGA_RNA-seq_limma分析

该博客详细介绍了如何利用R包limma进行RNA-seq数据的差异表达分析。首先,从读取计数数据开始,然后进行数据过滤和标准化。接着,利用经验贝叶斯方法估计dispersion,并构建广义线性模型。在模型拟合后,执行差异表达分析,通过设置不同的阈值来筛选显著差异表达的基因。整个过程包括了数据预处理、模型建立、统计检验和结果筛选等多个关键步骤。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

limma分析的大致流程

  • 导入read count, 保存为专门的对象用于后续分析
  • 原始数据过滤,根据标准化read count 或者 raw count 作为筛选标准
  • raw read count 标准化
  • 通过各种算法(如经验贝叶斯,EM)预测dispersion离散值
  • 广义线性模型拟合数据
  • 差异分析,也就是统计检验部分
rm(list=ls())
library(limma)
genesymbol<-read.delim("/Users/Desktop/ESCC_RNA_count_data/ESCC_count_cancer_expression.txt",header = TRUE,check.names = F)
condition<-read.delim("/Users/Desktop/ESCC_RNA_count_data/clinical merge.txt",header=TRUE)
genesymbol<-avereps(genesymbol[,-1],ID=genesymbol$id)
condition<-condition[,c("submitter_id","alcohol_history")]
for(i in 1:ncol(genesymbol)){
  colnames(genesymbol)[i]<-substr(colnames(genesymbol)[i],1,12)
}
condition<-condition[!condition$alcohol_history=="not_reported",]
genesymbol_alcohol<-genesymbol[,colnames(genesymbol)%in%condition$submitter_id]

#将0,1转换成not,yes,将0,1转化为not与yes是为了构建差异比较矩阵
condition$alcohol_history_text=as.factor(ifelse(condition$alcohol_history==1,"Yes","Not"))
condition$alcohol_history<-NULL
condition[,2]<-as.factor(condition[,2])

#构建比较矩阵
design<-model.matrix(~0+condition$alcohol_history_text)
colnames(design)=levels(condition$alcohol_history_text)
rownames(design)=colnames(genesymbol_alcohol)

#一共需要三个矩阵,分别是表达矩阵,样本矩阵和差异比较矩阵
#创建差异比较矩阵

condition$alcohol_history_text=relevel(condition$alcohol_history_text,ref="Not")
constrasts=paste(rev(levels(condition$alcohol_history_text)),collapse = "-")
cont.matrix<-makeContrasts(contrasts=constrasts,levels=design)

#建模前需要进行归一化
#limma输入数据一定要进行log2转化
#一种方法是通过查看boxplot查看数据整齐与否再决定是否归一化
#另一种直接用voom方法进行归一化

dge<-DGEList(counts=genesymbol_alcohol,group=condition$alcohol_history_text)

#过滤基因表达量(至少在两个样本中表达量大于1)

keep <- rowSums(cpm(dge)>0.5) >= 2
table(keep)

#重置文库大小,标化数据

dge<-dge[keep,,keep.lib.sizes=FALSE]
dge<-calcNormFactors(dge) ##TMM方法进行标化

#差异分析

#使用“limma-voom"进行差异分析
v<-voom(dge,design,plot=TRUE)#会自动计算log(cpm)值
#拟合线性模型
fit<-lmFit(v,design)
#针对给定的对比计算估计系数和标准误差
fit3=contrasts.fit(fit,cont.matrix)
#计算出t统计量,F统计量和差异表达倍数的对数
fit3<-eBayes(fit3)
DEG3=topTable(fit3,coef=constrasts,n=Inf)

#不同阈值输出结果

logFC_cutoff_1<-c(1)
DEG3$change_1=as.factor(
  ifelse(DEG3$PValue<0.05&abs(DEG3$logFC)>logFC_cutoff_1,
         ifelse(DEG3$logFC>logFC_cutoff_1,"up","down"),"NOT")
)

logFC_cutoff_2<-c(2)
DEG3$change_2=as.factor(
  ifelse(DEG3$PValue<0.05&abs(DEG3$logFC)>logFC_cutoff_2,
         ifelse(DEG3$logFC>logFC_cutoff_2,"up","down"),"NOT")
)

logFC_cutoff_sd<-with(DEG3,mean(abs(logFC))+2*sd(abs(logFC)))
DEG3$change_sd=as.factor(
  ifelse(DEG3$PValue<0.05&abs(DEG3$logFC)>logFC_cutoff_sd,
         ifelse(DEG3$logFC>logFC_cutoff_sd,"up","down"),"NOT")
)
write.table(DEG3,'/Users/Desktop/ESCC_RNA_count_data/DE_limmaR.txt',quote = F,sep = '\t')
# 加载必要R包 suppressMessages({ library(data.table) library(tidyverse) library(ggsignif) library(RColorBrewer) library(limma) library(ggplot2) library(ggpubr) library(beepr) library(gplots) library(pheatmap) }) # 设置桌面路径(适用于Windows) desktop_path <- "C:/Users/86183/Desktop" # Windows 示例 # ------------------------- 模拟数据生成 ------------------------- set.seed(123) gene_ids <- paste0("Gene", 1:100) # 定义正常和肿瘤样本 normal_samples <- paste0("TCGA-XX-XXXX-11A", 1:5) tumor_samples <- paste0("TCGA-XX-XXXX-01A", 1:5) # 模拟FPKM数据 fpkm_data <- matrix(rlnorm(100*10, meanlog = 1, sdlog = 0.5), nrow = 100) fpkm_data[1:5, 6:10] <- fpkm_data[1:5, 6:10] + 4 combined_data <- data.frame(GeneID = gene_ids, fpkm_data) colnames(combined_data) <- c("GeneID", normal_samples, tumor_samples) # 写入文件 write.table(combined_data, file = file.path(desktop_path, "combined_RNAseq_FPKM.txt"), sep = "\t", quote = FALSE, row.names = FALSE) # ------------------------- 数据预处理 ------------------------- rt = read.table(file.path(desktop_path, "combined_RNAseq_FPKM.txt"), sep = "\t", header = T, check.names = F) rt = as.matrix(rt) rownames(rt) = rt[, 1] exp = rt[, 2:ncol(rt)] dimnames = list(rownames(exp), colnames(exp)) data = matrix(as.numeric(as.matrix(exp)), nrow = nrow(exp), dimnames = dimnames) data = avereps(data) data = data[rowMeans(data) > 0, ] data2 = as.data.frame(data) # ------------------------- 样本分组 ------------------------- exp_data_T = data2 %>% dplyr::select(str_which(colnames(.), "-01A")) exp_data_N = data2 %>% dplyr::select(str_which(colnames(.), "-11A")) rt = cbind(exp_data_N, exp_data_T) conNum = ncol(exp_data_N) treatNum = ncol(exp_data_T) # ------------------------- 数据标准化 ------------------------- pdf(file = file.path(desktop_path, "rawBox.pdf")) boxplot(rt, col = "blue", xaxt = "n", outline = F) dev.off() rt = normalizeBetweenArrays(rt) pdf(file = file.path(desktop_path, "normalBox.pdf")) boxplot(rt, col = "red", xaxt = "n", outline = F) dev.off() # ------------------------- 差异表达分析 ------------------------- group1 = sapply(strsplit(colnames(rt), "-"), "[", 4) group1 = sapply(strsplit(group1, ""), "[", 1) group1 = gsub("2", "1", group1) class <- c(rep("con", conNum), rep("treat", treatNum)) design <- model.matrix(~factor(class) + 0) colnames(design) <- c("con", "treat") fit <- lmFit(rt, design) fit <- contrasts.fit(fit, makeContrasts(con - treat, levels = design)) fit2 <- eBayes(fit) allDiff = topTable(fit2, adjust = 'fdr', n = Inf) write.table(allDiff, file = file.path(desktop_path, "limmaTab.xls"), sep = "\t", quote = F) # ------------------------- 筛选显著差异基因 ------------------------- diffLab = allDiff[with(allDiff, ((logFC > 1 | logFC < (-1)) & adj.P.Val < 0.05)), ] write.table(diffLab, file = file.path(desktop_path, "diffExp.xls"), sep = "\t", quote = F) diffExpLevel = rt[rownames(diffLab), ] write.table(diffExpLevel, file = file.path(desktop_path, "diffExpLevel.xls"), sep = "\t", quote = F) # ------------------------- 可视化单个基因 ------------------------- gene = "Gene1" data = t(rt[gene, , drop = F]) Type = c(rep(1, conNum), rep(2, treatNum)) exp = as.data.frame(cbind(data, Type)) colnames(exp) = c("gene", "Type") exp$Type = ifelse(exp$Type == 1, "Normal", "Tumor") exp$gene = log2(exp$gene + 1) exp$Type = factor(exp$Type, levels = c("Normal", "Tumor")) my_comparisons = list(c("Normal", "Tumor")) plot = ggboxplot(exp, x = "Type", y = "gene", color = "Type", xlab = "", ylab = paste0(gene, " expression (log2(FPKM+1))"), legend.title = "Type", palette = c("blue", "red"), add = "jitter") + stat_compare_means(comparisons = my_comparisons, label = "p.signif", symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns"))) pdf(file = file.path(desktop_path, paste0(gene, ".diff.pdf")), width = 5, height = 4.5) print(plot) dev.off() # ------------------------- 分析完成提示音 ------------------------- beepr::beep()
最新发布
07-09
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值