返回生信工具
序列分析

limma包分析TCGA数据库中某种癌症的差异基因表达

1.作用

  • 从TCGA中挖掘癌组织和正常组织的差异基因的表达，从而为我们探索特异基因寻找思路。
  • 从TCGA数据库下载好数据之后就可对里面的基因进行差异表达分析（我们使用的是limma包，因此最好使用FPKM数据）。

2.具体代码

这段代码是R语言中用于数据处理的，主要执行如下几个步骤： 使用read.table函数将名为"combined_RNAseq_FPKM.txt"的文件导入到R中，并将其保存为rt对象。该文件必须为制表符（"\t"）分隔的文本文件，并且包含表头。 将rt对象转化为矩阵格式。 使用rownames函数设置矩阵的行名，行名就是rt的第一列。 从第二列开始，全都赋值给新的exp矩阵。 创建一个新的矩阵data，并将exp矩阵的数字转化为数值格式。该矩阵的行名和列名与exp矩阵相同。 使用avereps函数处理data中的数据：对重复的行（即同一行名的多行数据）进行求平均值。 删除平均值为0的行。 将最后得到的数据矩阵data转换为数据框data2。 以上就是这段代码的基本逻辑，主要用于对RNA测序数据的预处理。

[!NOTE] 具体的实现代码

remove(list = ls()) ##清空当前环境 
setwd("E:/微生物-赵建组/课题进展/生信挖掘相关基因/03.TCGA胃癌数据/02.使用limma分析差异基因") ##💚💚💚💚💚设置路径

#  1.包的准备
#biocLite("limma")
#install.packages("gplots")
library(data.table)
library(tidyverse)
library(ggsignif) 
library(RColorBrewer)
library(limma)
library(ggplot2)
library(ggpubr)
library(beepr)
library(gplots)
library(pheatmap)



# 2.删除无表达基因，并对多个基因表达量取均值
###limma包分析TCGA数据要用FPKM，不然矫正后的P值会很大
###limma包分析TCGA数据要用FPKM，不然矫正后的P值会很大
###limma包分析TCGA数据要用FPKM，不然矫正后的P值会很大
###limma包分析TCGA数据要用FPKM，不然矫正后的P值会很大

rt=read.table("E:/微生物-赵建组/课题进展/生信挖掘相关基因/03.TCGA胃癌数据/02.使用limma分析差异基因/combined_RNAseq_FPKM.txt",sep="\t",header=T,check.names=F)    ##💚💚💚💚💚设置路径，要确保我们这个路径下存在这个文件，不然肯定运行不起来💚💚💚💚💚此外我们需要关注我们矩阵的间隔方式sep="\t"，就是我们是空格还是逗号❗❗❗❗❗❗❗这个之前我在这里报错下标出界，就是因为这里明明是空格，但是之前的源代码这里显示的是逗号
#多行取平均值
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]####❗❗❗❗❗❗❗这里原来是从第2列到之后的，但是我看了第2列不是基因表达量，它写的是蛋白基因表达量什么的
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)

#以01A和11A分组，正常放前面，肿瘤放后面
exp_data_T = data2%>% dplyr::select(str_which(colnames(.), "-01A")) # 匹配列名或用下示写法
nT = ncol(exp_data_T) 
exp_data_N = data2%>% dplyr::select(str_which(colnames(.), "-11A"))
nN = ncol(exp_data_N) 
rt= cbind(exp_data_N, exp_data_T)
#write.table(rt,file = "groupout.txt",sep="\t",quote=F)



#加载limma包，用于校正和比较差异
#rt=read.table("groupout.txt",sep="\t",header=T,check.names=F,row.names = 1)

