代码（1）：热图

写在前面

大概就是找差异基因，这部分差异基因随着不同组别值表达上升，然后这部分差异基因还得是免疫相关基因。最后用热图展示出来

免疫基因

免疫基因可以在immport数据库中下载

通路可以在kegg上寻找

代码

# load data
options(stringsAsFactors = F)
rm(list = ls())
pwd = "E:\\process\\"
setwd(pwd)
library(readr)
library(dplyr)
data = read.csv("exp.txt",header = T,sep = "\t",
                fill=T,stringsAsFactors=F,skip = 1)

data = data[,c(1,7:15)]

# 1.transform probe 
library("clusterProfiler")
library(dplyr)
gene = data$Geneid# %>% as.data.frame()
df = data.frame(gene)
#BiocManager::install("AnnotationDbi")
library("AnnotationDbi")
library("org.Mm.eg.db")
df$symbol <- mapIds(org.Mm.eg.db,
                    keys=gene,
                    column="SYMBOL",
                    keytype="ENSEMBL",
                    multiVals="first")
df = na.omit(df)


# 2.subset  the exp 
library(limma)
## 2.1 subset the probe correspondant  to gene_symbols 
data = subset(data,Geneid%in%df$gene)
rownames(data) = data$Geneid 
df = df[match(df$gene,rownames(data)),]
data = data[,-1]
## 2.2 a gene correspondant different probe,you can choose mean or max 
data_finally = apply(data,2,function(x){
  tapply(x,df$symbol,max)
})

# 3.differentiall analyse 
library(edgeR)
dge <- DGEList(counts=data_finally)

group_list = c(rep("P",3),rep("y",3),rep("y",3))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))

## 3.1 filter some low express counts and transform to log2
keep <- filterByExpr(dge, design)
dge <- dge[keep,,keep.lib.sizes=FALSE]
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
rownames(design)=colnames(dge)

## 3.2 construct contrast matrix and search degs  
contrast_class=c("y-P","y-P","y-yin")

design
search_limma_degs=function(data_subset,group_list,contrast_class){
    design <- model.matrix(~0+factor(group_list))
    colnames(design)=levels(factor(group_list))
    head(design)
    exprSet=as.matrix(data_subset)
    mode(exprSet)="numeric"
    rownames(design)=colnames(exprSet)
    head(design)
    
    # 2.2 contrast.matrix
    
    contrast.matrix<-makeContrasts(contrasts = contrast_class,
                                   levels = design)
    
    # 3.search the degs 
    # 3.1 degs 
    ##step1,lmfit起过滤的手续
    fit <- lmFit(exprSet,design)
    ##step2
    fit2 <- contrasts.fit(fit, contrast.matrix)
    fit2 <- eBayes(fit2)  ## default no trend !!!
    ##eBayes() with trend=TRUE
    
    data_degs=list()
    for (i in 1:length(contrast_class)){
      print(i)
      tempOutput1 = topTable(fit2, coef=i, number =Inf)
      data_degs[[i]] = na.omit(tempOutput1) 
    }
    return(data_degs)
  }
data_degs=search_limma_degs(data_subset = logCPM,group_list,contrast_class)
deg_yangxing_PBS = data_degs[[1]]
deg_yinxing_PBS = data_degs[[2]]
deg_yangxing_yinxing = data_degs[[3]]

