cmap的大部分转录相关数据都是保存为GCT/GCTx格式的，后者是前者的二进制版本。本篇将了解，学习具体格式内容，以及读取方法

1、GCT 格式

• GCT (Gene Cluster Text)：简单来说就是行为基因，列为样本的表达矩阵；同时分别储存行与列的metadata注释信息。
• 如下图是储存有6个样本的10个基因的表达数据的矩阵，同时还包含有许多其它类型数据
  （1）第1行：GCT格式的版本号1.3
  （2）第二行：四个数字分别表示有10个基因，6个样本；基因有2列注释信息，样本有5行注释信息。
  （3）rids与cids标识区域分别表示基因名、样本名；
  （4）rids与chds标识区域分别表示基因注释信息类型、样本注释信息类型。
  （5）data matrix、row metadata、column metadata则储存具体的表达信息与注释信息。
  

2、R包cmapR处理GCT格式数据

• 这个包虽然叫cmapR，但实际与cmap数据库没有太大关联，可以认为是一个处理GCT格式的工具包

# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("cmapR")
library(cmapR)
# my_ds = ds
# write_gct(my_ds, "my_ds")
# write_gctx(my_ds, "my_ds")

• （1）解析、读入本地的GCT/GCTx格式文件

my_ds = parse_gctx("my_ds_n272x978.gct")
class(my_ds)
# [1] "GCT"
# attr(,"package")
# [1] "cmapR"

• （2）"GCT"对象主要结构组成

#mat：表达矩阵
my_ds@mat[1:4,1:4]
#           CPC006_A549_6H:BRD-U88459701-000-01-8:10 CPC020_A375_6H:BRD-A82307304-001-01-8:10
# 200814_at                                    -0.5236                                   0.4525
# 222103_at                                    -0.3721                                  -0.0186
# 201453_x_at                                   0.5784                                   0.3282
# 204131_s_at                                   0.5688                                   1.4957
#           CPC020_HT29_6H:BRD-A82307304-001-01-8:10 CPC020_PC3_24H:BRD-A82307304-001-01-8:10
# 200814_at                                    -0.1483                                  -0.6449
# 222103_at                                    -0.3630                                  -0.2050
# 201453_x_at                                  -0.2848                                  -0.4270
# 204131_s_at                                   0.2366                                   0.6773

# rid：基因名
str(my_ds@rid)
#chr [1:978] "200814_at" "222103_at" "201453_x_at" "204131_s_at" "200059_s_at" "205067_at" "213702_x_at" "214435_x_at" ...

# cid：样本名
str(my_ds@cid)
#chr [1:272] "CPC006_A549_6H:BRD-U88459701-000-01-8:10" "CPC020_A375_6H:BRD-A82307304-001-01-8:10" ...

# rdesc：基因注释
head(my_ds@rdesc[,1:4])
#             is_bing is_lm pr_gene_id pr_gene_symbol
# 200814_at         1     1       5720          PSME1
# 222103_at         1     1        466           ATF1
# 201453_x_at       1     1       6009           RHEB
# 204131_s_at       1     1       2309          FOXO3
# 200059_s_at       1     1        387           RHOA
# 205067_at         1     1       3553           IL1B

# cdesc：样本注释
head(my_ds@cdesc[,1:4])
#                                             brew_prefix cell_id distil_cc_q75 distil_nsample
# CPC006_A549_6H:BRD-U88459701-000-01-8:10   CPC006_A549_6H    A549          0.18              4
# CPC020_A375_6H:BRD-A82307304-001-01-8:10   CPC020_A375_6H    A375          0.46              5
# CPC020_HT29_6H:BRD-A82307304-001-01-8:10   CPC020_HT29_6H    HT29          0.14              4
# CPC020_PC3_24H:BRD-A82307304-001-01-8:10   CPC020_PC3_24H     PC3          0.57              5
# CPC006_A375_6H:BRD-U88459701-000-01-8:10   CPC006_A375_6H    A375          0.20              7
# CPC006_A375_24H:BRD-U88459701-000-01-8:10 CPC006_A375_24H    A375          0.75              5

• （3）解析、读入本地的GCT/GCTx格式文件子集。适用于处理大体积的GCTx二进制文件

#仅读入colmeta信息
col_meta <- read_gctx_meta("my_ds_n272x978.gctx", dim="col")
dim(col_meta) #272 16
#仅读入rowmeta信息
row_meta <- read_gctx_meta("my_ds_n272x978.gctx", dim="row")
dim(row_meta) #978 6

#仅读取前10个样本的相关信息
my_ds_10_columns <- parse_gctx("my_ds_n272x978.gctx", cid=1:10)
dim(my_ds_10_columns@mat)
# [1] 978  10
my_ds_100_rows <- parse_gctx("my_ds_n272x978.gctx", rid=1:100)
dim(my_ds_100_rows@mat)
# [1] 100 272

• （4）取GCT对象子集

my_ds = parse_gctx("my_ds_n272x978.gct")
# [1] 978 272
vemurafenib_ds3 <- subset_gct(my_ds,
                              cid=which(col_meta$pert_iname=="vemurafenib"))
# [1] 978 140

• （5）从头创建GCT对象

m <- matrix(stats::rnorm(100), ncol=10)
rownames(m) <- letters[1:10]
colnames(m) <- LETTERS[1:10]
(my_new_ds <- new("GCT", mat=m))
rdesc <- data.frame(id=letters[1:10], type=rep(c(1, 2), each=5))
cdesc <- data.frame(id=LETTERS[1:10], type=rep(c(3, 4), each=5))
(my_new_ds <- new("GCT", mat=m, rdesc=rdesc, cdesc=cdesc))

• （6）储存GCT对象为GCT/GCTx格式文件

# GCT file
write_gct(my_ds, "my_ds")

# GCTx file
write_gctx(my_ds, "my_ds")

# 转为SummarizedExperiment 对象格式
(se <- as(ds, "SummarizedExperiment"))

cmapR还提供由对应的docker版本，具体使用方法可见下次docker专门学习的笔记~

3、grp、gmt、gmx数据格式

3.1 GRP format

#Read a GRP file and return a vector of its contents
cmapR::parse_grp()

#Write a vector to a GRP file
cmapR::write_grp()

• 本质是只有一列元素的文本文件，没有行名/列名，以换行符为分隔符；
• 元素可以是基因名，样本名...

#genelist为基因名字符串·
write.table(genelist,sep = "\n",file = "myrid.grp",
            row.names = F, col.names = F, quote = F)

3.2 GMT format

#Read a GMT file and return a list
cmapR::parse_gmt()

#Read a GMT file and return a list
cmapR::write_gmt()

• 用于储存基因集的文本格式；
• 一行代表对一个基因集的记录，以\t为分隔符；
• （从左往右）每个第一个元素为基因集名，第二个元素为对该基因集的描述，第三个元素及之后为该基因集里的所有基因
• clusterProfiler包的read.gmt()函数可支持读取，解析gmt文件
  

3.3 GMX format

#Read a GMX file and return a list
cmapR::parse_gmx()

• 同样是用于储存基因集的文本格式；相当于GMT格式的转置过的~
• 一列代表对一个基因集的记录，以\t为分隔符；
• （从上往下）每个第一个元素为基因集名，第二个元素为对该基因集的描述，第三个元素及之后为该基因集里的所有基因
• clusterProfiler包的read.gmt()函数可支持读取，解析gmt文件