#normalize以下为对数据的校正
pdf(file="rawBox.pdf")
#画箱线图
boxplot(rt,col = "blue",xaxt = "n",outline = F)
dev.off()
#校正
rt=normalizeBetweenArrays(rt)
pdf(file="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)
conNum=length(group1[group1==1])       #正常组样品数目
treatNum=length(group1[group1==0])     #肿瘤组样品数目
#differential差异分析
class <- c(rep("con",conNum),rep("treat",treatNum))  
design <- model.matrix(~factor(class)+0)
colnames(design) <- c("con","treat")
#算方差
df.fit <- lmFit(rt,design)
df.matrix<- makeContrasts(con - treat,levels=design)
fit<- contrasts.fit(df.fit,df.matrix)
#贝叶斯检验
fit2 <- eBayes(fit)
#输出基因
allDiff=topTable(fit2,adjust='fdr',n=Inf) 
#写入表格
write.table(allDiff,file="limmaTab.xls",sep="\t",quote=F)

#找出差异两倍以上，pvalue小于0.05，写入表格
diffLab <- allDiff[with(allDiff, ((logFC> 1 | logFC< (-1)) & adj.P.Val < 0.05 )), ]
write.table(diffLab,file="diffExp.xls",sep="\t",quote=F)

#差异基因表达水平，用于共表达
diffExpLevel <- rt[rownames(diffLab),]
write.table(diffExpLevel,file="diffExpLevel.xls",sep="\t",quote=F)

#可视化
gene="THBS2"   ##💚💚💚💚💚设置我们想要探究的基因
data=t(rt[gene,,drop=F])
Type=c(rep(1,conNum), rep(2,treatNum))
exp=cbind(data, Type)
exp=as.data.frame(exp)
colnames(exp)=c("gene", "Type")
exp$Type=ifelse(exp$Type==1, "Normal", "Tumor")
exp$gene=log2(exp$gene+1)
group=levels(factor(exp$Type))
exp$Type=factor(exp$Type, levels=group)
comp=combn(group,2)
my_comparisons=list()
for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
boxplot=ggboxplot(exp, x="Type", y="gene", color="Type",
                  xlab="",
                  ylab=paste0(gene, " expression"),
                  legend.title="Type",
                  palette = c("blue","red"),
                  add = "jitter")+ 
  stat_compare_means(comparisons=my_comparisons,symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif")
#输出图片
pdf(file=paste0(gene,".diff.pdf"), width=5, height=4.5)
print(boxplot)
dev.off()

3.最终效果

文件夹效果

[!NOTE]- 文件夹中多出来的图片，当然其中那个最初的FPKM文件是从一开始分析TCGA数据库中拿出来的 image.png

R语言中的效果

[!NOTE]- 最终出三张图片 image.png

4.遇到的问题

[!NOTE]- 原始代码，跑完dimnames和之后的矩阵全部都是NA

> library(pheatmap)
> rt=read.table("D:/生信代码复现/2.差异分析/2.差异分析/combined_RNAseq_FPKM.txt",sep="\t",header=T,check.names=F)    ##💚💚💚💚💚设置路径，要确保我们这个路径下存在这个文件，不然肯定运行不起来💚💚💚💚💚此外我们需要关注我们矩阵的间隔方式sep="\t"，就是我们是空格还是逗号
> #多行取平均值
> 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)
Warning message:
In matrix(as.numeric(as.matrix(exp)), nrow = nrow(exp), dimnames = dimnames) :
  强制改变过程中产生了NA

image.png

问的人工智能，这段代码注释

# 将rt数据框转换为矩阵格式
rt <- as.matrix(rt)

# 将rt矩阵的第一列设为行名
rownames(rt) <- rt[, 1]

# 从rt矩阵中提取第二列到最后一列的数据作为表达量数据 ###❗❗❗❗❗❗❗这个原来这里是以第2列到最后一列作为表达数据的，但是我看了一下就是我们最初获得的这个矩阵，第一列是基因symbol，但是第2列是基因的类型，就是蛋白编码基因什么的，所以我想从这里改成第3列及之后的作为基因表达量
exp <- rt[, 2:ncol(rt)]

# 创建包含行名和列名的列表，用于后面创建新矩阵时保持标签信息
dimnames <- list(rownames(exp), colnames(exp))

# 将表达量数据转换为数值，并且创建一个新的矩阵data
# nrow=nrow(exp)用于指定矩阵的行数
# dimnames=dimnames用于指定矩阵的维度名
data <- matrix(as.numeric(as.matrix(exp)), nrow = nrow(exp), dimnames = dimnames)

# 执行上述代码时出现的警告信息，意味着在数据转换过程中某些值无法转换为数值而变成了NA
# 这里没有包含代码处理这个问题的部分

解决办法 将这个“exp <- rt[, 2:ncol(rt)]”改成这个“exp <- rt[, 3:ncol(rt)]”