RNA-seq，R语言部分

代码如下：谢谢指导

#1-相关性分析
library('corrplot')
exprSet <- read.table('all_id.txt',header = T,sep = '\t',fill = T)
exprSet <- exprSet[,c(1,7,8,10,11)]
corrplot(data_cor,type="upper")
#2-相关性热图
require(pheatmap)
pheatmap(scale(cor(log2(exprSet+1))))
png('heatmap.png')
pheatmap(scale(cor(log2(exprSet+1))))
dev.off()
#3-差异性分析
#3-1处理一下读入的all-id。countd
library(stringr)
exprSet$Geneid <- substr(exprSet$Geneid,1,15)  #把基因ID后面的小数点去掉
exprSet=exprSet[!duplicated(exprSet$Geneid),]  #去掉重复的基因ID。不知道是直接去掉，还是两个相同的ID取平均值？
row.names(exprSet) <- exprSet$Geneid
exprSet <- exprSet[,-1]
#3-2构建样本矩阵
condition <- factor(c("control","control","treat","treat"))
colData <- data.frame(row.names=colnames(exprSet), condition)
#3-3 DESeq2差异性分析
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ condition)
dds <- DESeq(dds)        
res <- results(dds, contrast=c("condition","control","treat"))
DEG <- as.data.frame(res)
DEG = DEG[order(DEG$pvalue),]
write.csv(DEG2,file="DEseq差异分析结果.csv")
#3-4 DESeq2结果可视化
#定义筛选
DEG[which(DEG$padj %in% NA),'sig'] <- 'no diff'
DEG[which(DEG$log2FoldChange >= 1 & DEG$padj < 0.05),'sig'] <- 'rich (p.adj < 0.05, log2FC >= 1)'
DEG[which(DEG$log2FoldChange <= -1 & DEG$padj < 0.05),'sig'] <- 'down (p.adj < 0.05, log2FC <= -1)'
DEG[which(abs(DEG$log2FoldChange) < 1 | DEG$padj >= 0.05),'sig'] <- 'no diff'
#画图
volcano_p <- ggplot(DEG, aes(log2FoldChange, -log(padj, 10))) +
  geom_point(aes(color = sig), alpha = 0.6, size = 1) +
  scale_color_manual(values = c('blue2', 'gray30', 'red2')) +
  theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.position = c(0.26, 0.92)) +
  theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent')) +
  geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.25) +
  geom_hline(yintercept = -log(0.05, 10), color = 'gray', size = 0.25) +
  labs(x = 'log2 Fold Change', y = '-log10 p-value', color = NA) +
  xlim(-5, 5)

ggsave('volcano_p.png', volcano_p, width = 5, height = 6)
#谢谢指导

结果图片

image.png image.png image.png