生信代码：绘制热图和火山图
一、加载预处理文件
setwd("---") #记得设置正确的设置默认路径
library("TCGAbiolinks") #加载TCGAbiolonks包
##1 加载前几期的数据，获取45配对的肿瘤样本和正常组织样本
TCGA_LIHC_data <- read.csv(file = "TCGA_LIHC_final.csv",
header = T,
row.names = 1,
check.names = FALSE #保证列名不发生自动更正)
#获取正常组织与肿瘤组织的barcodes
samplesNT <- TCGAquery_SampleTypes(colnames(TCGA_LIHC_data), typesample = c("NT"))
samplesTP <- TCGAquery_SampleTypes(colnames(TCGA_LIHC_data), typesample = c("TP"))
#获取配对正常组织与肿瘤组织的barcodes
paired <- intersect(substr(samplesNT,1,12),substr(samplesTP,1,12))
length(paired)
#[1] 45
#用前12位barcodes为桥梁，建立肿瘤与正常样本barcodes的数据框
NT <- data.frame(NT1=substr(samplesNT,1,12),NT2=samplesNT)
TP <- data.frame(TP1 =substr(samplesTP,1,12),TP2=samplesTP)
TP_NT <- merge(TP,NT,by.x = "TP1",by.y = "NT1")
head(TP_NT,3)
# TP1 TP2 NT2
# 1 TCGA-BC-A10Q TCGA-BC-A10Q-01A-11R-A131-07 TCGA-BC-A10Q-11A-11R-A131-07
# 2 TCGA-BC-A10R TCGA-BC-A10R-01A-11R-A131-07 TCGA-BC-A10R-11A-11R-A131-07
# 3 TCGA-BC-A10T TCGA-BC-A10T-01A-11R-A131-07 TCGA-BC-A10T-11A-11R-A131-07
#1.1获取配对正常组织的barcodes:
TP <- TP_NT$TP2
#1.2获取配对肿瘤组织的barcodes:
NT <- TP_NT$NT2二、配对肿瘤组织与正常组织数据的下载与预处理
#2 参照前面几期进行数据下载和数据预处理:
queryDown <- GDCquery(project = "TCGA-LIHC",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts",
barcode = c(TP, NT))
GDCdownload(queryDown,method = "api", directory = "GDCdata",
files.per.chunk = 10)
#2.1 将数据准备成R语言可处理的形式
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename =
"LIHC_case2.rda")
#2.2数据预处理：根据样本与样本之间的spearman相关系数去掉离群值
dataPrep2 <- TCGAanalyze_Preprocessing(object= dataPrep1,
cor.cut = 0.6,
datatype = "HTSeq - Counts")
#2.3选择肿瘤纯度大于60%的肿瘤样本（因为这是基于前面两期获得的肿瘤样本barcodes,故这里获得的肿瘤样本的肿瘤纯度均大于60%，可不执行此步骤）
# purityDATA <- TCGAtumor_purity(colnames(dataPrep1), 0, 0, 0, 0, 0.6)
# # filtered 为被过滤的数据， pure_barcodes是我们要的肿瘤数据
# Purity.LIHC<-purityDATA$pure_barcodes
# normal.LIHC<-purityDATA$filtered
#
# #获取肿瘤纯度大于60%的样本+正常组织样本,共计90个样本
# puried_data <-dataPrep2[,c(Purity.LIHC,normal.LIHC)]
# #有45个肿瘤纯度大于60%的样本
#2.4基因注释
library("SummarizedExperiment")
rowData(dataPrep1)
# DataFrame with 56512 rows and 3 columns
# ensembl_gene_id external_gene_name original_ensembl_gene_id
# <character> <character> <character>
# ENSG00000000003 ENSG00000000003 TSPAN6 ENSG00000000003.13
# ENSG00000000005 ENSG00000000005 TNMD ENSG00000000005.5
rownames(puried_data)<-rowData(dataPrep1)$external_gene_name
#2.5使用EDAseq进行文库大小和GC丰度标准化
dataNorm <- TCGAanalyze_Normalization(tabDF = puried_data,
geneInfo = geneInfo,
method = "gcContent")
#2.6过滤低count的基因,并将结果输出
dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,
method = "quantile",
qnt.cut = 0.25)
write.csv(dataFilt,file = "paired_TCGA_LIHC_final.csv",quote = FALSE)
#保留的是90个样本（前45肿瘤，后45正常组织）三、差异表达分析
##3 差异表达分析:（前）45个肿瘤样本 vs （后）45个正常组织样本
TCGA_LIHC_data <- read.csv(file = "paired_TCGA_LIHC_final.csv",header = T,row.names = 1,check.names = FALSE)
mat1 <- TCGA_LIHC_data[,1:45]
mat2 <- TCGA_LIHC_data[,46:90]
DEG.LIHC.edgeR <- TCGAanalyze_DEA(mat1 = mat1,
mat2 = mat2,
pipeline = "edgeR",
batch.factors = c("TSS"), Cond1type = "tumor",
Cond2type = "normal",
voom = FALSE, ##设置了paired时，会出错(paired = TRUE)，故此处未设置
method = "glmLRT",
contrast.formula = "Mycontrast =tumor -normal",
fdr.cut = 0.01, #设置过滤参数1，保留FDR＜0.01的基因
logFC.cut = 1 #设置过滤参数2，保留logFC＞1的基因)
write.csv(DEG.LIHC.edgeR,file = "paired_DEG_by_edgeR.csv")四、增加不同分组条件下的gene平均表达量
TCGAanalyze_LevelTab（）将差异表达基因在正常和肿瘤组织中的表达量数据添加到差异表达分析结果中的主要用法：
TCGAanalyze_LevelTab(FC_FDR_table_mRNA, typeCond1, typeCond2, TableCond1,
TableCond2, typeOrder = TRUE)参数详解:
主要参数 | 用法 |
|---|---|
FC_FDR_table_mRNA | 通过LogFC绝对值≥1过滤的差异分析结果数据 |
typeCond1 | 条件1的分类标签，如对照组 |
typeCond2 | 条件2的分类标签，如试验组 |
TableCond1 | 条件1对应的表达矩阵，行代表样本名，列代表基因名 |
TableCond2 | 条件2对应的表达矩阵，行代表样本名，列代表基因名 |
typeOrder | typeOrder |
R中具体示例：
#4.1 TCGAquery_SampleTypes（）用于获取特定组织对应的barcodes,如肿瘤组织（TP）、正常组织（NT）
#设置typesample =NT,获取正常组织对应的barcodes
samplesNT <- TCGAquery_SampleTypes(colnames(dataFilt.LIHC.final), typesample = c("NT"))
#设置typesample=TP,获取肿瘤组织对应的barcodes
samplesTP <- TCGAquery_SampleTypes(colnames(dataFilt.LIHC.final), typesample = c("TP"))
dataDEGsFilt <- DEG.LIHC.edgeR[abs(DEG.LIHC.edgeR$logFC) >= 1,]
str(dataDEGsFilt)
# 'data.frame': 2647 obs. of 5 variables:
#获取肿瘤和正常组织对应的表达矩阵
dataTP <- dataFilt.LIHC.final[,samplesTP]
dataTN <- dataFilt.LIHC.final[,samplesNT]
#4.2 将需要的参数传入dataDEGsFiltLevel （）函数中
dataDEGsFiltLevel <- TCGAanalyze_LevelTab(FC_FDR_table_mRNA = dataDEGsFilt,
typeCond1 = "Normal",
typeCond2 = "Tumor",
TableCond1 = dataTN,
TableCond2 = dataTP)
#4.3 查看结果
head(dataDEGsFiltLevel,2)
# mRNA logFC FDR Normal Tumor Delta
#HP HP 1.184681 7.568768e-06 1087091.5 490505.1 1287856.1
#APOA2 APOA2 -1.155957 6.493243e-06 308058.6 707503.1 356102.5差异表达分析的结果如下：
添加两种条件下基因的平均表达水平后的结果如下:
五、PCA主成分分析
TCGAvisualize_PCA()实现主成分分析的主要用法:
TCGAvisualize_PCA(dataFilt, dataDEGsFiltLevel, ntopgenes, group1, group2)参数详解：
主要参数 | 用法 |
|---|---|
dataFilt | TCGAanalyzeFiltering（）过滤预处理的数据，行代表样本，列代表基因 |
dataDEGsFiltLevel | TCGAanalyzeLevelTab（）输出的结果，具体内容可参见上方的输出结果截图 |
ntopgenes | 在PCA中绘制的差异基因数目，如200 |
group1 | 条件1对应的样本barcodes列表 |
group2 | 条件2对应的样本barcodes列表 |
R中具体示例：
#由于在TCGAanalyze_LevelTab（）中，我们已经得到了一些参数，故可将参数直接带入主成分分析的函数中。
pca <- TCGAvisualize_PCA(dataFilt= dataFilt.LIHC.final,
dataDEGsFiltLevel=dataDEGsFiltLevel,
ntopgenes = 100,
group1=samplesNT,
group2=samplesTP)PCA主成分分析输出结果如下：
六、绘制差异表达基因的热图
TCGAvisualize_Heatmap()绘制热图的主要用法：等号后面对应的为默认参数。
TCGAvisualize_Heatmap(data, col.metadata, row.metadata,
col.colors = NULL, row.colors = NULL, show_column_names = FALSE,
show_row_names = FALSE, cluster_rows = FALSE,
cluster_columns = FALSE, sortCol, extrems = NULL,
rownames.size = 12, title = NULL, color.levels = NULL,
values.label = NULL, filename = "heatmap.pdf", width = 10,
height = 10, type = "expression", scale = "none",
heatmap.legend.color.bar = "continuous")参数详解：
主要参数 | 用法 |
|---|---|
data | 用于绘制热图的举证，如基因表达矩阵或甲基化矩阵 |
col.metadata、row.metadata | 行和（或）列的补充信息，可作为行或列的注释信息 |
col.colors、row.colors | 设置行、列注释注释信息的颜色参数 |
showcolumnnames、showrownames | 是否展示行或（列）的注释信息 |
clusterrows、clustercolumns | 是否根据行或（列）进行聚类信息 |
sortCol | 用于列排序的列名 |
extrems | 颜色的极端值 |
rownames.size | 行名的大小 |
color.levels | 设置不同表达水平的颜色（对应的表达水平分别为：low level, middle level, high level） |
title | 热图的标题 |
filename | 设置保存时的文件名，默认为“heatmap.pdf” |
width、height | 图片的宽和高 |
type | 设置热图中值的颜色，有“expression”和“methylation”，默认为“expression” |
scale | 是否使用z分数标准化，如果想要比较基因间的差异，对样本进行标准化;如果想要比较样本间的差异，对基因进行标准化。可选参数，“row","col","none"。 |
heatmap.legend.color.bar | 设置图例的类型，continuous（连续型）或者disctrete（离散型） |
R中具体示例：
#6.1 获取差异表达基因的表达水平
datDEGs <- dataFilt.LIHC.final[match(rownames(DEG.LIHC.edgeR),rownames(dataFilt.LIHC.final)),]
str(datDEGs)
#'data.frame': 2647 obs. of 90 variables:
#前45是肿瘤样本，后45是正常组织样本
#6.2 根据临床信息构造patient的metadata信息;或（和）基因的相关信息（此处不添加基因的注释信息）
#获取每一个患者barcode（barcode的前12位代表的patient信息，前16位代表的是sample信息）对应的临床信息，但是其barcodes与datDEGs_test顺序不一致
query <- GDCquery(project = "TCGA-LIHC",
data.category = "Clinical",
file.type = "xml",
barcode = substr(colnames(datDEGs),1,12))
GDCdownload(query)
clinical <- GDCprepare_clinic(query,"patient") #48个样本
##根据表达矩阵中的样本barcodes对样本临床信息匹配
datDEGs_test_barcodes <- as.data.frame（substr(colnames(datDEGs),1,12), ncol = 1 )
colnames(datDEGs_test_barcodes) <-"LIHC_patient_barcode"
m <- clinical[match(datDEGs_test_barcodes[,1], clinical[ , 1 ]),]
str(m)
#'data.frame': 90 obs. of 82 variables:
table(duplicated(m))
#FALSE TRUE
#45 45
#使用table(duplicated()查看m矩阵中是否有重复数据。
#这里的重复数据来源（肿瘤组合和癌旁正常组织来源于同一患者）由于使用的是配对正常样本和肿瘤组织，其对应的患者12位barcodes是一致的，在使用TCGAbiolinks包自带的热图绘制函数时会出现样本信息匹配错误，故小编在这里使用pheatmap()绘制热图。
library(pheatmap)
#基本做图结果
pheatmap::pheatmap(datDEGs,scale = "row",show_rownames = F,show_colnames = F)
上图为未设置其他参数时的基本如图结果。为增加图片的信息量，可增加metadata信息(即行注释和列注释信息)，注意pheatmap()的注释信息的列名必须与表达矩阵的行名、列名一致，演示如下：
#增加metadata信息
col.mdat <- data.frame(Sex=m$gender,
status=m$vital_status,
group=c(rep("tumor",45),rep("normal",45)))
rownames(col.mdat) <- colnames(datDEGs)
#保证列注释信息的行名与样本名（对应列）一致
#设置图例的范围
bk <- c(seq(-1,6,by=0.01))
#绘制热图
pheatmap::pheatmap(datDEGs,scale = "row",show_rownames = F,show_colnames = F,
annotation_col = col.mdat3,
border_color=NA,
main = "Heatmap by pheatmap(edgeR)",
filename = "Heatmap_by_pheatmap.pdf",
color =c(colorRampPalette(colors = c("blue","white"))(length(bk)/2),
colorRampPalette(colors = c("white","red"))(length(bk)/2)) #设置图例的颜色,
legend_breaks=seq(-1,6,2),
breaks=bk )
以上为热图的输出结果，我们可以看到按照行（样本）进行聚类，基本上能够把肿瘤组织与正常组织分类开，说明两种组织的基因表达是具有差异的。相反，在不同存活状态和性别中，暂时未能发现于基因差异表达的相关性。
七、绘制差异差异表达分析结果的火山图
TCGAVisualize_volcano（）绘制火山图的主要用法：
TCGAVisualize_volcano(x, y, filename = "volcano.pdf",
ylab = expression(paste(-Log[10], " (FDR corrected -P values)")),
xlab = NULL, title = "Volcano plot", legend = NULL, label = NULL,
xlim = NULL, ylim = NULL, color = c("black", "red", "green"),
names = NULL, names.fill = TRUE, show.names = "significant",
x.cut = 0, y.cut = 0.01, height = 5, width = 10,
highlight = NULL, highlight.color = "orange", names.size = 4,
dpi = 300)参数详解：
主要参数 | 用法 |
|---|---|
x | X轴对应的数据 |
y | y轴对应的数据 |
filename | 设置保存时的文件名，默认为"volcano.pdf" |
ylab、xlab | y轴、x轴的标题 |
title | 图片的标题 |
legend | 设置图例 |
label | 图片中的注释文本。比如"NotSignificant","Hypermethylated in group1"，"Hypomethylated in group1" |
xlim、ylim | x、y轴坐标轴范围 |
color | 设置图片中使用的颜色 |
names | 是否在图中标记具有显著性差异的基因名称 |
names.fill | 是否将具有显著性差异的基因名称写入方框内 |
show.names | 展示哪种基因的名称，可设置的选项："significant"（具有显著性差异差异基因）、"highlighted"（突出显示的基因）或者"both"（以上两种类型的基因名称都显示）。 |
x.cut | x轴的阈值，默认为0.0。如0.2，那么阈值为±0.2；如c（-0.3，-0.4）,则范围为（-0.3，-0.4） |
y.cut | p值的阈值 |
height、width | 图片的高、宽 |
highlight | 需要突出显示的gene或探针列表 |
hight.color | 突出显示的gene的颜色 |
name.size | 设置为“significant”或highlighted”名称的字体的大小 |
R中具体示例：
#为了做图的需要，突出显示FC≥8的gene名称DEG.LIHC.filt<-DEG.LIHC.edgeR[which(abs(DEG.LIHC.edgeR$logFC) >= 8), ]
str(DEG.LIHC.filt)
#'data.frame': 29 obs. of 5 variables:
#共有29个基因满足|logFC|≥8
TCGAVisualize_volcano(DEG.LIHC.edgeR$logFC,
DEG.LIHC.edgeR$FDR,
filename = "TumorvsNormal_FC8.edgeR.pdf",
xlab = "logFC",
names = rownames(DEG.LIHC.edgeR),
show.names = "highlighted",
x.cut = 1,
y.cut = 0.01,
highlight = rownames(DEG.LIHC.edgeR)[which(abs(DEG.LIHC.edgeR$logFC) >= 8)],
highlight.color = "orange",
title = "volcano plot by edgeR")
title = "volcano plot by limma")如下为此次分析的火山图结果，通过查看图片，可以发现一些基因在肿瘤组织中表达量升高较高，而一些基因在肿瘤组织的表达量低于正常组织中，具体它有什么含义，就需要查阅文献明确。
九、结语
今天的热图和火山图就暂告一段落。但在实际过程中应该结合自己的数据，调整一些参数和分组，以得出更有意义的结论，为科研助力......接下来我们将使用TCGAbiolinks包继续演示TCGA数据中甲基化分析，我们一起努力哦~~~
免责声明：本文转载仅仅是出于传播信息的需要，并不意味着代表本公众号观点或证实其内容的真实性；如其他媒体、网站或个人从本公众号转载使用，须保留本公众号注明的“来源”，并自负版权等法律责任；作者如果不希望被转载，请与我们接洽。版权属于原作者
腾讯云开发者
