iCluster分析是基于多组学数据的整合分类算法。更多知识分享请到 https://zouhua.top/。

Rscript

  1. #!/usr/bin/R
  2. library(argparser)
  3. library(dplyr)
  4. library(tibble)
  5. library(ggplot2)
  6. library(data.table)
  7. library(iCluster)
  9. rm(list = ls())
  10. options(stringsAsFactors = F)
  11. options(future.globals.maxSize = 1000 * 1024^2)
  13. grp <- c("S1", "S2")
  14. grp.col <- c("#6C326C", "#77A2D1")
  16. # parameter input
  17. parser <- arg_parser("iCluster function") %>%
  18.     add_argument("-p", "--phen", 
  19.         help = "phenotype") %>%
  20.     add_argument("-c", "--copyNumber", 
  21.         help = "copyNumber matrix") %>%
  22.     add_argument("-g", "--geneExp", 
  23.         help = "gene Expression matrix") %>%
  24.     add_argument("-m", "--methylation", 
  25.         help = "DNA methylation matrix")  %>%
  26.     add_argument("-r", "--RPPA", 
  27.         help = "protein_RPPA")  %>%
  28.     add_argument("-o", "--out", 
  29.         help = "result with director", default = "./")
  31. args <- parse_args(parser)
  33. # prepare for function 
  34. phen <- fread(args$p)                                  
  35. copyNumber <- fread(args$c)
  36. geneExp <- fread(args$g)
  37. methylation <- fread(args$m)
  38. protein_RPPA <- fread(args$r)
  39. out <- args$o    
  41. get_profile <- function(dataset = copyNumber,
  42.                         metadata = phen, 
  43.                         tag = "copyNumber"){
  45.   # dataset = protein_RPPA
  46.   # metadata = phen
  47.   # tag = "protein_RPPA"
  49.   sid <- intersect(phen$Barcode, colnames(dataset)) 
  50.   res <- dataset %>% dplyr::select(c("V1", sid)) %>% 
  51.     mutate(Type=tag) %>%
  52.     mutate(Name=paste(V1, Type, sep = "_")) %>%
  53.     dplyr::select(Name, V1, Type, everything()) %>%
  54.     dplyr::select(-c("V1", "Type")) %>%
  55.     column_to_rownames("Name") %>%
  56.     t()
  57.   return(res)
  58. }
  60. datasets <- list(copyNumber=get_profile(dataset = copyNumber, tag = "copyNumber"),
  61.                  geneExp=get_profile(dataset = geneExp, tag = "geneExp"),
  62.                  methylation=get_profile(dataset = methylation, tag = "methylation"),
  63.                  protein_RPPA=get_profile(dataset = protein_RPPA, tag = "protein_RPPA"))
  65. print(names(datasets))
  68. # icluster 
  69. fit <- iCluster(datasets = datasets, k=2, lambda=rep(0.2, 4), max.iter = 50, epsilon = 1e-3)
  70. #plotiCluster(fit=fit, label=rownames(datasets[[2]]))
  71. #compute.pod(fit)
  73. phen_new_icluster <- inner_join(phen,
  74.               data.frame(SampleID=rownames(datasets[[2]]), Cluster=fit$clusters), 
  75.              by = c("Barcode"="SampleID")) %>% 
  76.   dplyr::select(Barcode, Cluster, everything()) %>%
  77.   mutate(Cluster=paste0("S", Cluster))
  79. name1 <- paste0(out, "phenotype_cluster_iCluster.csv")
  80. write.csv(phen_new_icluster, file = name1, row.names = F)
  83. fit2 <- iCluster2(datasets = datasets, k=2, lambda=list(0.2, 0.2, 0.2, 0.2),
  84.                  max.iter = 50, verbose = TRUE)
  86. if(0){
  87.     phen_new_icluster2 <- inner_join(phen,
  88.                 data.frame(SampleID=rownames(datasets[[2]]), Cluster=fit2$clusters), 
  89.                 by = c("Barcode"="SampleID")) %>% 
  90.     dplyr::select(Barcode, Cluster, everything()) %>%
  91.     mutate(Cluster=paste0("S", Cluster))
  93.     name2 <- paste(out, "phenotype_cluster_iCluster2.csv")
  94.     write.csv(phen_new_icluster2, file = name2, row.names = F)
  95. }
  97. save(fit, fit2, file = "iCluster_fit.RData")

Run

Rscript iCluster.R -p ../../Result/phenotype/common_survival_data.tsv \
 -c ../../Result/profile/copyNumber_filter.tsv \
 -g ../../Result/profile/geneExp_filter.tsv \
 -m ../../Result/profile/methylation_filter.tsv \
 -r ../../Result/profile/protein_RPPA_filter.tsv \
 -o ./

Reference

  1. iCluster package

参考文章如引起任何侵权问题，可以与我联系，谢谢。