前言交代1、学习参考之前参加了生信技能树花花老师的TCGA数据挖掘试讲课,收获很多,最近整理一下上课笔记,同时参考了老师的简书相关教程。生信入门的朋友也可微信加入生信星球公众号,个人觉得很好的一个学习平台。一起来学习吧~ 2、笔记内容(1)在TCGA公共数据库,下载某一类癌症的某一种RNA的count数据,以及相应的病人临床信息;
3、学习目标TCGA的大致流程;R操作..... 第一步:数据下载与整理1、下载1.1 官网下载“脚本”文件
mkdir mirna clinical head gdc_manifest.2020-05-02.clinical.txt ./gdc-client download -m gdc_manifest.2020-05-02.mirna.txt -d mirna/ # -d表示储存路径 ./gdc-client download -m gdc_manifest.2020-05-02.clinical.txt -d clinical/ # 切换到R length(dir("./clinical/")) #统计子目录数 length(dir("./mirna/"))
2、整理2.1 整理临床信息library(XML) #clinical临床信息为xml格式,我们目的转换成df result <- xmlParse("./clinical/0f133012-23ef-4237-acfd-47b132b99775/nationwidechildrens.org_clinical.TCGA-W5-AA2Q.xml") rootnode <- xmlRoot(result) rootsize <- xmlSize(rootnode) #结果显示有两个节点 print(rootnode[2]) #病人信息在第二部分 xmldataframe <- xmlToDataFrame(rootnode[2]) #转换成df,发现是长长的一行 head(t(xmlToDataFrame(rootnode[2]))) 临床信息很多 xmls = dir("clinical/",pattern = "*.xml$",recursive = T) #目的是取出所有的xml文件,pattern = "*.xml$"用正则表达式匹配;recursive = T表示递归模式。 head(xmls) #是带路径的 head(xmls) td = function(x){ result <- xmlParse(file.path("clinical/",x)) #补成全路径 rootnode <- xmlRoot(result) xmldataframe <- xmlToDataFrame(rootnode[2])[c( 'bcr_patient_barcode', 'vital_status', 'days_to_death', 'days_to_last_followup', 'race_list', 'days_to_birth', 'gender' , 'stage_event' )] return(t(xmldataframe)) #返回每个病人的信息df } #编写小函数,批量处理 cl = lapply(xmls,td) # xmls即为td函数的x,lapply逐个处理,结果储存在cl列表中 do.call(cbind,cl) cl_df <- t(do.call(cbind,cl)) #将所有病人信息合并 cl_df[1:3,1:3] class(cl_df) #为矩阵 clinical = data.frame(cl_df) rownames(clinical) <- stringr::str_to_upper(clinical$bcr_patient_barcode) clinical <- clinical[,-1] dim(clinical) clinical[1:4,1:4] 最终的clinical数据 2.2 整理表达矩阵length(dir("./mirna/")) #157 x = read.table("mirna/01b6ce76-61a8-4bc4-b0d7-f69d4ce187d6/3fa67c62-6856-42fe-a9f1-fa6831f42b53.mirbase21.mirnas.quantification.txt" ,header = TRUE) head(x) 原始每个样本表达矩阵信息 mis = dir("mirna/",pattern = "*tification.txt$",recursive = T) mis length(mis) #同样先取名 ex = function(x){ result <- read.table(file.path("mirna/",x),sep = "\t",header = TRUE)[,1:2] return(result) } #批量读取函数 mi = lapply(mis,ex) #结果是第一列为基因名,第二列为count数的两列数据 mi_df <- t(do.call(cbind,mi)) #合并并转置,因为45个病人所以90行 dim(mi_df) tmp = mi_df[1:4,1:4] tmp tmp colnames(mi_df) <- mi_df[1,] #添加列名 mi_df[1:4,1:5] #奇数列(基因名)是重复多余的,只保留偶数列 mi_df <- mi_df[seq(2,nrow(mi_df),2),] dim(mi_df) mi_df[1:4,1:4] #转为数值型 mi_df <- apply(mi_df, 2, as.numeric) mi_df[1:4,1:4] mi_df
#接下来利用json文件给表达矩阵加上行名---- dim(mi_df) meta <- jsonlite::fromJSON("metadata.cart.2020-05-02.157.json") #观察到associated_entities列中每一格都是一个df,里面存放着病人的各种相关ID; #其中entity_submitter_id即我们想要的,meta有157行,就对上了 meta$associated_entities[[1]] meta里的id entity <- meta$associated_entities #取第四列,为列表,包含4个元素,每个元素为一个df class(entity[157]) entity[[157]][4] class(entity) jh = function(x){ as.character(x[4]) } #取df中的entity_submitter_id,即我们最终想给矩阵加上的样本名 jh(entity[[1]]) ID = sapply(entity,jh) #取得了所有的病人ID head(ID) image.png options(stringsAsFactors = F) file2id = data.frame(file_name = meta$file_name, #meta中file_name列即为我们下的mrna的文件名字 ID = ID) #file2id表格157行2列,即储存着我们想要的文件名与ID的对应关系 head(mis) #mis文件是我们下载mrna的gz文件的递归路径,注意是有顺序的 head(mis) #注意mis的样本顺序与之前表达矩阵样本顺序相同 mis2 = stringr::str_split(mis,"/",simplify = T)[,2] #取文件名 mis2[1] %in% file2id$file_name #序列匹配与否 #接下来就要把file2id列表按mis2排序,再把病人ID加到表达矩阵上 #match(A,B) head(match(mis2,file2id$file_name)) #排序可以验证一下,match返回值的第一个是38,意思mis2的第一个元素是file2id$file_name的第38个元素。 mis2[1] == file2id$file_name[38] row_tcga = file2id[match(mis2,file2id$file_name),] #调整file2id行顺序以适应我们下载的miran数据 #如上match(mis2,file2id$file_name)第一个值为4;就把第四行放到第一行 rownames(mi_df) = row_tcga$ID #最后一步 #给mi_df添加行名 mi_df[1:4,1:4] mi_df[1:4,1:4] expr = t(mi_df) dim(expr) expr = expr[apply(expr, 1, function(x) { sum(x > 1) > 9 #过滤掉低表达基因 }), ] dim(expr) group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,'tumor','normal') # 01转为tumor,大于10转为normal table(group_list) #只有tumor样本 # group_list # tumor # 157 group_list <- factor(group_list,levels = c("normal","tumor")) #分组因子
save(expr,clinical,group_list,file = "gdc.Rdata")
|
|
来自: 新用户3677sdB0 > 《数据库》