在 Posit Cloud（以前称为 RStudio）上以节省内存的方式合并大型 TCGA .tsv 文件

Question

我有 448 个 .tsv 文件，其中包含从癌症基因组图谱基因组数据共享 (GDC) 门户网站下载的基因表达数据 (RNAseq)。

这些文件有 60666 行和 9 列。然后有一个元数据文件，第一列有

case_id

，第二列有

file_name

。

所以我的

merge_RNA

函数所做的是将元数据文件和所有448个

.tsv

文件所在的目录作为输入。然后如果文件的

basename

与

file_name

匹配，则按前三列“gene_id”、gene_name、“gene_type”和第八列“fpkm_unstranded”合并文件。

所以合并后的矩阵应该有 60661 行和 451 列，其中 gene_id 是行，case_id 是列。

但是，我遇到了内存问题，因为

read.csv

函数将整个 .tsv 文件读入内存。

注意：我使用的是 Posit Cloud（以前称为 RStudio），所以我的计算能力有限（只有 16 GB 和 4 个 CPU）。

# Define the function to merge all RNAseq quantification files into one dataframe

merge_rna <-function(metadata, fdir){
  filelist <- list.files(fdir, pattern="*.tsv$", 
                     recursive = TRUE, full.names=TRUE)
  for (i in 1:length(filelist)){
    iname <- basename(filelist[i])
    isamplename <- metadata[metadata$file_name==iname, "sample"]
    idf <- read.csv(filelist[i], sep="\t", skip=1, header=TRUE)
    # remove first 4 rows
    remove <- 1:4
    idf_subset <- idf[-remove, c("gene_id", "gene_name", "gene_type", "fpkm_unstranded")]
    rm(idf)
    names(idf_subset)[4] <- isamplename
    #print(dim(idf_subset))
    if (i==1){
      combined_df <- idf_subset
      rm(idf_subset)
    } else {
      combined_df <- merge(combined_df, idf_subset, by.x='gene_name', by.y="gene_name", all=TRUE)
      rm(idf_subset)
    }
  }

# use gene_id as row names and remove gene_id column
rownames(combined_df) <- combined_df$gene_id
combined_df <- combined_df[,-which(names(combined_df) %in% c("gene_id"))]
return(combined_df)
}

rnaCounts <-  merge_rna(metaMatrix, "/home/r1816512/TSV_subset")

我提供的代码似乎可以工作，但它占用了大量内存。请告诉我如何使用更少的内存。也许我们不必将整个文件读入内存？

我尝试了 ChatGPT，但尽管我尽了最大努力对其进行微调以重写我的代码并给我一个使用更少内存的版本，但它给我的所有代码都产生了错误。关于 ChatGPT 在不到 5 分钟内生成代码但调试该代码需要 24 小时的模因确实是真的。

ChatGPT 建议减少内存使用的一种方法是分块读取 TSV 文件，而不是一次全部读取。这可以使用

readr

包来完成，特别是

read_tsv_chunked()

功能。

ChatGPT 还建议使用较少内存的另一种方法是分块处理数据，而不是使用

data.table

包将整个文件读入内存并分块处理数据。

但是，它提供给我的同时使用

readr

和

data.table

函数的代码都不起作用。

Answer 1

使用

data.table::fread()

拉入每个tsv，只选择你需要的列，使用

rbindlist()

，合并一次回到

metaMatrix

以获得样本ID，然后（可选）dcast

files = list.files("../gdc", pattern="*.tsv", full.names = T, recursive = T)
targ_cols = c("gene_id", "gene_name", "gene_type","fpkm_unstranded")
result = rbindlist(
  setNames(
    lapply(files, \(f) fread(f,header=F, skip=6, select = c(1,2,3,8),col.names = targ_cols)),
    basename(files)
),idcol = "file_name")
dcast(
  result[metaMatrix, on="file_name"][,file_name:=NULL],
  gene_id+gene_name+gene_type~sample, value.var = "fpkm_unstranded"
)

输出：

Key: <gene_id, gene_name, gene_type>
                  gene_id  gene_name      gene_type      S1     S10      S2      S3      S4      S5      S6      S7      S8      S9
                   <char>     <char>         <char>   <num>   <num>   <num>   <num>   <num>   <num>   <num>   <num>   <num>   <num>
    1: ENSG00000000003.15     TSPAN6 protein_coding  2.3871 22.5869  6.4368 15.3576 15.7226  8.5207 12.0493  7.3320  4.6332  4.0896
    2:  ENSG00000000005.6       TNMD protein_coding  0.0148  0.0262  0.0104  0.0000  0.0131  0.0431  0.0000  0.0000  0.0000  0.0274
    3: ENSG00000000419.13       DPM1 protein_coding 47.2737 35.4827 37.7895 18.0376 31.9880 21.8047 33.9062 28.4055 35.2305 38.7990
    4: ENSG00000000457.14      SCYL3 protein_coding  1.0667  1.5429  1.0221  2.5562  1.7876  1.3778  1.8601  1.4731  1.7868  1.0997
    5: ENSG00000000460.17   C1orf112 protein_coding  1.0218  1.4619  1.2092  3.6312  3.0915  1.0475  0.8116  1.1096  3.5854  1.0718
   ---                                                                                                                             
60656:  ENSG00000288669.1 AC008763.4 protein_coding  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0203  0.0000  0.0000
60657:  ENSG00000288670.1 AL592295.6         lncRNA  2.9298  3.4835  1.7472  1.7875  1.8250  1.1271  2.2591  2.1847  1.9592  1.3517
60658:  ENSG00000288671.1 AC006486.3 protein_coding  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
60659:  ENSG00000288674.1 AL391628.1 protein_coding  0.0023  0.0000  0.0048  0.0016  0.0060  0.0050  0.0107  0.0084  0.0085  0.0042
60660:  ENSG00000288675.1 AP006621.6 protein_coding  0.1816  0.2758  0.1551  0.3328  0.3099  0.3220  0.2454  0.2414  0.1583  0.3484

输入（元矩阵）

structure(list(file_name = c("6d1ffd0c-ba13-463e-b5c4-20d2dab6294d.rna_seq.augmented_star_gene_counts.tsv", 
"054bcc02-0463-4a2f-94a7-b268c15a47d4.rna_seq.augmented_star_gene_counts.tsv", 
"5bf0c22a-d47f-4b0d-9d43-760fcaa1b42c.rna_seq.augmented_star_gene_counts.tsv", 
"d1c3a993-a9b3-4938-949a-a38793413139.rna_seq.augmented_star_gene_counts.tsv", 
"d1d934ef-2e6d-4d28-8418-9a0b5ef66c0c.rna_seq.augmented_star_gene_counts.tsv", 
"eb102931-4d5d-469d-82bb-187c6e81d6fb.rna_seq.augmented_star_gene_counts.tsv", 
"0269349e-aab4-459d-8a8b-99534c01c90e.rna_seq.augmented_star_gene_counts.tsv", 
"963cb78b-6480-453f-9f62-d92c205f034a.rna_seq.augmented_star_gene_counts.tsv", 
"de53dd9c-416a-492f-89bf-a79866fd3f2d.rna_seq.augmented_star_gene_counts.tsv", 
"b92a6b85-da4b-4ce0-bc1e-fd8908f1edba.rna_seq.augmented_star_gene_counts.tsv"
), sample = c("S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", 
"S9", "S10")), class = "data.frame", row.names = c(NA, -10L))