无参转录组差异表达分析（重点火山图）

拿到表达矩阵之后，首先构建dds。

library(limma)
library(ggplot2)
library(DESeq2)
library(pheatmap)
options(stringsAsFactors = F)
raw_count <- read.table('matrix.counts.matrix', header = T,sep = '\t')
a<- raw_count[,-1]    #删掉第一列
count_data <- a[,-4:-6]  #我删掉了第4到6列，根据自己的数据来
exprSet=apply(count_data,2, as.integer)   #把所有数据转化为整数，因为他只认整数
rownames(exprSet)=raw_count[,1]
condition <- factor(c("trt","trt","trt","untrt","untrt","untrt"))  
col_data <- data.frame(row.names = colnames(exprSet), condition)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = col_data,design = ~ condition)
dds_filter <- dds[ rowSums(counts(dds))>1, ]   #去掉表达量小于1的基因
dds_out <- DESeq(dds_filter)   #计算差异表达的基因数目
res <- results(dds_out)
summary(res)  #查看结果

差异表达结果

可以看到有2053的上调基因，有2129的下调基因以及140624的低表达的基因。这三个数据就是我们后续画火山图需要的前景基因和背景基因。

以上代码画其他图都会用到，主要是构建dds，挑选差异表达基因，供后续分析。
一：画几个图直观的看一下差异表达情况

rld <- rlogTransformation(dds_out)
exprSet_new=assay(rld)
par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(2,2))
boxplot(exprSet, col = cols,main="exprSet_expression value",las=2)
boxplot(exprSet_new, col = cols,main="exprSet_new_expression value",las=2)
hist(exprSet)
hist(exprSet_new)

箱线图

二：MAplot

plotMA(res, ylim = c(-10,10))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
  
})

MAplot

三：pheatmap
1.0

table(res$padj<0.05)    ##看一下p小于0.05的有多少
res_deseq <- res[order(res$padj),]
res_deseq=as.data.frame(res_deseq)
nrDEG=res_deseq
nrDEG=na.omit(res_deseq)   ##去掉NA
choose_gene=head(row.names(nrDEG),50)  ##画top50的gene
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix,cluster_row = T)

pheatmap

2.0

rld <- rlogTransformation(dds_out,blind = F)
write.csv(assay(rld),file="JSM.DESeq2.pseudo.counts.csv")
topVarGene <- head(order(rowVars(assay(rld)),decreasing = TRUE),50)
mat  <- assay(rld)[ topVarGene, ]
anno <- as.data.frame(colData(rld)[,c("condition","sizeFactor")])
pheatmap(mat, annotation_col = anno, cluster_cols = T)

pheatmap

四：volcano-ggplot

resLFC <- lfcShrink(dds_out,coef = 2,res=res)
write.csv(resLFC, file ="mm_deseq2_resLFS.csv" )

到这一步以及生产了我们花火山图所需要的差异表达基因文件（mm_deseq2_resLFS.csv）,查看一下这个文件：

csv文件

---

可以看到这个csv文件是包含7列，但是我们画画火山图需要上调下调和既不上调也不下调的基因，所以我们就需要这个文件中的1(gene_id),3(log2FoldChange),7(padj)列，并且要把7列中的NA删掉,至于这三列代表的什么意思，自己去查。
到Linux中运行如下代码：

##########################切记##########################
##########################切记##########################
你要看你的csv文件的分割符是什么，我也是看了好久才明白下面这串代码:

1. 要是你的csv文件时以","分割，这下面这串代码的-F 要用‘，’.我的就是以‘，’分割的。
2. 要是你的csv文件是以‘\t’分割的话，切记-F要用‘\t’.不然这串代码只会删掉表头，数据不会删掉的哦，因为他没法识别你是以什么分割的。

perl -F',' -alne '$F[5]=~s/\r//;if(/baseMean/){print "gene\tlog2FoldChange\tpadj\tsignificant"} else{unless(/NA/){if($F[2] >=1 && $F[5]<0.05){print "$F[0]\t$F[2]\t$F[5]\tup"} elsif($F[2]<=-1 && $F[5]<0.05){print "$F[0]\t$F[2]\t$F[5]\tdown"} else{print "$F[0]\t$F[2]\t$F[5]\tno"}}}' mm_deseq2_resLFS.csv >volcano.txt

最后就会生成如下的包含4列的文件，这个文件就是我们画火山图的东西。

火山图所需要的差异表达基因

再到R中跑如下代码：

data <-read.table(file="volcano.txt",header = TRUE,sep = '\t')
volcano <-ggplot(data = data,aes(x=log2FoldChange,y= -1*log10(padj)))+geom_point(aes(color=significant))+scale_color_manual(values = c("red","grey","blue")) + labs(title="Volcano_Plot",x=expression((log[2](FC)), y=expression(-log[10](padj)) ))+geom_hline(yintercept=1.3,linetype=4)+geom_vline(xintercept=c(-1,1),linetype=4)
volcano

最后看一下生成的图片：

volcano-ggplot

是不是还阔以？？