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. 将转置的结果重新转变为数据框