1. 整理表达矩阵
下载的文件是按样本存放的,每个tsv文件中都记录着一个样本的基因表达量,需要将所有tsv文件合并,得到所有样本的基因表达量的表格。
准备
setwd("D:/R/CHOL") if(!require("rjson"))install.packages("rjson",update = F,ask = F) library(rjson) library(tidyverse)
1. 设置下载数据路径(需修改)
2. 判断式安装“rjson”包
获取每个样品id和相应tsv文件名的对应关系:
样品id和相应tsv文件的对应关系存放在metadata文件中。
json <- jsonlite::fromJSON("metadata.cart.2022-09-26.json") View(json) sample_id <- sapply(json$associated_entities,function(x){x[,1]}) file_sample <- data.frame(sample_id,file_name=json$file_name)
1. 读取Metadata文件中每个样品的信息(需修改matadata文件的名称)
3. 提取所有样品id,如:TCGA-W5-AA2X-01A-11R-A41I-07
4. 获取样品id和相应tsv文件名的对应关系:
合并所有样品id的基因表达量
每个tsv文件中,gene_id的顺序是都一致的,可以直接cbind,不会导致顺序错乱。
count_file <- list.files("expdata",pattern = "*.tsv",recursive = TRUE) count_file_name <- strsplit(count_file,split="/") %>% sapply(function(x){x[2]}) matrix = data.frame(matrix(nrow = 60660,ncol = 0)) for (i in 1:length(count_file)){ data<- paste0("expdata/",count_file[i]) %>% read.delim(fill = TRUE,header = FALSE,row.names = 1) colnames(data) <- data[2,] data <- data[-c(1:6),] data <- data[3] colnames(data) <- file_sample$sample_id[which(file_sample$file_name == count_file_name[i])] matrix <- cbind(matrix,data) } write.csv(matrix,"COUNT_matrix.csv",row.names = TRUE)
1. 获取expdata文件夹下的所有tsv表达文件的“路径+文件名”
2-3. 在count_file中分割出文件名
4. 新建一个空白数据框,行数等于tsv文件中的基因数
5-13. 循环读取每一个tsv文件
10. 取出unstranded列(第3列),即count数据,对应其它数据
14. 输出表达矩阵文件:(行名为样品id,列名为Ensembl ID,值为基因的表达量)
设置Gene Symbol为表达矩阵的列名(上一步得到的列名是Ensembl ID)
因为“matrix”和“gene_name”均来源于tsv文件,所以可以直接cbind,不会导致顺序错乱。
data <- paste0("expdata/",count_file[1]) %>% read.delim(fill = TRUE,header = FALSE,row.names = 1) %>% as.matrix() gene_name <-data[-c(1:6),1] matrix0 <- cbind(gene_name,matrix) matrix0 <- aggregate( . ~ gene_name,data=matrix0, max) rownames(matrix0) <- matrix0[,1] matrix0 <- matrix0[,-1]
1-3. 读取一个tsv文件(因为所有的tsv文件中“gene_name”的顺序都相同,读取一个即可)
4. 提取该tsv文件中的第一列“gene_name”
5. 将“gene_name”这一列插入到表达矩阵“matrix”中
6. 将“gene_name”列去除重复的基因,保留每个基因最大表达量结果。aggregate函数中,“. ~ gene_name”表示以gene_name为分类,计算其他列的最大值。
7. 将“gene_name”列设为行名
8. 去掉“gene_name”列
过滤掉在很多样本中表达量都为0的基因
matrix0 = matrix0[apply(matrix0,1,function(x)sum(x > 1) >= 20), ]
1. 判断“matrix0”中的每个基因,是否在至少20个样本中>1。保留符合条件的基因,剔除不符合条件的基因。
样本数“>=20”根据总样本量修改,可以是数值,比如大于正常样本数“>=40”,也可以是百分比,如大于总样本数的75%“>=0.75*(ncol(matrix0))”
将样本分为normal组和tumor组
sample <- colnames(matrix0) normal <- c() tumor <- c() for (i in 1:length(sample)){ if((substring(colnames(matrix0)[i],14,15) > 10)){ normal <- append(normal,sample[i]) } else { tumor <- append(tumor,sample[i]) } } tumor_matrix <- matrix0[,tumor] normal_matrix <- matrix0[,normal] group_list = ifelse(as.numeric(str_sub(colnames(matrix0),14,15)) < 10,'tumor','normal') %>% factor(levels = c("normal","tumor")) table(group_list)
1. 从表达矩阵中提取列名
2-3. 创建“normal”和“tumor”空白向量
4-10. 样本id中14、15位置大于10的为normal样本,否则为tumor样本
11-12. 将总的表达矩阵分为normal和tumor矩阵
13. 将正常样品和肿瘤样品分开,并设置因子,便于后续分析
14. 查看正常样品和肿瘤样品的数量
2. 整理临床信息
临床信息中没有分期grade,看后续分析怎么解决
https://www.jingege.wang/2022/07/16/%e6%96%b0%e7%89%88tcga%e6%95%b0%e6%8d%ae%e4%b8%8b%e8%bd%bd%e5%8f%8ar%e8%af%ad%e8%a8%80%e6%95%b4%e5%90%88%e4%bb%a3%e7%a0%81/
先对单个文件进行检查:
library(XML) result <- xmlParse("./clinical/0f133012-23ef-4237-acfd-47b132b99775/nationwidechildrens.org_clinical.TCGA-W5-AA2Q.xml") rootnode <- xmlRoot(result) rootsize <- xmlSize(rootnode) print(rootnode[1]) print(rootnode[2]) xmldataframe <- xmlToDataFrame(rootnode[2]) head(t(xmlToDataFrame(rootnode[2])))
2. 读取单个xml文件,需修改文件目录和文件名
3. 将根节点从 xml 文件中移除。
4. 查找根中的节点数。
5. 查看第一个节点
6. 查看第二个节点,确认第二个节点为临床信息的保存位置
7. 将xml中第二个节点转换为数据框
8. 转置并查看该数据框
读取所有临床数据文件:
获得所有临床数据文件路径
xmls = dir("clinical/",pattern = "*.xml$",recursive = T)
读取目录为“clinical/”,读取模式为以“*.xml”结尾的文件,recursive表示列表是否应重新进入目录
构建读取函数:
td = function(x){ result <- xmlParse(file.path("clinical/",x)) rootnode <- xmlRoot(result) xmldataframe <- xmlToDataFrame(rootnode[2]) return(t(xmldataframe)) }
读取并整合所有临床信息
cl = lapply(xmls,td) cl_df <- t(do.call(cbind,cl)) cl_df[1:3,1:3] clinical = data.frame(cl_df) clinical[1:4,1:4]
1. 对于xmls中的每一个文件,执行td函数,输出每一个病人的临床信息,并整合到列表中
2. 合并列表中的每一个元素为数据框,并将数据框转置(转置的输出结果为矩阵)
4. 将转置的结果重新转变为数据框