如何正确格式化以下R代码以读取GZ文件?

cvxl0en2  于 2023-05-26  发布在  其他
关注(0)|答案(2)|浏览(218)

我试图打开在此链接上找到的两个文件(在V7下):https://gtexportal.org/home/datasets,通过R命名为“基因TPM”和“转录TPM”,这样我就可以对它们进行分析。下面是我的代码:

geneTPM = read.table(gzfile("C:/Users/ual-laptop/Downloads/GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz"))
transcriptTPM = read.table(gzfile("C:/Users/ual-laptop/Downloads/GTEx_Analysis_2016-01-15_v7_RSEMv1.2.22_transcript_tpm.txt.gz"))

但是当我运行代码时,第一行给出了一个错误消息,说“第1行没有11690个元素”,第二行只是给出了一个空行(甚至没有一个“+”或错误消息),使它看起来像是在运行,但它只是保持这种状态几个小时。我知道这些文件很长,所以预计需要一段时间,但我认为我在这里做错了一些事情(第一行肯定是这样,因为有一个错误消息,但idk如何解决它)。

pxq42qpu

pxq42qpu1#

可以使用data.table::fread()直接读取gz文件,但需要安装R.utils包。

data.table::fread(
  "GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz",
  nrows = 10,
  select = 1:6)
)

输出:

Name  Description GTEX-1117F-0226-SM-5GZZ7 GTEX-111CU-1826-SM-5GZYN GTEX-111FC-0226-SM-5N9B8
               <char>       <char>                    <num>                    <num>                    <num>
 1: ENSG00000223972.4      DDX11L1                  0.10820                  0.11580                  0.02104
 2: ENSG00000227232.4       WASH7P                 21.40000                 11.03000                 16.75000
 3: ENSG00000243485.2   MIR1302-11                  0.16020                  0.06433                  0.04674
 4: ENSG00000237613.2      FAM138A                  0.05045                  0.00000                  0.02945
 5: ENSG00000268020.2       OR4G4P                  0.00000                  0.00000                  0.00000
 6: ENSG00000240361.1      OR4G11P                  0.00000                  0.10510                  0.07637
 7: ENSG00000186092.4        OR4F5                  0.00000                  0.00000                  0.00000
 8: ENSG00000238009.2 RP11-34P13.7                  0.23920                  0.07206                  0.20940
 9: ENSG00000233750.3       CICP27                  0.18690                  0.39590                  0.08927
10: ENSG00000237683.5   AL627309.1                  9.46700                  8.64300                  2.38600
    GTEX-111VG-2326-SM-5N9BK
                       <num>
 1:                  0.02329
 2:                  8.17200
 3:                  0.00000
 4:                  0.03260
 5:                  0.00000
 6:                  0.00000
 7:                  0.00000
 8:                  0.23180
 9:                  0.19770
10:                 31.61000

(在完整的数据集中,实际上有56K行和11.7K列)

k = data.table::fread(
  "GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz"
)
dim(k)

[1] 56202 11690
ivqmmu1c

ivqmmu1c2#

这些是相对较大的文件(transcripts.tpm.gz是~2Gb,gene.tpm.gz是~ 800 Mb压缩),并且它们在文件的前两行中包含元数据,这可能会导致将数据导入R的问题(尽管我在langtang的答案中看到data.table::fread()处理元数据行没有问题,+1)。
一个可能的解决方案是使用readr package中的read_tsv()来解压缩文件,指定分隔符(制表符)并跳过元数据行,例如

library(readr)
library(tictoc)

Sys.setenv("VROOM_CONNECTION_SIZE" = 2500000) # increase buffer size

tic() # see how long the process takes

geneTPM <- read_tsv("~/Downloads/GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz",
                    skip = 2, progress = TRUE)
#> indexed 17.00B in  0s,
#> 2.16MB/s indexed 30.00MB in  0s,
#> 138.31MB/s indexed 32.50MB in  0s,
#> 140.64MB/s indexed 35.00MB in  0s,
#> 138.09MB/s indexed 37.50MB in  0s,
#> 138.83MB/s indexed 40.00MB in  0s,
#> 140.06MB/s indexed 42.50MB in  0s,
#> 141.16MB/s indexed 45.00MB in  0s, 
...
#> 106.88MB/s indexed 1.01GB in  9s,
#> 106.97MB/s indexed 1.01GB in  9s,
#> 75.94MB/s indexed 2.81GB in 37s,
#> 75.95MB/s indexed 2.82GB in 37s,
#> 75.93MB/sindexed 1.00TB in 37s,
#> 26.97GB/s                                                                              
#> Rows: 56202 Columns: 11690
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr     (2): Name, Description
#> dbl (11688): GTEX-1117F-0226-SM-5GZZ7, GTEX-111CU-1826-SM-5GZYN, GTEX-111FC-...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

str(geneTPM)
#> spc_tbl_ [56,202 × 11,690] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
#>  $ Name                          : chr [1:56202] "ENSG00000223972.4" "ENSG00000227232.4" "ENSG00000243485.2" "ENSG00000237613.2" ...
#>  $ Description                   : chr [1:56202] "DDX11L1" "WASH7P" "MIR1302-11" "FAM138A" ...
#>  $ GTEX-1117F-0226-SM-5GZZ7      : num [1:56202] 0.1082 21.4 0.1602 0.0505 0 ...
...
#>   [list output truncated]
#>  - attr(*, "spec")=
#>   .. cols(
#>   ..   Name = col_character(),
#>   ..   Description = col_character(),
#>   ..   `GTEX-1117F-0226-SM-5GZZ7` = col_double(),
...
#>   ..   `GTEX-ZXG5-0005-SM-57WCN` = col_double()
#>   .. )
#>  - attr(*, "problems")=<externalptr>

toc() # how long did it take?
#> 234.907 sec elapsed

创建于2023-05-23带有reprex v2.0.2

相关问题