美文网首页生信分析流程
TCGA基因差异表达分析流程(1)

TCGA基因差异表达分析流程(1)

作者: ZJCLucasC22 | 来源:发表于2019-05-07 22:34 被阅读84次

TCGA简介

1.TCGA数据简介

TCGA全称为The Caner Genome Atlas(肿瘤基因图谱)。该数据库从2005年开始由美国进行主持开发至今。
TCGA里面包含了研究的癌症共33种,其中包括稀有类型癌症10种,统共收录了大约11,000个病人7种类型的数据。
这7种数据主要包括,基因组,转录组,甲基化,miRNA数据,临床信息数据等。全部数据约占2.5PB大小,(≈2560TB≈2621440GB)

1.2TCGA数据管理规范

一般来说,TCGA的数据是以Barcode地形式进行记录的。
Barcode的一般形式如下:
TCGA-02-0001-01B-02-0182-01
那么这个编码是怎么来的呢?通常是经过如下流程:
1.不同的研究机构获得病人的组织,进行编号。(TCGA-02)
2.不同的研究机构把病人的组织统一送到BCR(某个机构)进行统一管理与编号。(TCGA-02-0001)
3.该机构又将病人的组织样本分成好几份,以作备份。(TCGA-02-0001-01B)
4.将备份的组织样本又细分成几小份,用以做不同的检测,比如SNP检测,CNV检测。
(TCGA-02-0001-01B-02)
5.把这些微量样本放在96孔板中进行检测。(TCGA-02-0001-01B-02-0182)
6.最后为不同的检测进行记录。(TCGA-02-0001-01B-02-0182-01)

下面再来具体介绍Barcode的含义,
以TCGA-02-0001-01B-02D-0182-01为例:

TCGA 02 0001 01 B 02 D 0182 01
project TSS Participant Sample Vial Portion Anylate Plate Center
image.png
image.png
Label Identifier for Value
Project Project name TCGA
TSS 组织来源,来自于哪个组织,每个组织都有固定的编号(Tissue source site) 02
Participant 病人编号 0001
Sample 样本类型(一般分为两类,正常NP和瘤TP) 01
Vial 样本序列中的样本次序 C
Portion 100~120毫克样品中的部分顺序 01
Analyte 分析类型 D
Plate 放在序号为0182的96孔板 0182
Center 该96孔板送去某个机构进行检测 01

1.3GDC网站介绍

下载TCGA的对外公开数据一般通过GDC网站进行。
https://portal.gdc.cancer.gov/repository
网站一览:

image.png
下面逐个介绍网站值得注意的部分:
Data Category:
image.png
1.SNV/SNP数据
2.转录谱数据
3.原始数据,(一般拿不到,需要权限)
4.生物样本信息
5.一般来说记录着基因组上的变化
6.临床数据
7.DNA甲基化信息
……
Data Type

这里保存着不同的数据类型,比如切片图片,BAM文件。
Program
image.png
不同机构的研究资料,TCGA只是其中一种。
Project
image.png
不同的癌症种类
https://www.jianshu.com/p/3c0f74e85825
该网站有所有的缩写的中文介绍。
想知道这些简写代表什么癌症,可以点击该网页进行查询。
接下来就不一一介绍了。
接下来要介绍如何选择需要的数据。
如果需要乳腺癌的数据,则选择Project中的TCGA-BRCA:
image.png
返回File中
image.png
选择相应的数据类型,
image.png
经过相应的选择后,就可以下载数据进行分析

2 TCGA 下载

2.1基因和lncRNA表达数据下载

本小节内容主要介绍利用R软件进行下载转录组数据,并从中提取基因和lncRNA表达量。
首先先打开RStudio。

# 安装TCGAbiolinks包
source("https://bioconductor.org/biocLite.R")
biocLite("TCGAbiolinks")
#  加载需要的包
library(SummarizedExperiment)
library(TCGAbiolinks)
#设置下载的地方
work_dir <- "/train/TCGA/lab/Download_Data/Gene_lncRNA"

