生信分析做制作火山图代码

#火山图

plot(tT$logFC,-log10(tT$P.Value),pch=16,main="Volcano map",cex.lab=1.5,

    xlab="log2FoldChange",ylab="-log10Pvalue")

points(tT[tT$logFC>1&tT$P.Value<0.05,]$logFC,pch=16,

      -log10(tT[tT$logFC>1&tT$P.Value<0.05,]$P.Value),

      col='red')

points(tT[tT$logFC< -1&tT$P.Value<0.05,]$logFC,pch=16,

      -log10(tT[tT$logFC< -1&tT$P.Value<0.05,]$P.Value),

      col='green')

abline(h=-log10(0.05),v=c(-1,1))

#signicant difference genes

interesting_genes=head(tT[order(tT$P.Value,decreasing = F),],10)

text(interesting_genes$logFC,-log10(interesting_genes$P.Value),

    interesting_genes$Gene.symbol,cex=0.8,pos=1)

library(pheatmap)

#筛选差异显著的基因

selected_gens=tT[tT$P.Value<0.05&abs(tT$logFC)>1,]

dim(selected_gens)

dim(tT)

##基因表达矩阵

head(ex)[,1:5]

##筛选差异表达的基因

selected_significant_genes=ex[rownames(ex) %in% rownames(selected_gens),]

dim(selected_significant_genes)

#分组信息

gsms <- paste0("11111111111111111111111111111111111111111111111111",

              "11111111111111111111000000000000000000000000000000",

              "0000000")

sml <- c()

for (i in 1:nchar(gsms)) { sml[i] <- substr(gsms,i,i) }

group=ifelse(sml=="G1","cancer","normal")

#70个cancer，37个normal

table(group)

#定义样本属性

annotation_col=data.frame(Sample=group)

rownames(annotation_col)=colnames(selected_significant_genes)

colnames(selected_significant_genes)=rownames(annotation_col)

pheatmap(selected_significant_genes,annotation_col = annotation_col,scale = 'row',

        show_rownames = F,color = colorRampPalette(c('green','black','red'))(50),

        clustering_distance_cols = 'maximum',treeheight_col=15,

        cluster_rows = T,treeheight_row = F,show_colnames = F,cellwidth = 5,

        annotation_colors = list(Sample=c(cancer='red',normal='green')))