TCGA_RNA-seq_limma分析

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')