## 3.3 valcano plot 
nrDEG = deg_yangxing_PBS
valcano_pathway = function(nrDEG,i,contrast_class){
df=nrDEG
attach(nrDEG)
df$y= -log10(P.Value)
df$state=ifelse(df$adj.P.Val>0.05,'NO',
                ifelse( df$logFC >2,'UP',
                        ifelse( df$logFC < -2,'DOWN','NO') )
)
a = table(df$state)
df$state1<-factor(df$state,
                      levels = names(a),
                      labels = paste(names(a),as.numeric(a)))
df$name=rownames(df)
head(df)
library(ggplot2)
library(ggpubr)



## 4 kegg pathway 
#BiocManager::install("KEGGREST")
library("KEGGREST") 
listDatabases()  
gs<-keggGet('mmu04064')
genes<-unlist(lapply(gs[[1]]$GENE,function(x) strsplit(x,';'))) 
genelist_NFKB <- genes[1:length(genes)%%3 ==2] 

gs<-keggGet('mmu04660')
genes<-unlist(lapply(gs[[1]]$GENE,function(x) strsplit(x,';'))) 
genelist_Tcell <- genes[1:length(genes)%%3 ==2] 

gs<-keggGet('mmu04630')
genes<-unlist(lapply(gs[[1]]$GENE,function(x) strsplit(x,';'))) 
genelist_Jak <- genes[1:length(genes)%%3 ==2] %>% unique()

logCPM_mean = apply(logCPM,1,function(x){
  tapply(x, group_list, mean)
}) %>% t() %>% as.data.frame()
logCPM_mean = logCPM_mean[,c("PBS","yinxing","yangxing")]
logCPM_mean_order = subset(logCPM_mean,yangxing>yinxing&yinxing>PBS)

#5.免疫相关基因
gene_list = read.table("GeneList.txt",header = T,sep = "\t",
                       fill=T,stringsAsFactors=F)
library(Hmisc) 
immune_gene = gene_list$Symbol %>% tolower() %>% capitalize() %>% unique()



# 6.都差异的基因
data_degs_1 = sapply(data_degs, function(x){
  rownames(subset(x,adj.P.Val<0.05))
})
con_gene = intersect(data_degs_1[[2]],data_degs_1[[3]]) %>% unique()


# 取特定通路顺序
kegg_order_jak = logCPM_mean_order[genelist_Jak,] %>% .[con_gene,] %>% .[immune_gene,] %>% 
  na.omit() %>% unique()

kegg_order_t = logCPM_mean_order[genelist_Tcell,] %>% .[con_gene,] %>% .[immune_gene,] %>% 
  na.omit() %>% unique()

kegg_order_NFKB = logCPM_mean_order[genelist_NFKB,] %>% .[con_gene,] %>% .[immune_gene,] %>% 
  na.omit() %>% unique()

# 画热图
#BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)
library(circlize)
data_plot = rbind(kegg_order_jak,kegg_order_NFKB,kegg_order_t)
ha2 <- rowAnnotation(df = data.frame(
  Pathway = c(rep("Jak/Stat signaling pathway",nrow(kegg_order_jak)),
              rep("NF-kappa B signaling pathway",nrow(kegg_order_NFKB)),
              rep("T cell receptor signaling pathway",nrow(kegg_order_t))
              )
),
col = list(Pathway = c("Jak/Stat signaling pathway" = "blue",
                    "NF-kappa B signaling pathway" = "tomato",
                    "T cell receptor signaling pathway" = "red"
                    )
           )
)
scaled_expr = scale(t(data_plot)) %>% t()
pdf(file="complexheatmap.pdf")
Heatmap(
  scaled_expr,
  name="Color_key",        # 图例的名字
  #top_annotation = ha1,    # 将注释放在最上面
  left_annotation =ha2 ,   # 将注释放在最左边
  cluster_rows = FALSE,  
  cluster_columns  = FALSE,# 不对行进行聚类
  #show_row_names = FALSE,  # 不显示行名
  column_names_gp = gpar(fontsize = 12),  # 列名的字体大小
  col=colorRamp2(
    c(min(scaled_expr), 0, max(scaled_expr)), 
    c("blue", "white", "red")
  )
)
dev.off()

png(file="complexheatmap.png",width = 520)
Heatmap(
  scaled_expr,
  name="Color_key",        # 图例的名字
  #top_annotation = ha1,    # 将注释放在最上面
  left_annotation =ha2 ,   # 将注释放在最左边
  cluster_rows = FALSE,  
  cluster_columns  = FALSE,# 不对行进行聚类
  #show_row_names = FALSE,  # 不显示行名
  column_names_gp = gpar(fontsize = 12),  # 列名的字体大小
  col=colorRamp2(
    c(min(scaled_expr), 0, max(scaled_expr)), 
    c("blue", "white", "red")
  )
)
dev.off()

结果