要学习如何更好地让AI写代码（别问我是怎么知道的）

#导入数据-------------------------------------------------------------------------------------------
rm(list=ls())
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]
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)
kidney_results <- rain(t(kidney),period = 24, deltat = 4, nr.series = 2,
                       peak.border = c(0.3,0.7), verbose = FALSE)
kidney_select <- row.names( kidney_results[order(kidney_results$pVal)[1:2000], ])
SCN_select <- row.names(SCN_results[order(SCN_results$pVal)[1:1000],])
fibro_select <- row.names(fibro_results[order(fibro_results$pVal)[1:1000],])
library(ggvenn)
vennlist1 <- list(kidney_rhythm = kidney_select,
                  SCN_rhythm = SCN_select,
                  fibroblast_rhythm = fibro_select)
p <- ggvenn(vennlist1,
            stroke_size=0.5,  
            set_name_size=3.3, 
            digits = 0, 
            show_percentage = F 
)
p       
x <- intersect(fibro_select,SCN_select)
y <- intersect(kidney_select,fibro_select)
screen <- merge(as.data.frame(setdiff(x,y)),SCN, by.x=1,by.y=0)
row.names(screen) <- screen$`setdiff(x, y)`
screen <- screen[,-1]
write.csv(screen,file = "screen.csv")
screen_cilia <- merge(cilia_mouse,screen,by.x = 1,by.y = 0)
screen_cilia <-rename(screen_cilia,gene_name=`gene_exchange$mmusculus_homolog_associated_gene_name`)
row.names(screen_cilia) <- screen_cilia$`gene_exchange$mmusculus_homolog_associated_gene_name`
screen_cilia <- screen_cilia[,-1]
write.csv(screen_cilia,file = "screen_cilia.csv")

#所筛选纤毛相关基因表达节律变化可视化--------------------------------------------------------
library(tidyr)
long_data <- gather(screen_cilia, 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---------------------------------------------------------------------
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")