Hello, everyone, I am trying to reproduce the lncRNA get from the paper "Portraying breast cancer with lncRNAs" by Olivier et. al. Here I came to the reproducing step:

![enter image description here][1]


I don't know whether my code is doing right, can someone give me your comment? 



    # load RAW files
    RAW_files = list.files(path = "./RAWsets", pattern = "_RAW.tar")
    # RNA.list = list of RAW file after processed by Frozen Robust Multiarray Analysis
    RNA.list = list()
    # iterate through the .CEL files
    for (RAWstr in RAW_files) {
      print(RAWstr)
      # extract those raw .tar files
      RAWdirectory = paste("./RAWsets", substr(RAWstr, 1, nchar(RAWstr) - 8), sep = "/")
      tardir = paste("./RAWsets", RAWstr, sep = "/")
      untar(tardir, exdir = RAWdirectory, verbose = F)
      # extract those extracted .CEL files
      cels = list.files(RAWdirectory, pattern = "CEL.gz", full.names = T)
      sapply(cels, gunzip)
      # read all cel files to data.raw variable
      data.raw = ReadAffy(filenames = substr(cels, 1, nchar(cels)-3), verbose = F, cdfname = "hgu133plus2cdf")
      # frozen robust multiarray analysis
      data.matrix = frma(data.raw)
      # add to RNA.list as described above
      tmp = list(data.matrix)
      names(tmp) = substr(RAWstr, 1, nchar(RAWstr) - 8)
      RNA.list = append(RNA.list, tmp)
    }
    
    ### Combine into 1 matrix of all data sets
    RNAdata = as.matrix(exprs(RNA.list[[1]]))
    for (i in 2:length(RNA.list)) {
      RNAdata = CombineMatrix(RNAdata, as.matrix(exprs(RNA.list[[i]])))
    }
    # normalize between arrays on a data of fRMA before
    RNAdata_scaled = normalizeBetweenArrays(RNAdata)
    # because samples names are GSMxxxxxx.CEL => erase .CEL extension
    colnames(RNAdata_scaled) = substr(colnames(RNAdata_scaled), 1, 9)
    
    # load data from table S1, S2, S3, S7
    load("tables_set.RData")
    # mutual index between total data and tableS1 which have clinical information
    tableS1 = tableS1[match(colnames(RNAdata_scaled), tableS1[, 1]), ]
    # compute batch data using tumor vs normal tissue as covariate
    batch = as.numeric(lapply(tableS1[, 3], function(str) if (str == "normal") return(1) else return(2)))
    # applied Combat algorithm in sva library with default parameters to adjust data for batch effects
    RNAdata_scaled = ComBat(dat = RNAdata_scaled, batch = batch)


Your help is really appreciated!

  [1]: https://image.ibb.co/h3Mvmk/Screenshot_from_2017_08_28_11_23_14.png
Hey, it seems that what you are setting as batch is in reality the covariate value (which is defined as the condition or "variable of interest" in the linear model). This is the code for applying ComBat if you specify the batches and the conditions that you have in your dataset

    sample <- Pheno$Sample.Name
    batch <- Pheno$Batch
    condition <- Pheno$Covariate.1
    pdata <- data.frame(sample, batch, condition)
    modmatrix <- model.matrix(~as.factor(condition), data=pdata)
    combat_data_matrix <- ComBat(dat=yourmatrix, batch=batch, mod=modmatrix)