转录组分析RNA-Seq（续）

if (!requireNamespace("BiocManager"))
  install.packages("BiocManager")
BiocManager::install("DEseq2")

# 加载DESseq2 R包
library(DESeq2)
# 读取基因表达矩阵数据
countData <- read.table("C:\\Users\\Jeff\\OneDrive\\桌面\\online_count.txt",row.names=1,header=F))
'''
BTW 如果你不想输入太长的绝对路径，你可以先设置工作路径，就像python中os模块中os.chdir()一样
setwd("C:\\Users\\Jeff\\OneDrive\\桌面")
countData <- read.table("online_count.txt",row.names=1,header=F))
'''
countData <- countData[,c(4:11)]
# 查看数据的头部
head(countData)
# 删除表达量均值小于1的基因
countData <- countData[rowMeans(countData)>1,]
# 构造一个处理组-列交联矩阵
condition <- factor(c(rep("control",4),rep("treatment",4)))
colData <- data.frame(row.names=colnames(countData), condition=condition)
# 构造DESeq2计算矩阵
dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~condition)
# 差异分析
dds <- DESeq(dds)
# 提取差异分析结果
res <- results(dds, contrast=c("condition","treatment","control"))
# 按照矫正后的p值进行结果的排序
res <- res[order(res$padj),]
# 提取所有的差异基因，其中阈值设置为p-adjust value < 0.01，log2FoldChange的绝对值大于1
# 另外如若想单独提取上下调基因的相关信息，修改subset的筛选条件即可：
#p-adjust value < 0.01，log2FoldChange > 1
#p-adjust value < 0.01，log2FoldChange <- 1
diff_gene_deseq2 <-subset(res,padj < 0.01 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene_deseq2 <- row.names(diff_gene_deseq2)
# 整合差异结果以及各文库标准化后的定量数值至同一数据框
resdata <-  merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
#得到csv格式的差异表达分析结果
write.csv(resdata,file= "DESEQ_result.csv",row.names = F)

#导入R包，如果未安装请按DESeq2的方法安装
library(RColorBrewer)
library(ggpubr)
library(ggthemes)
#在这里我们定义一个函数对接DESeq2产出的结果
#函数需要传入res和pdf_path两个参数
VolcanoPlot=function(res,pdf_path)
{
  res$logP <- -log10(res$padj) #padj取-log值生成新列
  res <- as.data.frame(res) #转为数据框形式
  res$group <- 'not-significant' #根据选定的阈值对基因添加label（上调？下调？还是没差异？）
  res$group[which((res$padj <=0.01) & (res$log2FoldChange >=1) )] = 'up-regulated'
  res$group[which((res$padj <=0.01) & (res$log2FoldChange <=-1) )] = 'down-regulated'
  table(res$group) #简单看下多少个上调基因多少个下调基因
  ggscatter(res, #画图咯
            x = 'log2FoldChange', #x轴数据
            y ='logP' , #y轴数据
            color = 'group', #用标签区分散点颜色为三个颜色
            size = 1, #散点大小
            palette = c("#2f5688","#BBBBBB","#CC0000"), #手动定义三个颜色应为什么色
            xlab = "log2FoldChange", #x轴标签
            ylab = "-log10(padj)" #y轴标签
           )+theme(text = element_text(size = 10),legend.position = "right")+ #theme中定义全局字体大小和图例的位置（在右边）
  geom_hline(yintercept = 2,linetype='dashed')+  #绘制y轴阈值线 即padj =0.01处
  geom_vline(xintercept = c(-1,1),linetype='dashed') #绘制x轴阈值线 即log2FoldChange分别为1和-1处
  ggsave(pdf_path, width = 6.4, height = 3.6 ,dpi = 300) #保存图片，指定pdf高度和宽短，dpi分辨率为300
  return(TRUE)
}
#在这里 我们运行一下我们定义的函数
VolcanoPlot(res,"tmp.pdf")

## PCA
library("FactoMineR")
library("factoextra")
count.pca <- counts(dds) %>% as.matrix %>% t %>% PCA(.,graph = F)
pdf("PCA.pdf")
fviz_pca_ind(count.pca, 
               col.ind=condition, 
               mean.point=T, 
               addEllipses = F, 
               legend.title="Groups")+ 
theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))
dev.off()
## Heatmap
library(pheatmap)
pdf("Heatmap.pdf")
normalized_counts <- counts(dds, normalized=TRUE)
normalized_counts_mad <- apply(normalized_counts, 1, mad)
normalized_counts <- normalized_counts[order(normalized_counts_mad, decreasing=T), ]
Heatmap_Data=normalized_counts[rownames(diff_gene_deseq2),]
annotation_col = data.frame(
    Conditions  = factor(c(rep("control", 4),rep("treatment", 4)) ))
colnames(Heatmap_Data)=factor(colnames(Heatmap_Data))
pheatmap(Heatmap_Data,
           scale = "row",
           show_rownames = F,
           cluster_cols = F,
           main = "DAR Heatmap",
           cutree_rows = 2,
           gaps_col  = 4,
           color = colorRampPalette(c("navy", "white", "firebrick3"))(50),
           ColSideColors = rep(c("pink", "purple"), each = 4))
dev.off()
### pearson
library(amap)
library(gplots)
rld <- rlog(dds, blind=TRUE)
rlogMat <- assay(rld)
rlogMat <- rlogMat[order(normalized_counts_mad, decreasing=T), ]
hmcol <- colorRampPalette(brewer.pal(9, "GnBu"))(100)
pearson_cor <- as.matrix(cor(rlogMat, method="pearson"))
hc <- hcluster(t(rlogMat), method="pearson")
heatmap.2(pearson_cor, Colv = as.dendrogram(hc), Rowv=as.dendrogram(hc), symm=T, trace="none",
          col=hmcol, main="The pearson correlation of each
sample")
dev.off()