TCGA及GEO数据挖掘/生物信息学

参考资料：

dplyr库：CRAN: Package dplyr

GEO官网:http://www.ncbi.nlm.nih.gov/geo/

R语言基础：R 基础语法 | 菜鸟教程

R语言常见语法

c(1,2)  : 将多个元素组合成一个向量

strsplit(count_file,split='/') ：字符串切割，返回向量

数据框变量$列名 ：提取某列

sapply(json$associated_entities,function(x){x[,1]}) ：提取某列，并应用函数。

x[,n]  :提取所有行的第n列

x[n,] :提取第n行的所有列

x(1:6)  : 向量的第1个到第6个

data[-c(1:6),]   选取 data 中除了第 1 到第 6 行之外的所有行

as.matrix() 是一个类型转换函数，用于将其他数据结构（如数据框、向量、表格等）转换为矩阵（matrix） 类型。若原数据包含多种类型（如数值 + 字符），转换时会统一为字符型

which :匹配表达式

矩阵：矩阵要求所有元素类型一致

数据框: 是一种特殊的“矩阵”，元素类型可以不一致。

subset：根据条件表达式筛选行。

Question：

啥是注释包？

Warning

好像默认没安装build tool（貌似还没发现有影响，先这样吧）

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
试开URL’https://cran.rstudio.com/bin/windows/contrib/4.3/rjson_0.2.23.zip'
Content type 'application/zip' length 433808 bytes (423 KB)
downloaded 423 KB

环境：

• R 4.3.1

安装相关

R包安装的几种方式

参考资料

1. https://zhuanlan.zhihu.com/p/21069119503
2. 通过Bioconductor/BiocManager安装生物r包详解(问题汇总)-CSDN博客

官方在线安装方式：

install.packages("包名")
library("包名")

通过BiocManager安装：

//这玩意是个啥，好像是判断是否有导入这个包，没有导入才去安装BiocManager包
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
//导入BiocManager包
library("BiocManager")
//BiocManager库里调用install函数来下载需要的包。
BiocManager::install("lima")
library("lima")

dplyr库安装

dplyr 是 R 语言中一个非常流行的数据操作包

• 使用在线安装：package("dplyr")
• 导入库 library(dplyr)

历史版本：Index of /src/contrib/Archive/dplyr

TCGA RNA-seq数据下载与整理

看不懂，大为震撼。

 TCGA 数据库上对应的癌种的全称以及缩写，参考：

(32 封私信 / 84 条消息) TCGA / 癌症简称 / 缩写 / TCGA癌症中英文对照 - 知乎

TCGA数据下载

TCGA官网：portal.gdc.cancer.gov

click builder -- jump -- >  toggle  TCGA Program and some project

click Repository  ， 选择类别和类型

Add all file to cart

download cart and download metadata

QAQ

metadata文件：记录了文件名称与样本名之间的映射关系的信息

cart文件：存储了基因表达的信息

示例

cart: e789877b-6eb5-434c-85b8-a62a99937b98

cart文件中缺少样本名称，通过meta文件可以去关联。

安装所需要的库

//本地有就不用再安装了
install.packages("rjson")
install.packages("tidyverse")

//如果不知道本地有没有，感觉可以Require函数来判断是否下载了。
//如果未安装才会自动安装它。
if (!require("rjson", quietly = TRUE))
    install.packages("rjson")
if (!require("tidyverse", quietly = TRUE))
    install.packages("tidyverse")

引入库

library(rjson)
library(tidyverse)

设置读写路径

setwd("E:\\BaiduNetdiskDownload\\3TCGA")
//setwd("文件的根路径")

读取metadata Json文件并生成表结构

//读到的内容，赋值给变量json
json <- jsonlite::fromJSON("metadata.cart.2023-06-13.json")

获取样本名称及文件名称

// 从JSON 数据的 associated_entities 字段中提取每个实体的第一列数据
//json$associated_entities 提取第一列向量
// function(x){x[,1]} 提取associated_entities列中所有行的矩阵的第一列向量
sample_id <- sapply(json$associated_entities,function(x){x[,1]})
  

