-
数据下载
TCGA官网
数据下载步骤
-
从网页选择数据，下载manifest文件
- 下载表达数据expdata_manifest文件步骤：
-
Files ——> transcriptome profiling(Data Category) ——>Gene Expression Quantification (Data Type) ——>RNAseq(Experimental Strategy) ——> HTSeq-Counts (Workflow Type)
具体 -
Cases ——> TCGA (program) ——> TCGA-CHOL (project)
cases
TCGA-CHOL
添加到cart
点击cart
下载manifest文件方式1
下载manifest文件方式2
manifest文件改名
下载临床clinical_manifest信息
流程
流程
致此下载好了clinical_manifest并改名
cl.manifest
下载好gdc_client.exe，放在工作目录中
打开R软件terminal终端，并设置好工作目录
image.png
#先在工作目录中，新建两个文件夹，分别存放表达数据和临床信息
options(stringsAsFactors = F)
library(stringr)
proj="TCGA-CHOL"
if(!dir.exists("clinical"))dir.create("clinical")
if(!dir.exists("expdata"))dir.create("expdata")
dir()
#下面两行命令在terminal完成
./gdc-client download -m gdc_manifest_cl.2021-06-30.txt -d clinical
./gdc-client download -m gdc_manifest_expdata.2021-06-30.txt -d expdata
##
length(dir("./clinical/"))
length(dir("./expdata/"))
-
信息整理
#3 临床信息整理
library(XML)
result <- xmlParse("./clinical/142aea0e-7a7b-4ac4-9dbb-0f62e2379599/nationwidechildrens.org_clinical.TCGA-W5-AA2O.xml")
rootnode <- xmlRoot(result)
xmlSize(rootnode)#有两个节点的文件
#第1段是头文件信息
print(rootnode[1])
#print(rootnode[2])
xmldataframe <- xmlToDataFrame(rootnode[2])
head(t(xmlToDataFrame(rootnode[2])))#xml导成dataframe
#批量读取所有
xmls = dir("clinical/",pattern = "*.xml$",recursive = T)#recursive = T递归读取
cl = list()
for(i in 1:length(xmls)){
result <- xmlParse(paste0("clinical/",xmls[[i]]))
rootnode <- xmlRoot(result)
cl[[i]] = xmlToDataFrame(rootnode[2])
}
#拼接起来
clinical <- do.call(rbind,cl)
clinical[1:3,1:3]
# 4.整理表达矩阵
#探索数据：先任选两个counts文件读取，并观察geneid的顺序是否一致。
options(stringsAsFactors = F)
x = read.table("expdata/03aee74e-4e37-4a58-a720-c90e807d2f40/be5bc6a0-9720-47ac-953e-fa8d0c32cd82.htseq.counts.gz",
row.names = 1,sep = "\t")
x2 = read.table("expdata/10d08172-48d2-49e7-b760-721163492cc1/c1071bcd-5a0c-4e09-a578-fc4b6dbe26ad.htseq.counts.gz",
row.names = 1,sep = "\t")
#两个样本counts的geneid顺序是一致的
identical(rownames(x),rownames(x2))
#由此可知，他们的geneid顺序是一致的，可以直接cbind，不会导致顺序错乱。
#批量读取所有的counts.gz文件。
count_files = dir("expdata/",pattern = "*.htseq.counts.gz$",recursive = T)
exp = list()
for(i in 1:length(count_files)){
exp[[i]] <- read.table(paste0("expdata/",count_files[[i]]),row.names = 1,sep = "\t")
}
exp <- do.call(cbind,exp)
dim(exp)
exp[1:4,1:4]
#发现问题：这样产生出来的表达矩阵没有列名。
#解决办法：找到一个文件名与样本ID一一对应的文件。cart-json文件。
meta <- jsonlite::fromJSON("metadata.cart.2021-06-30.json")
colnames(meta)
ids <- meta$associated_entities;class(ids)
ids[[1]]
ids[[1]][,1]
#可以看到，meta$associated_entities是个列表，这个列表里包含数据框，数据框的第一列内容就是tcga样本id。
#注意，换了数据需要自己探索存放在哪一列。不一定是完全一样的，需要确认清楚。
ID = sapply(ids,function(x){x[,1]})#对里面的数据取第一列，将ID提取出来
file2id = data.frame(file_name = meta$file_name,
ID = ID)
#文件名与TCGA样本ID的对应关系已经得到，接下来是将其添加到表达矩阵中，成为行名。需要找到读取文件的顺序，一一对应修改。
head(file2id$file_name,2)
head(count_files,2)
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
table(count_files2 %in% file2id$file_name)
#count_files2的顺序就是列名的顺序，根据它来调整file2id的顺序。此处需要再次理解一下match函数。
file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
exp[1:4,1:4]
#表达矩阵整理完成,需要过滤一下那些在很多样本里表达量都为0的基因。过滤标准不唯一
dim(exp)
#exp = exp[rowSums(exp)>0,]
#至少在9个样本里的表达大于1
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
dim(exp)
exp[1:4,1:4]
exp = as.matrix(exp)
# 分组信息
#根据样本ID的第14-15位，给样本分组（tumor和normal）
table(str_sub(colnames(exp),14,15))
Group = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
Group = factor(Group,levels = c("normal","tumor"))
table(Group)
save(exp,clinical,Group,proj,file = paste0(proj,"_gdc.Rdata"))
-
差异分析
见链接差异分析
4、id转换
#TCGA的RNA-seq数据使用的geneid是ensembl id,两个常见的需求：
#1.差异分析结果中每个ensembl id对应的symbol和类型(mRNA/lncRNA或其它)
#2.将行名从ensembl id 转换为symbol
### 02.思路
#1.找到TCGA数据对应的参考基因组注释版本。
#2.下载该版本的参考基因组注释文件，提取ensembl id 与symbol的对应关系及每个基因的gene type信息。
#3.可以将symbol和gene type 用merge添加到差异分析结果中，也可以在差异分析前先转换矩阵的行名。
options(stringsAsFactors = F)
if(!file.exists("gtf_gene.Rdata")){
#step1:读取并探索gtf文件----
#BiocManager::install("rtracklayer")
library(rtracklayer)
gtf = rtracklayer::import("gencode.v22.annotation.gtf")
class(gtf)
gtf = as.data.frame(gtf);dim(gtf)
colnames(gtf)
table(gtf$type)
#step2:先筛选出gene对应的行
gtf_gene = gtf[gtf$type=="gene",]
save(gtf_gene,file = "gtf_gene.Rdata")
}
load("gtf_gene.Rdata")
load("TCGA-CHOL_DEG.Rdata")#如果自己准备话，就将数据做完差异分析后再进行ID转换
deg = DESeq2_DEG
table(rownames(deg) %in% gtf_gene$gene_id)
an = gtf_gene[,c("gene_name","gene_id","gene_type")]
deg = merge(deg,an,by.x = "row.names",by.y = "gene_id")
# mRNA和lncRNA总共有多少个？
lnc = c("3prime_overlapping_ncRNA",
"antisense",
"bidirectional_promoter_lncRNA",
"lincRNA",
"macro_lncRNA",
"non_coding",
"processed_transcript",
"sense_intronic" ,
"sense_overlapping")
k1 = gtf_gene$gene_type %in% lnc;table(k1)
k2 = gtf_gene$gene_type == "protein_coding";table(k2)
# deg中有多少mRNA和lncRNA?
k3 = deg$gene_type %in% lnc;table(k3)
k4 = deg$gene_type =="protein_coding";table(k4)
# 差异的 mRNA和lncRNA 各有多少
k5 = deg$change !="NOT"
table(k3&k5)
table(k4&k5)
### 表达矩阵的行名id转换
rm(list = ls())
load("TCGA-CHOL_gdc.Rdata")
# 一个简化函数，目前仅转换mRNA和lncRNA，其他类型未考虑。
library(tinyarray)
result = trans_exp(exp)
load("gtf_gene.Rdata")
an = gtf_gene[,c("gene_name","gene_id","gene_type")]
exp = exp[rownames(exp) %in% an$gene_id,]
an = an[match(rownames(exp),an$gene_id),]
identical(an$gene_id,rownames(exp))
k = !duplicated(an$gene_name);table(k)
an = an[k,]
exp = exp[k,]
rownames(exp) = an$gene_name
# 最终得到的结果
exp[1:2,1:2]
save(exp,file = paste0(proj,"_symbol_exp.Rdata"))
