同理心

在小丫畫圖交付的一個代碼項目中，需要先從XENA下載一個表達量數(shù)據(jù):https://toil.xenahubs.net/download/tcga_RSEM_gene_tpm.gz

樣本大概是10,5,35個, 考慮到人類的基因大概有2w多個，那么這就是一個10000 X 20000的大樣本數(shù)據(jù)，鑒于這還是一個TPM，數(shù)據(jù)類型是浮點型，文件解壓縮之后就是4.61G, 如果全部加載到R語言中，大部分的電腦估計都受不了

library(pryr)
test <- data.table::fread("./tcga_RSEM_gene_tpm.gz")
object_size(test)
# 5.11 GB

考慮到并非所有數(shù)據(jù)都是我們所需要的，是否可以只讀取部分的數(shù)據(jù)呢？原作者的解決方案是通過R調(diào)用命令行的方式，提取部分數(shù)據(jù)，然后讓R語言進行加載。

system命令

可是大部分人的操作系統(tǒng)都是Windows，所有運行的時候就會報錯，能不能就用戶R語言解決這個問題呢？當然可以，只要你認真讀過read.table的那么多參數(shù)，你就會知道他的那么多參數(shù)并不是裝飾用的。

讀取前幾行

讓我們先學習一個簡單的參數(shù)nrows, 他的作用就是讀取前N行，知道它之后，那就不需要去調(diào)用head了

headtcga <- read.table("./tcga_RSEM_gene_tpm",
                       sep = "\t",
                       stringsAsFactors = FALSE,
                       nrow = 1)

效果就是讀取第一行，構(gòu)建一個數(shù)據(jù)框，然后將其轉(zhuǎn)成向量。但既然目標是向量，其實還有另一種實現(xiàn)方案，readLines讀取的就是一個字符串，然后將其分隔成向量即可。

headtcga <- readLines("tcga_RSEM_gene_tpm", n =1)
headtcga <- strsplit(headtcga, split="\t")[[1]]

讀取指定列

讀取指定列會稍微困難一些，因為colClasses不太好理解。R語言在用read.table讀取數(shù)據(jù)的時候其實做了很多事情，有一件事情就是負責確認每一列的數(shù)據(jù)類型，R語言需要根據(jù)不同數(shù)據(jù)類型進行內(nèi)存分配。

如果你想實現(xiàn)讀取指定列，那么你就得自己去設置每一列的數(shù)據(jù)類型。如果哪些列不需要，就將其它的數(shù)據(jù)類型定義為NULL，R語言就會忽略它。

讀取代碼如下：

cat(paste0("Begin at ", Sys.time(),"\n"))
first_5_rows <- read.table("./tcga_RSEM_gene_tpm", nrows = 5,
                           stringsAsFactors = FALSE, 
                           header = FALSE,
                           skip = 1,
                           check.names = FALSE)
classes <- sapply(first_5_rows, class)
 # targetnum 你需要讀取的列
classes[-targetnum] <- rep("NULL", length(classes) - length(targetnum)) #將非目標列定義為NULL
classes[1] <- "character" # 加上第一列
# 讀取文件(跳過第一行)
targetCancerTPM <- read.table("tcga_RSEM_gene_tpm",  
                   sep= "\t", 
                   skip = 1,
                   colClasses = classes)
colnames(targetCancerTPM) <- tcgasample[targetnum]
targetCancerTPM[1:3, 1:3]
cat(paste0("End at ", Sys.time(),"\n"))

如果僅讀取我們需要的列的話，最終只消耗了500M的內(nèi)存，相對于之前的5G內(nèi)存，減少了將近10倍。

讀取指定行和指定列

這就是需要對文件進行逐行讀取解析了，我用readLines造了一個輪子，函數(shù)名為read_part,目前能用的參數(shù)為

• file: 輸入的文件路徑，支持.gz文件
• rows: 讀取指定行, 比如說1:100, 就是前100行。當為-1時則是讀取所有行
• rows: 讀取指定列, 比如說c(1,3,4,5,6), 就是1,3,4,5,6列。當為-1時則是讀取所有列
• comment.char = "#", 會把"#"開頭的行忽略掉，這個參數(shù)我還需要考慮下是否保留。

# 函數(shù)目標:
# 讀取文件中的指定行和指定列
# 不包括注釋行
read_part <- function(file, rows = 1, columns = -1, sep = "\t",
                      stringsAsFactors = FALSE,
                      header = FALSE,
                      check.names = FALSE, 
                      comment.char = "#", ...){
  dfl <- list()
  if (grepl("gz$", file)){
    con <- gzfile(file, open = "rb")
  } else{
    con <- file(file, open = "r")
  }
  
  i <- 0
  j <- 1
  repeat{
    
    rec <- readLines(con, 1)
    if (length(rec) == 0) break
    i <- i + 1
 
    # 當rows = -1時, 會讀取所有行 
    # 超過目標行時停止讀取
    if (i > max(rows) & rows != -1) break  
    # 不考慮注釋行
    if (grepl(comment.char, rec )) next
    if ( ! i %in% rows & rows != -1) next
    
    items <- strsplit(rec, split = sep, fixed = TRUE)[[1]]
    if ( columns == -1){
      select_cols <- items
    } else{
      select_cols <- items[columns]
    }
    #print(select_cols)
    dfl[[j]] <- select_cols
    j <- j + 1
    
    
  }
  close(con) 
  df <- do.call(rbind, dfl)
  return(df)
}