# импортироватьданныенабор1
data1 <- read.csv("PRJNA423456_FPKM.csv", header = TRUE) #Набор данных здесь представляет собой смоделированный набор данных. Измените его на свой собственный набор данных для анализа.
# Импорт данных. Эпизод 2.
данные2 <- read.csv("PRJNA777728_FPKM.csv", header = TRUE) #Набор данных здесь представляет собой смоделированный набор данных. Измените его на свой собственный набор данных для анализа.
# импортироватьданныенабор3
data3 <- read.csv("PRJNA6878999_FPKM.csv", header = TRUE) #Набор данных здесь представляет собой смоделированный набор данных. Измените его на свой собственный набор данных для анализа.
# импортироватьданныенабор4
data4 <- read.csv("PRJNA590000_FPKM.csv", header = TRUE) #Набор данных здесь представляет собой смоделированный набор данных. Измените его на свой собственный набор данных для анализа.
#### слитьданныенабор
library("sva") #Нетпрограммное программное обеспечение, пожалуйста, сначала установите программное обеспечение обеспечение
library("tidyverse") #Нетпрограммное программное обеспечение, пожалуйста, сначала установите программное обеспечение обеспечение
merge_eset1 <- inner_join(data1, data2, by = "Gene.ID", suffix = c(".data1", ".data2")) #X — это имя первого столбца. Имя гена должно быть помещено в первый столбец, а не стать именем строки.
merge_eset1 
merge_eset2 <- inner_join(data3, data4, by = "Gene.ID", suffix = c(".data3", ".data4"))
merge_eset <- inner_join(merge_eset1, merge_eset2, by = "Gene.ID", suffix = c(".merge_eset1", ".merge_eset2"))
merge_eset
rownames(merge_eset) <- merge_eset$Gene.ID
merge_eset <- merge_eset[,-1]  #Помещаем первый столбец с именами генов в качестве имен строк
dim(merge_eset)    #Viewданные измерения
exp <- as.matrix(merge_eset)
dimnames <- list(rownames(exp),colnames(exp))
data <- matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
dim(data)
class(data) #Вернем «матрицу» "array"
write.csv(data,"mrna_nocombat.csv")
#Отфильтровать гены с низкой экспрессией
Expr <- data[rowSums(data)>1,]
Expr=log2(Expr+1)

Анализ #PCA визуализирует эффект кластеризации отдельных наборов данных перед их разделением на пакеты.
# Здесь мы сначала получаем информацию о группе
group_list <- data.frame(
  sample = colnames(exp_all), c(rep("PRJNAX423456_A",3),rep("PRJNA423456_B",3),rep("PRJNA777728_A",3),rep("PRJNA777728_C",3),
                                rep("PRJNA777728_B",3),rep("PRJNA6878999_A",10),rep("PRJNA6878999_B",9),rep("PRJNA590000_A",3),rep("PRJNA590000_C",3),
                                rep("PRJNA590000_B",3)))
rownames(group_list) <- group_list$sample
colnames(group_list)[2] <- "dataset"
group <- factor(group_list$dataset)
View(group_list)
##Загрузка информации об образце, включая партии
data2 <- read.csv("All_data.csv",header = T)
data2[,2] <- as.factor(data2$Type )##type представляет различные биологические процессы
data2[,3] <- as.factor(data2$Batch )##batch представляет различную информацию о пакете
row.names(data2) <- data2$Sample
View(data2)
data2

библиотека(FactoMineR)##Если она недоступна, сначала установите ее
library(factoextra)
pdf(file = "PCA_before1.pdf",width = 7,height = 6)
pre.pca <- PCA(t(exp_all),graph = FALSE)
fviz_pca_ind(pre.pca,
             geom= "point",
             col.ind = data2$Batch,
             addEllipses = TRUE,
             legend.title="Group"
              )
dev.off()

#------------------Распаковка------------------#
#### Расчет пакетных эффектов с помощью пакета sva
library(sva)
exp_all_combat <- ComBat(exp_all, batch = group_list$dataset) # партия — это информация о партии
# Ознакомьтесь с результатами удаления пакетных эффектов
library(tinyarray)
draw_pca(exp = exp_all_combat, group_list = group)

