#文件需要：所有基因的logFC数据all.txt，所有基因的表达量数据allExp.csv，所有差异显著基因的logFC，diff_AllGene.txt，免疫基因的名字

#第一个小目标：提取差异表达的免疫基因的logFC数据及表达量数据

setwd("C:/Users/仙女/Documents/liver cancer")

rt=read.table("all.txt",sep="\t",header = T,check.names = F,row.names = 1)

allExp=read.csv("allExp.csv")

row.names(allExp)=allExp[,1]

allExp=allExp[,-1]

diffExp=read.table("diff_AllGene.txt",sep="\t",header = T,check.names = F,row.names = 1)

gene=read.table("immune_gene.txt",sep="\t",header = F)

immuneDiffALL=rt[intersect(gene[,1],row.names(rt)),]

#提取免疫基因相关的所有差异表达数据

immuneDiffGene=intersect(gene[,1],row.names(diffExp))

#从差异表达的基因中提取与免疫有关的基因,这里只提取了基因这一列

#用intersect函数将gene的第一列与diffExp的行名取交集

hmExp=allExp[immuneDiffGene,]

immuneDiffResult=immuneDiffALL[immuneDiffGene,]

#将差异免疫基因从免疫相关所有差异表达数据中提取出来

#输出免疫基因的差异结果

immuneDiffResult=rbind(ID=colnames(immuneDiffResult),immuneDiffResult)

write.csv(immuneDiffResult,file="immune_Diff.csv")

#输出差异免疫基因的表达量

immuneGeneExp=rbind(ID=colnames(hmExp),hmExp)

write.csv(immuneGeneExp,file="immuneGeneExp.csv")

write.csv(immuneGeneExp,file="TFGeneExp.csv")

第二个小目标：绘制火山图

library(pheatmap)

#绘制火山图

#横坐标是logFC，纵坐标是-log10（fdr）fdr是调整过得Pvalue，也可以是pvalue

pdf(file="vol.pdf",height = 5,width=5)

xMax=max(abs(as.numeric(as.vector(immuneDiffALL$logFC))))

yMax=max(-log10(immuneDiffALL$fdr))+1

a=as.numeric(as.vector(immuneDiffALL$logFC))

b=-log10(immuneDiffALL$fdr)

plot(a,b,xlab = "logFC",ylab = "-log10(fdr)",

    main="Volcano",ylim=c(0,yMax),xlim=c(-xMax,xMax),yaxs="i",pch=20,cex=0.8)

diffSub=subset(immuneDiffALL,fdr<fdrFilter & as.numeric(as.vector(logFC))>logFCfilter)

points(as.numeric(as.vector(diffSub$logFC)),-log10(diffSub$fdr),pch=20,col="red",cex=0.8)

diffSub=subset(immuneDiffALL,fdr<fdrFilter & as.numeric(as.vector(logFC))<(-logFCfilter))

points(as.numeric(as.vector(diffSub$logFC)),-log10(diffSub$fdr),pch=20,col="green",cex=0.8)

abline(v=0,lty=2,lwd=3)

dev.off()

第三个小目标：绘制热图

#绘制差异基因热图

hmExp=log2(hmExp+0.001)

#这个0.001是随便加的，hmExp是差异免疫基因的logFC等数据

Type=c(rep("N",conNum),rep("C",treaNum))

names(Type)=colnames(hmExp)

Type=as.data.frame(Type)

pdf(file="heatmap.pdf",height = 12,width = 15)

pheatmap(hmExp,

        annotation = Type,#分组

        color = colorRampPalette(c("green","black","red"))(50),

        cluster_cols = F,

        show_colnames = F,

        show_rownames = F,

        fontsize = 12,

        fontsize_row = 3,

        fontsize_col = 10)

dev.off()