# 设置程序参数
project <- "TCGA-CHOL"#癌症类型:胆管癌
data_category <- "Transcriptome Profiling"#数据类别为转录组
data_type <- "Gene Expression Quantification"#基因定量表达
workflow_type <- "HTSeq - Counts"
legacy <- FALSE #使用的是hg38,TURE为hg19
image.png
# 设置工作目录
setwd(work_dir)
getwd()
image.png
# 数据下载在本目录的GDC文件下
DataDirectory <- paste0(work_dir,"/GDC/",gsub("-","_",project))
FileNameData <- paste0(DataDirectory, "_","RNAseq_HTSeq_Counts",".rda")
#gsub()函数的作用是将 "-" 替换成 "_" 
#"RNAseq_HTSeq_Counts"是文件名
image.png
# 查询可以下载的数据,相当于在TCGA网站上选取相应的模块进行数据筛选
query <- GDCquery(project = project,
                  data.category = data_category,
                  data.type = data_type, 
                  workflow.type = workflow_type,
                  legacy = legacy)
image.png

成功结果如图,有时候服务器不稳定,需要连接多几次。


image.png
# 该癌症的总样品数量
samplesDown <- getResults(query,cols=c("cases"))
cat("Total sample to download:", length(samplesDown))
image.png
# TP 样品数量
dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown, 
typesample = "TP")
cat("Total TP samples to down:", length(dataSmTP))
5
# NT 样本数量
dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown,
typesample = "NT")
cat("Total NT samples to down:", length(dataSmNT))
image.png
# 下载数据, 数据比较大,耐心等待
GDCdownload(query = query,
            directory = DataDirectory,files.per.chunk=6, 
method='client')

真的等了挺久,大概二十分钟……


image.png
# 保存结果,方便后面使用
data <- GDCprepare(query = query, 
                   save = TRUE, 
                   directory =  DataDirectory,
                   save.filename = FileNameData)
image.png

接下来在电脑端查看下载好的数据


image.png
save.filename = FileNameData
image.png

接下来是基因表达量的提取

# 表达量提取
data_expr <- assay(data)
dim(data_expr)#查看基因
image.png

说明有56830个基因,45个样本

#存放数据到文件"All_HTSeq_Counts"
expr_file <- paste0(DataDirectory, "_","All_HTSeq_Counts",".txt")
write.table(data_expr, file = expr_file, 
sep="\t", row.names =T, quote = F)
image.png

文件内容


image.png
image.png

2.2lncRNA表达量提取

本小节介绍TCGA所采用的GTF注释文件中lncRNA信息的提取,如何从下载的转录组数据中提取lncRNA的表达量
https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/#rna-seq-alignment-command-line-parameters

image.png
该注释(V22版本)文件会告诉我们哪些是基因,哪些是lncRNA,然后从中可以选择提出来。
https://www.gencodegenes.org/human/
网页一览:
从此处可以下载GTF注释文件,当前版本已经更新到第30版
image.png
这个网站介绍了GTF如何将基因区分,主要分为四类:
基因编码蛋白,假基因,Long noncoding,Short noncoding
http://asia.ensembl.org/Help/Faq?id=468
image.png
注意:GTF文件未经过处理的话占用的内存很大,利用已经得到的GTF文件,还有Perl语言脚本处理(本章节不详细介绍),可以将GTF文件处理成更小规模的gene_info.txt,和lncRNA_info.txt。得到这两个文件后才能进行表达量矩阵提取
gene_info.txt内容一览:
image.png
# 提取Gene表达量矩阵
gene_info = read.table("./ref/gene_info.txt",header = 1)
cat("Total Gene annotation:", dim(gene_info)[1])
#选择Gene_info.txt中的Gene_id列
gene_selected <- row.names(data_expr) %in% gene_info$Gene_id
gene_expr <- data_expr[gene_selected,]
cat("Total Gene Selected:", dim(gene_expr)[1])

筛选之后会有些减少,但是影响不大

#文件保存
gene_expr_file <- paste0(DataDirectory, "_","Gene_HTSeq_Counts",".txt")
write.table(gene_expr, file = 
gene_expr_file, sep="\t", row.names =T, quote = F)

lncRNA表达量提取也如出一辙

# 提取lncRNA表达量矩阵
lncRNA_info = read.table("./ref/lncRNA_info.txt",header = 1)
cat("Total LncRNA annotation:", dim(lncRNA_info)[1])
image.png
lncRNA_selected <- row.names(data_expr) %in% lncRNA_info$Gene_id
lncRNA_expr <- data_expr[lncRNA_selected,]
cat("Total LncRNA Selected:", dim(lncRNA_expr)[1])
image.png
# 将 基因ID 转换成基因的名称
name_map <- lncRNA_info[,c('Gene_id','gene_name')]
rownames(name_map) <- name_map$Gene_id
rownames(lncRNA_expr) <- name_map[rownames(lncRNA_expr),]$gene_name

