1.导入数据
利用htseq计数得到的表达矩阵，下载到Windows中（再拷贝一个 .Rproj到下面文件夹中，直接双击R Project进入到当前路径下的R中进行后续操作）
做R分析的文件夹
rm(list = ls())
options(stringsAsFactors = T)
library(DESeq2)
#载入数据，只有1个对照，2个实验组（重复）
control <- read.table("SRR3589956_matrix.count",
sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("SRR3589957_matrix.count",
sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("SRR3589958_matrix.count",
sep="\t",col.names = c("gene_id","rep2"))
#将3个数据根据gene_id合并到一起
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count
raw_count
#删除raw_count的前5行
raw_count_filt <- raw_count[-1:-5,]
#gsub匹配raw_count_filt的gene_id行，生成ENSEMBL（ENSEMBL是一个字符串格式的向量）
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
#将raw_count_filt的行名替换为ENSEMBL（猜想应该是想要去掉raw_count_filt行名（即其gene_id）的小数点后面的数字）
row.names(raw_count_filt) <- ENSEMBL
#解决control没有重复的问题，自己创造一个control2#注意：创造数据不是无根据捏造。要考虑到高通量测序的read是默认符合泊松分布的，于是可以这样做：（1）计算KD重复组的均值差，作为泊松分布的均值；（2）使用概率函数rpois()随机产生一个数值，前一步的均值作为lambda；（3）对一些read count 低于均值的直接加上对应KD重复组之间的差值delta_mean <- abs(mean(raw_count_filt$rep1) - mean(raw_count_filt$rep2))sampleNum <- length(raw_count_filt$control)sampleMean <- mean(raw_count_filt$control)control2 <- integer(sampleNum)for (i in 1:sampleNum){ if(raw_count_filt$control[i] < sampleMean){ control2[i] <- raw_count_filt$control[i] + abs(raw_count_filt$rep1[i] - raw_count_filt$rep2[i]) } else{ control2[i] <- raw_count_filt$control[i] + rpois(1,delta_mean) }}#将创造好的control2添加到矩阵中去raw_count_filt$control2 <- control2
2.使用DESeq2进行差异基因分析
DESeq2要求的数据是raw count， 没必要进行FPKM/TPM/RPFKM/TMM标准化
#构建DESeq2所需的DESeqDataSet对象library(DESeq2)countData <- raw_count_filt[,2:5] #countData为基因表达矩阵condition <- factor(c("control","KD","KD","control")) #condition为处理条件dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
#过滤low_count的数据
nrow(dds)
dds <- dds[ rowSums(counts(dds)) > 1, ]
nrow(dds) #看下图，过滤了将近一半的数据
过滤前后基因的数量
#使用DESeq进行差异表达分析
dds <- DESeq(dds) #出现如下红色字符，表示运行成功
出现图上红色字段，表示运行成功
#查看结果
res <- results(dds)
mcols(res, use.names = TRUE)
summary(res)
res每一列表达的含义
基因差异表达的结果
可以看到上调和下调的基因占比，low_count的比率比较高，所以大部分基因表达量不高，主要集中在0附近（log2(1)=0,也就是变化1倍）。
3.画MA图
M代表： log fold change，能衡量基因表达量变化，上调或下调。
A代表：每个基因的count的均值。
#做一个没有经过 statistical moderation平缓log2 fold changes的MA图
library(ggplot2)
#画个MA图
plotMA(res, ylim = c(-5,5))
#标注出MA图中p值最小的基因
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")
})
MA图，蓝色圆圈标注的为p值最小的基因
#如果经过lfcShrink 收缩log2 fold change， 结果会好看很多
res.shrink <- lfcShrink(dds, contrast = c("condition","KD","control"), res=res)
plotMA(res.shrink, ylim = c(-5,5))
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")
})
MA图，蓝色圆圈标注的为p值最小的基因
4.提取差异基因分析的结果
##提取差异基因分析的结果
table(res$padj<0.05) #计算p值小于0.05的基因的个数
res_deseq <- res[order(res$padj),] #根据res的padj值进行排序，并赋值给res_deseq
diff_gene_deseq2 <- subset(res_deseq, padj<0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) #将res_deseq过滤：要求是p值小于0.05且log2FoldChange在1~-1范围内的
diff_gene_deseq2 <- row.names(diff_gene_deseq2)
res_diff_data <- merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
write.csv(res_diff_data,file = "11_30_humen_data.csv",row.names = F)
##得到的11_30_humen_data.csv为差异表达分析的结果。
5.火山图
#需要先将差异表达分析得到的结果11_30_humen_data.csv在Linux中进行处理
#第三列log2foldchange大于1为上调，小于-1下调，其他的不显著（需要构建significant列：up,down,no）
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"}}}' 11_30_humen_data.csv >volcano.txt
#火山图
library(ggplot2)
data <-read.table(file="volcano.txt",header = TRUE, row.names =1,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
火山图
6.gene聚类热图
library(genefilter)
library(pheatmap)
rld <- rlogTransformation(dds,blind = F)
#生成一个mm.DESeq2.pseudo.counts.csv矩阵
write.csv(assay(rld),file="mm.DESeq2.pseudo.counts.csv")
#head了前20个差异较大的基因
topVarGene <- head(order(rowVars(assay(rld)),decreasing = TRUE),20)
mat <- assay(rld)[ topVarGene, ]
# mat <- mat - rowMeans(mat) 减去一个平均值，让数值更加集中。得到的是第二个图。
anno <- as.data.frame(colData(rld)[,c("condition","sizeFactor")])
pheatmap(mat, annotation_col = anno)
聚类图之一
聚类图之二
做了2个聚类的图，都只做了前20个基因的聚类。第二幅图中的基因更加的集中。