输出第一列数据1到10条

sample_id[1:10]

创建一个数据框file_sample，以sample_id向量为列，metadata的file_name向量为列填充。

file_sample <- data.frame(sample_id,file_name=json$file_name)

最后就得到了，每个基因表达cart文件与样本的映射关系。

获取所有cart文件名

列出指定目录下所有文件包括子目录下文件

count_file <- list.files('gdc_download_20230613_075345.926514/',
                         pattern = '*.tsv',recursive = TRUE)

子文件夹下的文件，返回的文件名会包括子文件路径。

因此需要根据/去拆分字符串，拿到第二列的文件名。

count_file_name <- strsplit(count_file,split='/')

要拿到第二列文件夹，从切割后得到字符串数组中提取第二列的所有行数据，重新赋值给count_file_name

count_file_name <- sapply(count_file_name,function(x){x[2]})

合并

#构建一个空的数据框
matrix = data.frame(matrix(nrow=60660,ncol=0))

#逐个读取及合并
for (i in 1:length(count_file)){
//拼接字符串，得到完整路径
  path = paste0('gdc_download_20230613_075345.926514//',count_file[i])   #Counts文件夹名
// read.delim用于读取分隔符分隔的文本文件,例如TSV 文件
  data<- read.delim(path,fill = TRUE,header = FALSE,row.names = 1)
  colnames(data)<-data[2,]
  data <-data[-c(1:6),]
  #3：unstranded,counts；4：stranded_first；5：stranded_second；6：tpm_unstranded；7：fpkm_unstranded；8：fpkm_uq_unstranded
  #data <- data[3]
  data <- data[6]
  #data <- data[7]
  colnames(data) <- file_sample$sample_id[which(file_sample$file_name==count_file_name[i])]
  matrix <- cbind(matrix,data)
}

read.delim接口

read.delim用于读取分隔符分隔的文本文件,例如TSV 文件

• path：要读取的文件路径（需要替换为实际的文件路径字符串）
• fill = TRUE：当行的长度不一致时，自动填充空白值以保持数据框结构整齐
• header = FALSE：表示文件没有表头（第一行不是列名）
• row.names = 1：指定使用文件的第 1 列作为数据框的行名（行索引）

colnames：会返回数据框的所有列名

colnames(data)<-data[2,] 数据框第 2 行数据设置为该数据框的列名。

data <-data[-c(1:6),] : 选取 data 中除了第 1 到第 6 行之外的所有行, 并重新赋值给data。

data <- data[6] 最后将提取了第6列向量作为data

which(file_sample$file_name == count_file_name[i]) 从file_sample矩阵中file_name列向量中找到与当前迭代到的文件名相匹配的那个行，并返回索引位置。

根据which查到的行号，可以提取到对应的样本名，并将样本名作为data列向量的列名。这样就可以得到列名为样本名，每行值为基因表达数据中的第6列的值。

colnames(data) <- file_sample$sample_id[which(file_sample$file_name==count_file_name[i])]

将每次迭代计算得到的data，按照横向列向量合并填充进矩阵中。

matrix <- cbind(matrix,data)

这样就能很清楚的看出每一份基因文件对应的样本。

视频里说是啥，想要以gene_name为行名。我听的云里雾里的，目的是啥不清楚。估计这玩意是个啥医学专有名词简称？更直观点？

这下面相当是把所有值都变成字符串了，并多添加了两列gene_name、gene_type

#转化为gene_symbol
path = paste0('gdc_download_20230613_075345.926514//',count_file[1])
data<- as.matrix(read.delim(path,fill = TRUE,header = FALSE,row.names = 1))
gene_name <-data[-c(1:6),1]
#gene_name[1:10]
matrix0 <- cbind(gene_name,matrix)
#获取基因类型
gene_type <- data[-c(1:6),2]
#gene_type[1:10]
matrix0 <- cbind(gene_type,matrix0)

