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)
## 网络图 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)