• 火山图
• 热图
• 热图
• GO和KEGG分析

火山图

library(tidyverse)
library(ggplot2)
library(ggrepel)
head(res)
DEG_all <- res %>% rownames_to_column("symbol") %>% 
  dplyr::select(symbol,logFC=log2FoldChange,Pvalue=padj)%>%
  mutate(direction = factor(ifelse(Pvalue < 0.01 & abs(logFC)>2,
                                  ifelse(logFC>1,"Up","Down"),
                                  "NS"),levels = c("Up","Down","NS"))) %>%
  na.omit()

ggplot(DEG_all,aes(x=logFC,y=-log10(Pvalue),colour=direction))+
  geom_point(alpha=0.6)  +
  scale_color_manual(values = c("#DC143C","#00008B","#808080"))+
  geom_text_repel(data=DEG_all %>%filter(Pvalue<0.01,abs(logFC)>9),
                  aes(label=symbol),
                  size=3,
                  segment.color="black",
                  show.legend = F)+
  xlab(expression(log[2]("Fold Change")))+
  ylab(expression(-log[10]("Adjusted P Value")))+
  ggtitle("Volcano Plot")+
  theme_bw()+
  theme(legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  geom_vline(xintercept = c(-2,2),lty=2,col="black",lwd=0.6)+
  geom_hline(yintercept = -log10(0.01),lty=2,col="black",lwd=0.6)

DEG_all %>%filter(symbol %in% crc_mrna_deg$gene_id)

ggplot(DEG_all,aes(x=logFC,y=-log10(Pvalue),colour=direction))+
  geom_point(alpha=0.6)  +
  scale_color_manual(values = c("#DC143C","#00008B","#808080"))+
  geom_text_repel(data=DEG_all %>%filter(symbol %in% ppi_gene$gene_id),
                  aes(label=symbol),
                  size=3,
                  segment.color="black",
                  show.legend = F)+
  xlab(expression(log[2]("Fold Change")))+
  ylab(expression(-log[10]("Adjusted P Value")))+
  ggtitle("Volcano Plot")+
  theme_bw()+
  theme(legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  geom_vline(xintercept = c(-2,2),lty=2,col="black",lwd=0.6)+
  geom_hline(yintercept = -log10(0.01),lty=2,col="black",lwd=0.6)

tiff(file="figures/crc_mRNA_Volcano.tiff",
     width = 35,height = 22,units = "cm",compression = "lzw",bg = "white",res=300)
ggplot(DEG_all,aes(x=logFC,y=-log10(Pvalue),colour=direction))+
  geom_point(alpha=0.6)  +
  scale_color_manual(values = c("red","green","#000000"))+
  xlab(expression(log[2]("Fold Change")))+
  ylab(expression(-log[10]("Adjusted P Value")))+
  ggtitle("mRNA Volcano")+
  theme_bw()+
  theme(legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5,size = 20, face= "bold"),
        axis.title =  element_text(size=18),
        legend.text= element_text(size=15),
        aspect.ratio=1)+
  geom_vline(xintercept = c(-2,2),lty=2,col="black",lwd=0.6)+
  geom_hline(yintercept = -log10(0.01),lty=2,col="black",lwd=0.6) +
  ggsave("figures/mRNA_Volcano_square.png")

dev.off()

热图

• 筛选标准差最大的基因绘制热图

 cg=names(tail(sort(apply(dat,1,sd)),100)) ##取表达量标准差最大的100行的行名

pheatmap(gset_expr_annot_20)
gset_expr_annot_20 <-  t(scale(t(gset_expr_annot_20)))
col <- colorRampPalette(c("green","#000000","red"))(256)
pheatmap(gset_expr_annot_20,col=col,fontsize=8,scale="row",annotation_col = gset_pd_f,filename = "figures/pheatmap_group.png")
pheatmap(gset_expr_annot_20,col=col,fontsize=8,scale="row",filename = "figures/pheatmap.png")

library(pheatmap)

vsd2 <- assay(vsd) 
genes_sig <- DEG_all %>%
  filter(Pvalue < 0.01 & abs(logFC)>2) %>%
  arrange(Pvalue) %>%
  dplyr::slice(1:500)
table(genes_sig$direction)
DEG_all %>% filter(symbol=="")

heatmap_df <- vsd2[genes_sig$symbol,]
write.csv(heatmap_df,file = "result/heatmap_df.csv")

#c("red", "white", "blue")
col <- colorRampPalette(c("red","#000000","green"))(256)
hmMat <- log10(heatmap_df+0.001)
pheatmap(hmMat,col=col,fontsize=3,main="hdd")
?pheatmap
#,filename = "figures/mRNA_pheatmap.png"
heatmap_df2 <-  crc_mRNA_count_order[genes_sig$symbol,]

crc_mRNA_count_order <- cbind(crc_mRNA_count[,rownames(crc_mRNA_sample[crc_mRNA_sample$group=="Norm",])],crc_mRNA_count[,rownames(crc_mRNA_sample[crc_mRNA_sample$group=="Tumm",])])
dim(crc_mRNA_sample)
dim(crc_mRNA_count)
table(crc_mRNA_sample$Group)
identical(rownames(crc_mRNA_sample),colnames(crc_mRNA_count_order))
DEG_df <- 



dds <- DESeqDataSetFromMatrix(countData = crc_mRNA_count_order,colData = crc_mRNA_sample,design = ~ Group)
## PCA
vsd <- vst(dds,blind=T) # 数据均一化
vsd_df <- assay(vsd)
head(vsd_df)
DESeq2::plotPCA(vsd,intgroup="Group")
keep <- rowSums(counts(dds)>0) > 44
dds_filt <- dds[keep,]
dim(dds_filt)
dim(dds)
dds2 <- DESeq(dds_filt)
sizeFactors(dds2)
resultsNames(dds2)
res <- results(dds2)
res <- as.data.frame(res)

热图

library(gplots)
heatmap.2(hmMat,col=bluered(100),trace = 'none', density.info = "none")
?heatmap.2

heatmap_df2 <-  crc_mRNA_count_order[genes_sig$symbol,]

write.csv(heatmap_df2,file = "result/heatmap_df2.csv")

hmMat <- scale(heatmap_df2)
log2(hmMat)
pheatmap(hmMat,col=col,scale = "row",fontsize=2.5,)
heatmap.2(df)
library(ComplexHeatmap)

pheatmap(dat)
n=t(scale(t(dat)))
n[n>2]=2 #限定上限，使表达量大于2的等于2
n[n< -2]= -2 #限定下限，使表达量小于-2的等于-2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=cluster)
rownames(ac)=colnames(n)
pheatmap(n,annotation_col = ac,
show_colnames =F,show_rownames = T)
n[n< -1]= -1 # 重新限定下限，使表达量小于-2的等于-2
n[1:4,1:4]
pheatmap(n,annotation_col = ac,
show_colnames =F,show_rownames = T)

GO和KEGG分析

## 网络图
cnetplot(egobp,circular=T,colorEdge=T)
## 热图
heatplot(egobp)
## 富集图
x2 <- pairwise_termsim(egobp)
emapplot(x2,layout = "kk",showCategory=20)
## Venn图
library(UpSetR)
upsetplot(egobp)

png(file="figures/kegg_cnetplot.png",
    width = 30,height = 22,units = "cm",bg = "white",res=300)
cnetplot(kk,circular=T,colorEdge=T)
dev.off()
## 热图
heatplot(kk)
## 富集图
x2 <- pairwise_termsim(kk)
emapplot(x2,layout = "kk",showCategory=20)
## Venn图
library(UpSetR)
upsetplot(kk)


dotplot(kk,showCategory = 15)
barplot(kk,showCategory = 10)