1. dir <- c("import", "outport", "Rscript", "Rdata", "figures")
    2. for (i in 1:length(dir)){
    3.   if(!dir.exists(dir[i])) dir.create(dir[i])
    4. }
    6. get_disease_info = function(diseasename){
    7.   disease <- diseasename
    9.   if(!dir.exists(paste0("import/", disease)))
    10.     dir.create(paste0("import/", disease))
    11.   if(!dir.exists(paste0("import/", disease, "/clinical")))
    12.     dir.create(paste0("import/", disease, "/clinical"))
    13.   if(!dir.exists(paste0("import/", disease, "/expdata")))
    14.     dir.create(paste0("import/", disease, "/expdata"))
    16.   # 下载临床信息和表达矩阵
    17.   system(paste0("./gdc-client.exe download -m ./import/",disease,"/gdc_manifest_cl.txt -d ./import/",disease,"/clinical/"))
    18.   system(paste0("./gdc-client.exe download -m ./import/",disease,"/gdc_manifest_exp.txt -d ./import/",disease,"/expdata/"))
    20.   dir_cl <- paste0("import/", disease, "/clinical/")
    21.   dir_exp <- paste0("import/", disease, "/expdata/")
    23.   ## 临床信息整理
    24.   library(XML)
    25.   xmls <- dir(dir_cl, pattern = "*.xml$", recursive = T)
    26.   cl <- list()
    27.   for(i in 1:length(xmls)){
    28.     result <- xmlParse(paste0(dir_cl,xmls[[i]]))
    29.     rootnode <- xmlRoot(result)
    30.     cl[[i]] <- xmlToDataFrame(rootnode[2])
    31.   }
    32.   clinical <- do.call(rbind, cl)
    33.   if(!dir.exists(paste0("Rdata/", disease)))
    34.     dir.create(paste0("Rdata/", disease))
    35.   save(clinical, file = paste0("Rdata/", disease, "/clinical.Rdata"))
    37.   ## 表达矩阵整理
    38.   counts <- dir(dir_exp, pattern = "*.htseq.counts.gz$", recursive = T) 
    39.   expdata <- list()
    40.   for(i in 1:length(counts)){
    41.     expdata[[i]] <- read.table(paste0(dir_exp,counts[[i]]),row.names = 1,sep = "\t")
    42.   }
    43.   expdata <- do.call(cbind,expdata)
    45.   # 添加列名
    46.   meta <- jsonlite::fromJSON(paste0("import/", disease, "/metadata.json"))
    47.   ids <- meta$associated_entities
    48.   ID <- sapply(ids,function(x){x[,1]})
    49.   file2id <- data.frame(file_name = meta$file_name,
    50.                         ID = ID)
    51.   counts <- stringr::str_split(counts,"/",simplify = T)[,2]
    52.   file2id <- file2id[match(counts, file2id$file_name),]
    53.   colnames(expdata) <- file2id$ID
    55.   # 将行名换成gene symbol
    56.   # TCGA有自己的id转换方法
    57.   library(rtracklayer)
    58.   gtf <- import("import/gencode.v22.annotation.gtf") 
    59.   gtf <- as.data.frame(gtf)
    60.   gtf_gene <- gtf[gtf$type=="gene",]
    61.   an <- gtf_gene[,c("gene_name","gene_id","gene_type")]
    62.   expdata <- expdata[rownames(expdata) %in% an$gene_id,]
    63.   an <- an[match(rownames(expdata),an$gene_id),]
    64.   k <- !duplicated(an$gene_name);table(k)
    65.   an <- an[k,]
    66.   expdata <- expdata[k,]
    67.   rownames(expdata) <- an$gene_name
    69.   Group <- ifelse(stringr::str_sub(colnames(expdata),14 ,15)<10, "tumor", "normal")
    71.   # 过滤低表达基因(玄学)
    72.   number <- 0.5*ncol(expdata)
    73.   expdata <-  expdata[apply(expdata, 1, function(x) sum(x > 10) >number), ]
    74.   save(expdata, Group, file = paste0("Rdata/", disease, "/expdata.Rdata"))
    75. }
    77. get_disease_info('ACC')