下载count文件
image.png
使用gdc-client下载count文件和xml文件
###需要下载的文件有两种，一种是.xml格式的临床数据，一种是表达矩阵的数据（还要一个json文件）
library(stringr)
proj="TCGA-BC"
#建好两类文件的文件夹，把下好的数据存到里面
if(!dir.exists("clinical"))dir.create("clinical")
if(!dir.exists("expdata"))dir.create("expdata")
dir()
调用gdc-client下载的命令在windows的terminal终端进行
image.png
整合数据
library(BiocManager)
BiocManager::install("XML")
library(XML)
#.xml的临床数据保存在各个目录下，将它们整合到一起
xmls<-dir("clinical/",pattern = "*.xml$",recursive = T)
cl<-list()
for (i in 1:length(xmls)) {
result<-xmlParse(paste0("clinical/",xmls[[i]]))
rootnode<-xmlRoot(result)
cl[[i]]<-xmlToDataFrame(rootnode[2])
}#cl里面现在存了若干个表格，还要将若干表格整合到一张表里
library(plyr)
information_needed<-c("bcr_patient_barcode","vital_status","days_to_last_followup","days_to_death")#需要的临床信息
#将若干表格中需要的临床信息整合到一张表里
c<-rbind(cl[[1]][information_needed],cl[[2]][information_needed])
for (i in 3:408) {
c<-rbind(c,cl[[i]][information_needed])
}
#这一步是针对生存信息（Alive or Dead）的状态进行一个判断
for (i in 1:408) {
if (c$vital_status[i]=="Alive") {c$Status[i]=1}
if (c$vital_status[i]=="Dead") {c$Status[i]=0}
}
#保存表型数据
write.csv(c,"phenotype_data.csv",quote = F)
#表达数据保存在各个目录下，将它们整合到一张表格中
count_files<-dir("expdata_BC/",pattern = "*.counts.gz$",recursive = T)
exp<-list()
for (i in 1:length(count_files)) {
exp[[i]]<-read.table(paste0("expdata_BC/",count_files[[i]]),row.names = 1,sep = "\t")
}
exp<-do.call(cbind,exp)
dim(exp)
#这样得到的表达矩阵没有列名，加上列名
meta<-jsonlite::fromJSON("metadata.cart.2022-01-06.json")#之前通过metadata下载的.json文件
ids<-meta$associated_entities
#得到file_name和sample_ID的对应矩阵
ID<-sapply(ids,function(x){x[,1]})
file2id<-data.frame(file_name=meta$file_name,ID=ID)
View(file2id)
#修改矩阵的列名
count_files2<-stringr::str_split(count_files,"/",simplify = T)[,2]
table(count_files2 %in% file2id$file_name)
#count_files2的顺序就是列名的顺序，根据它来调整file2id的顺序。
file2id<-file2id[match(count_files2,file2id$file_name),]
colnames(exp)<-file2id$ID
#ENSEMBL号转换成Gene symbol名
BiocManager::install("AnnotationDbi")
BiocManager::install("org.Hs.eg.db")
# 加载包
library("AnnotationDbi")
library("org.Hs.eg.db")
#将ENSG号提取出来
ENSG_name<-rownames(exp)
head(ENSG_name)
#要去掉"ENSG00000000003.13"小数点后面的数字，不然没法转换
ENSG_name_2<-c()
for (i in 1:length(ENSG_name)) {
ENSG_name_2[i]<-substr(ENSG_name[i],1,15)
}
head(ENSG_name_2)
#ENSG转换成基因名
ENSG_symbol<- mapIds(org.Hs.eg.db,
keys=ENSG_name_2,
column="SYMBOL",
keytype="ENSEMBL",
multiVals="first")
symbol_names<-c()
for (i in 1:length(ENSG_name_2)) {
symbol_names[i]<-ENSG_symbol[[i]]
}
head(symbol_names)
#将symbol名添加到exp矩阵中
exp$symbol<-symbol_names
#有一些行没有对应的基因名要去掉
number<-c()
for (i in 1:length(ENSG_name_2)) {
if (is.na(exp[i,430])) { #这里要改一下列数
number<-c(number,i)
}
}
exp1<-exp
exp1<-exp1[-number,]
dim(exp1)
#对于重复的基因取均值最大的一行
rowMeans<-apply(exp1,1,function(x) mean(as.numeric(x),na.rm=T))
exp2<-exp1[order(rowMeans,decreasing = T),]
exp3<-exp2[!duplicated(exp2[,dim(exp2)[2]]),]#express的最后一列为gene name
exp4<-na.omit(exp3)
rownames(exp4)<-exp4[,dim(exp4)[2]]
exp4<-exp4[,-dim(exp4)[2]]
View(exp4)
#保存bulk数据集
write.csv(exp4,"BC_bulk_data.csv",quote = F)
对count文件归一化，计算FPKM
#标准化
group_information<- data.frame(TCGA_id =colnames(exp4))
table(substring(group_information$TCGA_id,14,15))#用table这个函数统计一下BC样本的分类
sample<-ifelse(substring(group_information$TCGA_id,14,15)=="01","cancer","tissue")#这里的factor是为了dds的需要
group_information$sample<-as.factor(sample)#这里的factor是为了dds的需要
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2",force = TRUE)
library(DESeq2)
mycounts<-data.frame("gene_id"=rownames(exp4))
mycounts<-data.frame(mycounts,exp4)
dds<-DESeqDataSetFromMatrix(countData=mycounts,
colData=group_information,
design=~sample,
tidy=TRUE)
dds <- DESeq(dds)
vsd <- vst(dds, blind = FALSE)
bulk_data<-as.data.frame(assay(vsd))
write.csv(exp4,"BC_bulk_data.csv",quote = F)
整合临床生存数据和bulk数据
image.png
image.png
#统一ID
index<-colnames(bulk_data)
index_1<-substring(index,1,4)
index_2<-substring(index,6,7)
index_3<-substring(index,9,12)
a<-list()
for (i in 1:length(index)) {
a[i]<-paste(index_1[i],index_2[i],index_3[i],sep = "-",collapse = NULL)
}
colnames(bulk_data)<-a
#在表达矩阵中找出那些有临床表型数据的列
b<-bulk_data[,colnames(bulk_data) %in% a]
#在临床表型数据中去除那些没有表达矩阵的列
pheno_3<-pheno_2[,colnames(pheno_2) %in% colnames(b)]
#将表达矩阵和临床表型数据整合到一起，列的顺序自然就一样了，再拆分，列名就是对应的了
c<-merge(b,pheno_3,by = colnames(b),all = T,sort = F)
#读取临床生存数据
bulk_survival<-read.csv("bulk_phenotype.csv")
bulk_data_unique<-bulk_data[,!duplicated(colnames(bulk_data))]#去掉bulk_data中重复的列
#在表达矩阵中找出那些有临床表型数据的列
b<-bulk_data_unique[,colnames(bulk_data_unique) %in% bulk_survival$TCGA_patient_barcode]
#bulk_survival$TCGA_patient_barcode的顺序就是列名的顺序，根据它来调整bulk_survival的顺序。
bulk_survival_sorted<-bulk_survival[match(colnames(bulk_data_unique),bulk_survival$TCGA_patient_barcode),]
