#文件需要：所有基因的logFC數(shù)據(jù)all.txt，所有基因的表達(dá)量數(shù)據(jù)allExp.csv，所有差異顯著基因的logFC，diff_AllGene.txt，免疫基因的名字

#第一個(gè)小目標(biāo)：提取差異表達(dá)的免疫基因的logFC數(shù)據(jù)及表達(dá)量數(shù)據(jù)

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)),]

#提取免疫基因相關(guān)的所有差異表達(dá)數(shù)據(jù)

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

#從差異表達(dá)的基因中提取與免疫有關(guān)的基因,這里只提取了基因這一列

#用intersect函數(shù)將gene的第一列與diffExp的行名取交集

hmExp=allExp[immuneDiffGene,]

immuneDiffResult=immuneDiffALL[immuneDiffGene,]

#將差異免疫基因從免疫相關(guān)所有差異表達(dá)數(shù)據(jù)中提取出來

#輸出免疫基因的差異結(jié)果

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

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

#輸出差異免疫基因的表達(dá)量

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

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

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

第二個(gè)小目標(biāo)：繪制火山圖

library(pheatmap)

#繪制火山圖

#橫坐標(biāo)是logFC，縱坐標(biāo)是-log10（fdr）fdr是調(diào)整過得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()

第三個(gè)小目標(biāo)：繪制熱圖

#繪制差異基因熱圖

hmExp=log2(hmExp+0.001)

#這個(gè)0.001是隨便加的，hmExp是差異免疫基因的logFC等數(shù)據(jù)

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