# 目的：通过starbase预测linrna下游mirna，人和鼠分别预测后，取交集。并预测mirna下游mrna。与gse27170差异基因取交集。
# 固定开头
rm(list = ls())
library(pacman)
pacman::p_load(GEOquery,stringr,dplyr,readxl,stringr,limma,ggplot2,cowplot,ggsci,ggpubr,styler,tibble,DESeq2,ggrepel,reshape2,VennDiagram,
               clusterProfiler,enrichplot) 
library(tidyverse)
# 设定工作目录
path <- "D:\\myresearch\\doctor_project\\starbase_data"
setwd(path)
## -------------------------读入starbase预测的结果-------------------------------------
# lincrna_hsa <- read_excel("ENCORI_hg19_CLIP-seq_miRNA-target_all_PVT1.xls")
lincrna_hsa <- read.table("ENCORI_hg19_CLIP-seq_miRNA-target_all_PVT1.txt",skipNul = T,fill = T,sep = "\t",header = T)
# 提取miRNA名称，去掉物种标记
lincrna_hsa <- lincrna_hsa %>% mutate(mirna = miRNAname %>% str_replace_all("hsa-",""))

lincrna_mmu <- read.table("ENCORI_mm10_CLIP-seq_miRNA-target_all_Pvt1.txt",skipNul = T,fill = T,sep = "\t",header = T)
lincrna_mmu <- lincrna_mmu %>% mutate(mirna = miRNAname %>% str_replace_all("mmu-",""))
# 指定lncrna
lincrna <- "NEAT1"
# 以下代码自动运行，包括鼠基因名大小写转换
lincrna_hsa <- read.table(str_c(c("ENCORI_hg19_CLIP-seq_miRNA-target_all_",lincrna,".txt"),collapse = "")
                          ,skipNul = T,
                          fill = T,
                          sep = "\t",
                          header = T)
# 提取miRNA名称，去掉物种标记
lincrna_hsa <- lincrna_hsa %>% mutate(mirna = miRNAname %>% str_replace_all("hsa-",""))
lincrna_lower <- str_c(c(str_sub(lincrna,1,1),str_sub(lincrna,2) %>% tolower()),collapse = "")
lincrna_mmu <- read.table(str_c(c("ENCORI_mm10_CLIP-seq_miRNA-target_all_",lincrna_lower,".txt"),collapse = ""),
                          skipNul = T,
                          fill = T,
                          sep = "\t",
                          header = T)
lincrna_mmu <- lincrna_mmu %>% mutate(mirna = miRNAname %>% str_replace_all("mmu-",""))

lincrna

lincrna_hsa$mirna
lincrna_mmu$mirna
# 取交集
lincrna_inter <- intersect(lincrna_hsa$mirna,lincrna_mmu$mirna)
lincrna_inter
# 以上为人鼠共同mirna
##------------------------------------预测下雨miRNA------------------------------------
library(multiMiR)
# 设定预测靶点，加入hsa-或者mmu-前缀
lincrna_inter <- lincrna_inter %>% sapply(function(x)str_c("hsa-",x,collapse = ""))
lincrna_inter
# 开始计算，全程联网，mirna,target决定预测方向，空值时作为目标

example3 <- get_multimir(org = "hsa",
                         mirna = lincrna_inter,          
                         # target = "GINS2",
                         table = "predicted",
                         predicted.cutoff = 35,
                         predicted.cutoff.type = "p",
                         predicted.site = "all")
# 提取结果
example3_result <- example3@data

example3_result

example3_sum <- example3@summary
example3_sum %>% head()

# 查看你各个数据库结果数
example3@data$database %>% table()

table(example3_result$mature_mirna_id)

diff_sig1 <- read.csv("D:\\myresearch\\GEO\\GSE27170\\tab\\diff_sig_1.csv",row.names = 1)


mrna_inter_diff <- list()
for(i in lincrna_inter){
  print(i)
  mrna_target <- example3_result[example3_result$mature_mirna_id == i,"target_symbol"]
  mrna_inter_diff[[i]] <- intersect(mrna_target,rownames(diff_sig1 %>% filter(sign == "up")))
  
}

mrna_inter_diff
saveRDS(mrna_inter_diff,file = "result_pvt1,rds")