昨日的改进

#导入数据-------------------------------------------------------------------------------------------
setwd("C:/Users/mrb06/Desktop")
kidney <-read.csv("annotated_kidney_transcriptome.csv")
rownames(kidney) <- kidney$X;kidney <- kidney[,-1]
SCN <-read.csv("annotated_SCN_transcriptome.csv")
rownames(SCN) <- SCN$X;SCN <- SCN[,-1]
fibro <-read.csv("annotated_fibroblast_transcriptome.csv")
rownames(fibro) <- fibro$X;fibro <- fibro[,-1]
pancrease <-read.csv("annotated_pancrease_transcriptome.csv")
rownames(pancrease) <- pancrease$X;pancrease <- pancrease[,-1]
#纤毛相关基因鼠同源基因获取-------------------------------------------------------------------------
cilia <- read.csv("CiliaCarta.csv")
library(biomaRt)
human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
gene_exchange <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name','mmusculus_homolog_associated_gene_name'), 
                       filters = 'ensembl_gene_id', 
                       values = cilia$Ensembl.Gene.ID,
                       mart = human)
cilia_mouse <- as.data.frame(gene_exchange$mmusculus_homolog_associated_gene_name)
#rhythm analysis: RAIN--------------------------------------------------------------------- 
library(rain)
SCN_results <- rain(t(SCN),period = 24, deltat = 4, 
                    peak.border = c(0.3,0.7), verbose = FALSE)
fibro_results <- rain(t(fibro),period = 24, deltat = 4, 
                      peak.border = c(0.3,0.7), verbose = FALSE)
pancrease_results <- rain(t(pancrease),period = 24, deltat = 4, 
                          peak.border = c(0.3,0.7), verbose = FALSE)
kidney_results <- rain(t(kidney),period = 24, deltat = 4, nr.series = 2,
                       peak.border = c(0.3,0.7), verbose = FALSE)
rm(list=c("kidney","SCN","fibro","pancrease","cilia"))
SCN_select <- row.names(SCN_results[order(SCN_results$pVal)[1:1500],])
fibro_select <- row.names(fibro_results[order(fibro_results$pVal)[1:1500],])
pancrease_select <- row.names(pancrease_results[order(pancrease_results$pVal)[1:1500],])
kidney_select <- row.names( kidney_results[order(kidney_results$pVal)[1:3000], ])
#venn图绘制----------------------------------------------------------------------------------------
library(ggvenn)
vennlist1 <- list(kidney_rhythm = kidney_select,
                  SCN_rhythm = SCN_select,
                  fibroblast_rhythm = fibro_select,
                  pancrease_rhythm = pancrease_select)
p <- ggvenn(vennlist1,
            stroke_size=0.5,  # 集合圆圈的线宽
            set_name_size=3.3, # 集合名称的文本大小
            digits = 0, # 小数点后保留位数
            show_percentage = F # 是否展示每一部分所占的百分比
)
p       
screen <- merge(
  as.data.frame(
    setdiff(
      setdiff(intersect(fibro_select,SCN_select),kidney_select),
      pancrease_select)),
  SCN, by.x=1,by.y=0)
write.csv(screen,file = "screen.csv")
screen_cilia <- merge(cilia_mouse,screen,by.x = 1,by.y = 1)
names(screen_cilia)[1] <- "gene_name"
write.csv(screen_cilia,file = "screen_cilia.csv")
#所筛选纤毛相关基因表达节律变化可视化--------------------------------------------------------
library(tidyr)
cilia_high <- screen_cilia[c(3,8),]
cilia_low_1 <- screen_cilia[-c(3,8),][c(1:4),]
cilia_low_2 <- screen_cilia[-c(3,8),][c(5:8),]
long_data <- gather(cilia_low_2, key = 'variable', value = 'value', -gene_name)
ggplot(long_data, 
       aes(x = variable, y = value, group = gene_name, color = gene_name)) + 
  geom_line(size = 1) +  
  geom_point(size = 2, shape = 21, fill = "white", color = "black") +
  labs(title = "Gene Expression Levels Over Time Points",
       x = "Time Points (Variable)",
       y = "Expression Value",
       color = "Gene Name") + 
  theme_minimal() +
  theme(legend.position="right")
#trend analysis: Mfuzz---------------------------------------------------------------------
row.names(screen) <- screen$`setdiff(setdiff(intersect(fibro_select, SCN_select), kidney_select), pancrease_select)`
screen <- screen[,-1]
screen <- screen[,c(1:6,1)]
require(Biobase)
bio.screen.Mfuzz <- ExpressionSet(as.matrix(screen))
require("Mfuzz")
bio.screen.Mfuzz <- standardise(
  filter.std(
    fill.NA(
      filter.NA(bio.screen.Mfuzz, thres = 0.25), 
      mode = "mean"), 
    min.std = 0)
)
m1 <- mestimate(bio.screen.Mfuzz)
c1 <- mfuzz(bio.screen.Mfuzz, c=4, m=m1)
mfuzz.plot(bio.screen.Mfuzz, cl=c1, 
           mfrow = c(2,2),time.labels = seq(0,24,4))
write.csv(as.data.frame(c1$cluster),file = "screen cluster.csv")