文件要求：正常与肿瘤组的mRNA表达量
本次前50列是正常组，后374是肿瘤样本
setwd("C:/Users/仙女/Documents/liver cancer")
library(limma)
inputFile="gene_mRNA_expression.txt"
fdrFilter=0.05#fdr临界值
logFCfilter=1#logFC临界值
conNum=50#normal组织样品数目
treaNum=374#cancer组织样品数目
#读取输入文件
outTab=data.frame()
grade=c(rep(1,conNum),rep(2,treaNum))
#将正常组织变为1，肿瘤组织为2
rt=read.table(inputFile,sep="\t",header = T,check.names = F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]#去除第一列基因名
head(exp)[,1:10]
head(rt)[,1:10]
#这里是要制作表达基因矩阵
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.2,]
#筛选数据,把低表达的基因去掉
#差异分析,Wilcoxon检验方法
for (i in row.names(data)) {
geneName=unlist(strsplit(i,"\\|",))[1]
geneName=gsub("\\/","_",geneName)
rt=rbind(expression=data[i,],grade=grade)
rt=as.matrix(t(rt))
wilcoxTest<-wilcox.test(expression~grade,data=rt)
conGeneMean=mean(data[i,1:conNum])
treatGeneMean=mean(data[i,(conNum+1):ncol(data)])
logFC=log2(treatGeneMean)-log2(conGeneMean)
pValue=wilcoxTest$p.value
conMed=median(data[i,1:conNum])
treatMed=median(data[i,(conNum+1):ncol(data)])
diffMed=treatMed-conMed
if(((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){
outTab=rbind(outTab,cbind(gene=i,conMean=conGeneMean,treaMean=treatGeneMean,logFC=logFC,pValue=pValue
)) }
}
pValue=outTab[,"pValue"]
fdr=p.adjust(as.numeric(as.vector(pValue)),method = "fdr")
outTab=cbind(outTab,fdr=fdr)
#输出所有基因的差异情况
write.table(outTab,file="all.txt",sep="\t",row.names = F,quote = F)
#输出经过平均值处理之后的所有基因表达情况
write.csv(data,file="allExp.csv")
#输出差异表格
差异表格有点出错，还是自己用excel解决吧
