本人使用的是R下载表达谱数据，从官网下载的临床数据。
2022年TCGA更新后，只能使用新版本的R进行下载，而且不能使用install.packages来安装包了。
使用这段代码配置包的安装环境
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(version = "3.14")
可以看到R都更新到4.3了，我还在用4.1，虽然也能用，但是会有提示，建议直接用4.3版本哈，别像我那么懒
成功之后再运行
library(TCGAbiolinks)
旧版本会有更新提示，好大一堆包😓，请在输入框中填“n”，不然要给你更新到天荒地老去
#利用这段代码可以查看所有TCGA支持的癌症种类的缩写
TCGAbiolinks:::getGDCprojects()$project_id
###下载表达谱数据
#设置下载路径，不知道现在工作路径的，请用“getwd()"查看路径
setwd('/Volumes/TCGA')
#指定下载肺腺癌的数据
cancer_type = 'TCGA-BRCA' #最好一次只下载一个类型，如果一次要下多个癌种，可以用cancer_type = c('TCGA-LUSC','TCGA-BRCA','TCGA-LUAD')
##下载RNA-seq表达谱数据
#设定任务
query <- GDCquery(project = cancer_type,
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "STAR - Counts", #现在只有这一种了，文件里面有FPKM
experimental.strategy = 'RNA-Seq')
#下载数据
GDCdownload(query, method = "api", files.per.chunk = 20) #下载经常会中断，所以喜欢分解成20-50个的小包，每次指定路径一样，下次重启任务可以识别到已下载的包，不会重复下载
#整理数据并存储为R对象
# expr = GDCprepare(query=query) #总是提示错误，最多就到70%，然后就死了，我也就不走这条路了
result = query[[1]][[1]] #用view看一看，id/file-id与文件夹名对应，cases是TCGA编号
write.csv(result,file = paste0(cancer_type,'_query_result.csv')) #保存到本地
###读取已经下载好的TCGA-LUAD数据
library(rjson)
setwd('/Volumes/Starry/NTAN1/TCGA')
#读取metadata文件，格式是json,文件中包含文件名与样本名的对应关系，用于后续与基因名的整合。
metadata <- jsonlite::fromJSON("metadata.cart.2023-09-27 (1).json")
#提取metadata中associated_entities列中的第一个元素为sample id
sample_id <- sapply(metadata$associated_entities,function(x){x[,1]})
#提取file_name
file_name <- metadata$file_id
#创建一个包含样本名和文件名对应关系的数据框
sample_file <- data.frame(sample_id,file_name)
write.csv(sample_file,'sample_id与file_name对应表.csv')
#可以查看前六行
head(sample_file)
#批量读取gdc_download文件夹下的所有gene_counts.tsv格式结尾的基因表达文件
count_file <- list.files("/Volumes/Starry/NTAN1/TCGA/GDCdata/TCGA-BRCA",pattern="*star_gene_counts.tsv",
recursive = TRUE)
#使用strsplit函数分割文件名，分割为文件夹名和tsv结尾的文件名
file_name <- strsplit(count_file,split="/")
#提取出file_name的第6个元素，即tsv文件
file_name <- sapply(file_name,function(x){x[4]}) #注意判断第几个/后是文件名
#创建一个60660行，0列的矩阵，并处理成数据框
data <- data.frame(matrix(nrow=60660,ncol=0))
#通过for循环一次递归读入文件。
for (i in 1:length(file_name)){ #从第一个文件读到最后一个文件
path = paste0("/Volumes/Starry/NTAN1/TCGA/GDCdata/TCGA-BRCA//",count_file[i]) #读取文件的路径
dat<- read.delim(path,fill = TRUE,header = FALSE,row.names = 1) #读取路径下的文件，设置第一列为行名。
colnames(dat)<-dat[2,] #设置第2行为列名
dat <-dat[-c(1:6),] #去掉前6行
dat <- dat[3] #取出表达量对应的这一列，这里选择的是counts类型的基因表达量
colnames(dat) <- sample_file$sample_id[which(sample_file$file_name==file_name[i])]
#将结果data中列名设置为sample_id
data <- cbind(dat,data) #将得到的data与创建好的dat合并
}
table(is.na(colnames(data))) #查看列名是否转换完整,如果有NA，就是metadata下载没对
write.csv(data,'乳腺癌1231例TCGA表达谱原始表格.csv',row.names = TRUE) #文件很大，没事儿不要轻易用excel打开
# setwd('/Users/luoliping/Documents/NTAN1')
# data = read.csv('肺癌939例TCGA原始表格.csv',header = T)
#得到的数据中基因名没有转换，所以要转换基因名，将ensemble ID 转换为gene symbol
#先读取路径下的第一个文件，获取gene symbol
rawdt<- as.matrix(read.delim(path,fill = TRUE,header = FALSE,row.names = 1))
#去掉前6行,并提取第一列为基因名
gene_name <- rawdt[-c(1:6),1]
#合并gene_name与data
rawdt1 = cbind(gene_name,data)
write.csv(rawdt1,'乳腺癌1231例TCGA表达谱原始表格含genename.csv',row.names = TRUE)
#以genename去除重复值，aggregate函数先按照基因名分类，重复的基因名选择保留表达量最大的。
rawdt2 <- aggregate(.~gene_name, data=rawdt1, max) #速度很慢，要等一会儿
write.csv(rawdt2,'乳腺癌1231例TCGA表达谱原始表格含genename-去重.csv',row.names = TRUE)
#将gene_name列设为行名
rownames(rawdt2) <- rawdt2[,1]
#去除第一列
rawdt3 <- rawdt2[,-1]
write.csv(rawdt3,'乳腺癌1231例TCGA表达谱原始表格不含genename-去重.csv',row.names = TRUE)
#----------------------------------分为normal和tumor矩阵------------------------------
sample <- colnames(rawdt2)
normal <- c()
tumor <- c()
for (i in 1:length(sample)){
if((substring(colnames(rawdt2)[i],14,15)>=10)){ #14、15位置大于等于10的为normal样本
normal <- append(normal,sample[i])
} else {
tumor <- append(tumor,sample[i])
}
}
tumor_matrix <- rawdt2[,tumor]
normal_matrix <- rawdt2[,normal]
normal_matrix = cbind(rawdt2$gene_name,normal_matrix)
#写入文件
write.csv(tumor_matrix,'breast_tumor_1119_genename.csv',row.names = TRUE)
write.csv(normal_matrix,'breast_normal_114_genename.csv',row.names = TRUE)
#----------------------------------筛选配对样本-----------------------------------
# 导入必要的包
library(dplyr)
library(stringr)
# setwd('/Users/luoliping/Downloads/NTAN1/TCGA')
# exp = read.csv('肺癌1788例TCGA表达谱原始表格含genename-去重.csv',header =T)
# exp = exp[,-1] #去掉gene_name列
exp = rawdt2 #导入数据
#第14-15位为样本类型标记，<10为肿瘤，>10为正常组织
table(stringr::str_sub(colnames(exp),14,15) < 10) #查看一下大概分组情况
#分组样本
exp_nor = exp[,str_sub(colnames(exp),14,15) > 10] #筛选正常样本
exp_tum = exp[,str_sub(colnames(exp),14,15) <= 10] #筛选肿瘤样本
#提取患者前12位编号，方便配对，两种样本前12位编号一致
patient_tum = str_sub(colnames(exp_tum),1,12) #肿瘤的
patient_nor = str_sub(colnames(exp_nor),1,12) #正常的
patient_mix = intersect(patient_tum, patient_nor) #取两个的交集
#截取exp_tum数据框中每个列名的前12个字符，并判断截取结果是否存在于patient向量
k_tum = str_sub(colnames(exp_tum),1,12) %in% patient_mix #肿瘤的
table(k_tum)
k_nor = str_sub(colnames(exp_nor),1,12) %in% patient_mix #正常的
table(k_nor)
#只取能匹配上的
exp_tum = exp_tum[,k_tum]
exp_nor = exp_nor[,k_nor]
#合并肿瘤和正常样本
exp2 = cbind(exp_tum,exp_nor)
# 随便挑一个id号，判断两列数据是否一致，如果不一致就要检查哪里错了
if(all(rawdt2$`TCGA-34-5232-01A-21R-1820-07` == exp2$`TCGA-34-5232-01A-21R-1820-07`)) {
print("两列数据一致")
} else {
print("两列数据不一致")
}
#按照列名的字母顺序进行排序，方便查看配对样本
sorted_exp2 <- exp2[, order(names(exp2))]
# 提取列名的前12位，忽略特殊符号
col_names <- gsub("[^a-zA-Z0-9]", "", substr(names(sorted_exp2), 1, 12))
# 统计重复的列名个数
table(col_names) #如果有次数为1的则不正确，>2说明采样不止一次，后面可以人工筛选去除
table(table(col_names))
sorted_exp = cbind(exp$gene_name,sorted_exp2) #加上gene_name列
write.csv(sorted_exp,'TCGA_breast_113例配对样本.csv')
#----------------------------------提取临床信息-----------------------------------
#提取临床信息
setwd('/Volumes/Starry/NTAN1/TCGA')
#从metadata数据中提取出case_id，用来匹配
library("rjson")
json <- jsonlite::fromJSON("metadata.cart.2023-09-27 (1).json")
View(json)
entity_submitter_id <- sapply(json$associated_entities,function(x){x[,1]})
case_id <- sapply(json$associated_entities,function(x){x[,3]})
sample_case <- t(rbind(entity_submitter_id,case_id)) #case_id可以匹配
#从clinical数据中提取出临床信息
library(data.table)
library(dplyr)
clinical = fread('/Users/luoliping/Downloads/clinical.cart.2023-09-27/clinical.tsv')
colnames(clinical)
dim(clinical)
# 删除值为 "--" 超过50%的列
clinical1 <- clinical %>% select_if(~ mean(. != "'--", na.rm = TRUE) > 0.5)
# 与metadata匹配
clinical_matrix <- merge(sample_case,clinical1,by="case_id",all.x=T)
table(clinical_matrix$project_id) #查看一下是否有NA
# 删除重复行
clinical_matrix1 <- distinct(clinical_matrix, entity_submitter_id, .keep_all = TRUE) #按entity_submitter_id
# clinical_matrix1 <- as.data.frame(clinical_matrix[duplicated(clinical_matrix$case_id),]) #按case_id
write.csv(clinical_matrix1,'clinical_matrix_0927.csv',row.names = TRUE)
