2022-07-25 RNA-Seq 火山图与PCA

normCout <- function(filtered){
  filtered <- as.matrix(filtered)
  df <- data.frame(condition=colnames(filtered))
  sample <- unique(df$condition)[2]
  dds <- DESeqDataSetFromMatrix(filtered,colData = df,design = ~condition)
  dds <- estimateSizeFactors(dds)
  dds_DE <- DESeq(dds)
  res_DE <-results(dds_DE,alpha = 0.05, contrast = c("condition",unique(df$condition)[2],"CTRL"))
  reOrdered <- res_DE[order(res_DE$padj),]
  resSig <- subset(reOrdered,padj < 0.05)
  resSigUp <- subset(resSig,log2FoldChange>1)
  resSigAll <- subset(resSig,(log2FoldChange>1|log2FoldChange<(-1)))
  # write.csv(resSigAll,file = paste0("3.RNA-Seq/2.featureCount/CTRL_",sample,".DEG.csv"),quote = F)
  # for_volcano <- data.frame("log2FoldChange"=res_DE$log2FoldChange,
  #                           "padj"=res_DE$padj,
  #                           "trend" = rep( "NoDiff", length(res_DE$log2FoldChange))) %>% na.omit()
  # up_sig_indices <- intersect(which(for_volcano$log2FoldChange >1), which(for_volcano$padj <0.05))
  # down_sig_indices <- intersect(which(for_volcano$log2FoldChange <(-1) ), which(for_volcano$padj < 0.05))
  # for_volcano[up_sig_indices, "trend" ]<- "Up"
  # for_volcano[down_sig_indices, "trend" ] <- "Down "
  # for_volcano$trend <- as.factor(for_volcano$trend)
  # for_volcano$padj <- -log10(for_volcano$padj)
  # ggplot(for_volcano,aes(x=log2FoldChange,y=padj,colour=trend) ) + geom_point(size=I(1.2))+
  #   scale_color_manual(values = c("blue","black","red"))+
  #   geom_vline(xintercept = c(1,-1),lty=2,size=I(0.4),colour = "grey11")+
  #   geom_hline(yintercept = c(-log(x=0.05, base=10)),lty=2, size=I(0.4), colour = "grey11")+
  #   my_clear_theme+
  #   theme(panel.background = element_rect(colour = "black", size=1,fill = "white"),panel.grid = element_blank())+
  #   labs(x="1og2(FoldChange)", y="-log10(Pvalue) ")
  # ggsave(paste0("3.RNA-Seq/2.featureCount/CTRL_",sample,".DEG.png"),width = 8,height = 6,units = "in",dpi=300)
  # rld <- rlog(dds)
  # pcaData <- plotPCA(rld, intgroup=c("condition"), returnData=T)
  # percentVar <- round(100*attr(pcaData, "percentVar"))
  # pcaData$condition <- factor(pcaData$condition,levels = c("CTRL",sample))
  # print(sample)
  # ggplot(pcaData, aes(PC1, PC2, color=condition, shape=condition)) + 
  #   geom_point(size=3) + 
  #   xlab(paste0("PC1: ", percentVar[1], "% variance")) + 
  #   ylab(paste0("PC2: ", percentVar[2], "% variance"))+
  #   geom_hline(aes(yintercept=0),linetype="dashed")+
  #   geom_vline(aes(xintercept = 0),linetype="dashed")+
  #   theme(panel.background = element_rect(colour = "black", size=1,fill = "white"),
  #         panel.grid = element_blank(),legend.title = element_blank(),
  #         legend.text = element_text(size = 25,face = "bold"),
  #         legend.key = element_blank(),
  #         legend.key.size = unit(0.6, 'inch'),legend.background = element_blank(),
  #         axis.text.x = element_text(face="bold", size=25,color = "black"),
  #         axis.text.y = element_text(face="bold",   size=25,color = "black"),
  #         axis.title.y = element_text(face = "bold",size = 28),
  #         axis.title.x = element_text(face = "bold",size = 28))
  # ggsave(paste0("3.RNA-Seq/2.featureCount/CTRL_",sample,".PCA.png"),width = 8,height = 6,units = "in",dpi=300)
  return(resSigAll)
}


count <- read.table("3.RNA-Seq/2.featureCount/All.name.txt",skip = 1,header = T,row.names = 1)
colnames(count)[6:14] <- c("CTRL","CTRL","CTRL","*","*","*","*","*","*")
CTRL_RAD <- count[,c(6:8,12:14)] %>% filter() %>% `colnames<-`(c("CTRL","CTRL","CTRL","*","*","*"))
CTRL_BLEO <- count[,c(6:11)] %>% filter() %>% `colnames<-`(c("CTRL","CTRL","CTRL","*","*","*"))