# 由于基因名称有重复,只保留一个基因名称对应的表达量
uniq_gene_name <- unique(rownames(lncRNA_expr))
lncRNA_expr <- lncRNA_expr[uniq_gene_name,]

lnc_expr_file <- paste0(DataDirectory, "_","LncRNA_HTSeq_Counts",
".txt")
write.table(lncRNA_expr, 
file = lnc_expr_file, sep="\t", row.names =T, quote = F)
image.png

2.3miRNA表达数据下载

library(SummarizedExperiment)
library(TCGAbiolinks)
#设置下载的地方
work_dir <- "D:/train/TCGA/lab/Download_Data/miRNA"

project <- "TCGA-CHOL"
data_category <- "Transcriptome Profiling"
data_type <- "miRNA Expression Quantification"#数据类型更改了
workflow_type <- "BCGSC miRNA Profiling"
legacy <- FALSE
# 设置工作目录
setwd(work_dir)
# 下载基因表达量,count数格式的结果
DataDirectory <- paste0(work_dir,"/GDC/",gsub("-","_",project))
FileNameData <- paste0(DataDirectory, "_","miRNA",".rda")
# 查询可以下载的数据,有时候查询的API不稳定,需要多试几次
query <- GDCquery(project = project,
                  data.category = data_category,
                  data.type = data_type, 
                  legacy = legacy)
# 该癌症总样品数量
samplesDown <- getResults(query,cols=c("cases"))
cat("Total sample to down:", length(samplesDown))
image.png
# TP 样品数量
dataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown, typesample = "TP")
cat("Total TP samples to down:", length(dataSmTP))

# NT 样本数量
dataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown, typesample = "NT")
cat("Total NT samples to down:", length(dataSmNT))
image.png
# 下载数据
GDCdownload(query = query,directory = DataDirectory,method='client',files.per.chunk=6)
image.png
# 数据整合,保存结果,方便后面使用
data <- GDCprepare(query = query, 
                   save = TRUE, 
                   directory =  DataDirectory,
                   save.filename = FileNameData)
image.png
# miRNA表达量,count提取
read_countData <-  colnames(data)[grep("read_count", colnames(data))]
selectCol <- c("miRNA_ID",read_countData)
data.miR <- data[,selectCol]
colnames(data.miR) <- gsub("read_count_","", colnames(data.miR))#删除后面的多余地方
miRNA_exp_file <- paste0(DataDirectory, "_","miRNA_Counts",".txt")
write.table(data.miR, file = miRNA_exp_file,sep="\t", row.names =F, quote = F)

2.4临床数据下载

library(SummarizedExperiment)
library(TCGAbiolinks)
#设置下载的地方
work_dir <- "D:/train/TCGA/lab/Download_Data/miRNA"

project <- "TCGA-CHOL"
data_category <- "Clinical"
data_type <- "Clinical Supplement"
legacy <- FALSE

# 设置工作目录
setwd(work_dir)
# 下载基因表达量,count数格式的结果
DataDirectory <- paste0(work_dir,"/GDC/",gsub("-","_",project))
FileNameData <- paste0(DataDirectory, "_","miRNA",".rda")

# 查询可以下载的数据,有时候查询的API不稳定,需要多试几次
query <- GDCquery(project = project,
                  data.category = data_category,
                  data.type = data_type, 
                  legacy = legacy)

# 该癌症总样品数量
samplesDown <- getResults(query,cols=c("cases"))
cat("Total Clinical sample to down:", length(samplesDown))

# 下载数据
GDCdownload(query = query,
            directory = DataDirectory,files.per.chunk=6, method='client')

# 整合下载好的数据
clinical <- GDCprepare_clinic(query, clinical.info = "patient",directory = DataDirectory)

# 将数据保存到文件
clinical_file <- paste0(DataDirectory, "_","clinical",".txt")
write.csv(clinical, file = clinical_file, row.names = F, quote = F)

相关文章

网友评论

    本文标题:TCGA基因差异表达分析流程(1)

    本文链接:https://www.haomeiwen.com/subject/rzbcoqtx.html