没完。gene_name列里有重复的值，重复的值不能作为行名。因此需要想办法剔除，仅保留一个。

它这里是留下样本值最大的那个。

matrix0 <- aggregate( . ~ gene_name,data=matrix0, max)   这块处理有点慢，要点时间

• aggregate() 用于数据聚合
• . ~ gene_name 聚合公式：
  • 左侧的 . 代表 "所有其他列"（除了分组列 gene_name 之外的列）
  • 右侧的 gene_name 是分组依据列（按基因名分组）
• data = matrix0 指定要处理的数据对象
• max 是聚合函数，表示对每组计算最大值

示例：

筛选gene_type值为protein_coding的那些行。

matrix0 <- subset(x = matrix0, gene_type == "protein_coding")

将gene_name列设为行名：

rownames(matrix0) <- matrix0[,1]

#转化为导出格式
matrix0 <- matrix0[,-c(1,2)]
matrix1 = data.frame(ID=rownames(matrix0),matrix0)
colnames(matrix1) = gsub('[.]', '-', colnames(matrix1))

#写入，导出
write.table(matrix1,'TCGA_LUSC_TPM.txt', sep="\t", quote=F, row.names = F)

GEO

GEO数据基因注释三种方法

• 获取TXT文档
  • 
• 注释信息Soft文件
  • 
• 注释包

示例

GEOquery：专门用于处理GEO数据包的库。

TXT包

GEO数据包（TXT包）下分三部分，分别是：

第一部分是样本信息：

第二部分上临床信息

第三部分是表达矩阵

获取GEO包中的表达矩阵：exprs(gset[[1]])

获取临床信息并转成表结构或矩阵：pData(gset[[1]])

注释文件是啥，为啥要读？

GEO 注释文件（通常指 GPL 平台注释文件）是用于解释基因表达数据中 “探针 ID” 或 “特征 ID” 含义的关键文件。包含了这些 ID 与基因名（如 TP53）、基因符号（Symbol）的对应关系，能将原始数据 “翻译” 成可理解的基因名称。

目标：从注释文件中提取出探针ID和基因符号。

ids=gpl[,c("ID","gene_assignment")] 只需关心这两个，那就剔除其他的。

很奇怪的是，基因符号在注释文件的gene_assignment列中的第二个位置，因此要想拿到基因符号，需要分割字符串成向量，从向量中取出基因符号。

ids$symbol=trimws(str_split(ids$symbol,'//',simplify = T)[,2])

1. 首先先对gene_assignment每行都进行字符串分割，并取出第二个元素（基因符号）
2. 再通过trimws操作来去除首尾空白的字符，最后得到基因符号。

有些基因探针没有对应的基因符号，会有异常值。需要剔除过滤掉。

ids=ids[ids$symbol != '',]
ids=ids[ids$symbol != '---',]
ids=ids[ids$symbol != '--- ',]

这段 R 代码是按条件筛选数据框行的语法。感觉也可以用which。ids <- ids[which(ids$symbol != ''), ]

最后，要判断平台文件的ID和矩阵（表达矩阵）中的ID是否匹配。ids=ids[ids$probe_id %in% rownames(dat),] %in%用于判断是否匹配。

啥是平台文件？看代码，平台文件应该就是指的注释文件。

用 ids$probe_id 作为行索引，筛选出 dat 中与这些探针 ID 匹配的行 dat=dat[ids$probe_id,]

#合并
dat <- cbind(ids,dat)

#删除重复基因，保留最大值
dat <- aggregate( . ~ symbol,data=dat, max)

Soft包

可以通过GPL@dataTable@table 来从 GPL Soft注释文件 中提取 注释信息表格。

#下载soft注释文件，目录下有会自动读入
GPL=getGEO(filename = 'GPL25318_family.soft.gz',destdir = "F:/3SP/6GEO/2GSE116918-GPL25318") 

#获取注释信息
gpl=GPL@dataTable@table