#Проведите калиброванный анализ PCA
pdf(file = "9_PCA_after.pdf",width = 7,height = 6)
pre.pca <- PCA(t(exp_all_combat),graph = FALSE)
fviz_pca_ind(pre.pca,
             geom= "point",
             col.ind = data2$Batch,
             addEllipses = TRUE,
             legend.title="Group")
dev.off()

setwd("F:\\RRA Analysis-2024-03-29\\4_RRA Analysis-2\\Anagen_Catagen")
#1. Считайте набор данных GSE.
files=c("PRJNA590000_A_C.DESeq2.csv",
        "PRJNA779755_A_C.DESeq2.csv")
upList=list()
downList=list()
allFCList=list()
for(i in 1:length(files)){
  inputFile=files[i]
  rt=read.csv(inputFile)
  header=unlist(strsplit(inputFile,"_"))
  downList[[header[1]]]=as.vector(rt[,1]) 
  upList[[header[1]]]=rev(as.vector(rt[,1]))
  fcCol=rt[,1:2]   colnames(fcCol)=c("Gene",header[[1]]) 
  allFCList[[header[1]]]=fcCol 
}
View(allFCList)

#2. Объедините все значения FC.
mergeLe=function(x,y){
  merge(x,y,by="Gene",all=FALSE)}
newTab=Reduce(mergeLe,allFCList)
View(newTab)
rownames(newTab)=newTab[,1]
newTab=newTab[,2:ncol(newTab)]
newTab[is.na(newTab)]=0
dim(newTab)

#3. Получить регулируемые гены.
library(RobustRankAggreg)
upMatrix = rankMatrix(upList)
View(upMatrix)
dim(upMatrix)
upAR = aggregateRanks(rmat=upMatrix)
View(upAR)
dim(upAR)
colnames(upAR)=c("Name","Pvalue")
upAdj=p.adjust(upAR$Pvalue,method="bonferroni")
#help(p.adjust)
upXls=cbind(upAR,adjPvalue=upAdj)
upFC=newTab[as.vector(upXls[,1]),]
upXls=cbind(upXls,logFC=rowMeans(upFC))
write.table(upXls,file="up.xls",sep="\t",quote=F,row.names=F)
upSig=upXls[(upXls$Pvalue<0.01 & upXls$logFC> log2(2)),]
dim(upSig)
write.table(upSig,file="upSig_2.xls",sep="\t",quote=F,row.names=F)
#3. Получить гены с пониженной регуляцией.
downMatrix = rankMatrix(downList)
downAR = aggregateRanks(rmat=downMatrix)
colnames(downAR)=c("Name","Pvalue")
downAdj=p.adjust(downAR$Pvalue,method="BH")
downXls=cbind(downAR,adjPvalue=downAdj)
downFC=newTab[as.vector(downXls[,1]),]
downXls=cbind(downXls,logFC=rowMeans(downFC))
write.table(downXls,file="down.xls",sep="\t",quote=F,row.names=F)
downSig=downXls[(downXls$Pvalue<0.01 & downXls$logFC< -log2(2)),]
dim(downSig)
write.table(downSig,file="downSig_2.xls",sep="\t",quote=F,row.names=F)
#Интегрируйте гены с повышающей и понижающей регуляцией
allSig = rbind(upSig,downSig)
colnames(allSig)
allSig = allSig[,c("Name","logFC")]
View(allSig)
write.table(allSig,file = 'allSign_2.xls',sep = '\t',quote = F)
#Отображение результатов 20 лучших генов с повышающей и понижающей регуляцией logFC.tiff
hminput=newTab[c(as.vector(upSig[1:20,1]),as.vector(downSig[c(1:20),1])),]
library(pheatmap)
tiff(file="logFC2_1_2.tiff",width = 15,height = 20,units ="cm",compression="lzw",bg="white",res=400)
pheatmap(hminput,display_numbers = TRUE,fontsize_row=10,fontsize_col=12,
         color = colorRampPalette(c("blue","white", "red"))(200),
         cluster_cols = FALSE,cluster_rows = FALSE, angle_col = 45)
dev.off()