生信代码：差异分析（TCGAbiolinks包）
引言：上一期（这里可到达上一期）我们利用得到的肝癌的数据，进行了预处理，得到了最终的表达矩阵TCGA_LIHC_final.csv，今天我们的主要任务就是进行差异表达分析。此外，还会顺带讲两个进行富集分析和聚类分析的函数。
基因差异表达分析
01
# 首先读入表达矩阵文件
dataFilt_LIHC_final <- read.csv("TCGA_LIHC_final.csv", header = T,check.names = FALSE）
# 定义行名
rownames(dataFilt_LIHC_final) <- dataFilt_LIHC_final[,1]
dataFilt_LIHC_final <- dataFilt_LIHC_final[,-1]# 先看一下矩阵长啥样，心里有个数：每一行是一个基因，每一列是一个样本
View(dataFilt_LIHC_final)
# 定义样本的分组，前340个是肿瘤样本，后50个是正常样本
# 定义肿瘤样本分组
mat1 <- dataFilt_LIHC_final[,1-340]
mat1 <- log(mat1+1)
# 定义正常组织样本分组
mat2 <- dataFilt_LIHC_final[,341-390]
mat2 <- log(mat2+1)# 然后就可以进行差异表达分析啦！
Data_DEGs <- TCGAanalyze_DEA(mat1 = mat1,
mat2 = mat2,
Cond1type = "Tumor",
Cond2type = "Normal",
pipeline="limma",
batch.factors = c("TSS"),
voom = TRUE,
contrast.formula = "Mycontrast=Tumor-Normal")# 差异表达分析结果：
View(Data_DEGs)
接下来我们来详细讲解一下 TCGAanalyze_DEA 函数的参数：
用法：
TCGAanalyze_DEA(mat1, mat2, metadata = TRUE, Cond1type, Cond2type,
pipeline = "edgeR", method = "exactTest", fdr.cut = 1,
logFC.cut = 0, elementsRatio = 30000, batch.factors = NULL,
ClinicalDF = data.frame(), paired = FALSE, log.trans = FALSE,
voom = FALSE, trend = FALSE, MAT = data.frame(),
contrast.formula = "", Condtypes = c())参数：
mat1 | 表达矩阵，行是基因，列是样本（如正常组织的表达矩阵） |
|---|---|
mat2 | 表达矩阵，行是基因，列是样本（如肿瘤组织的表达矩阵） |
metadata | 添加 metadata |
Cond1type | mat1中样品的分组信息（如对照组） |
Cond2type | mat2中样品的分组信息（如病例组） |
pipeline | 指明使用哪个R包："limma" or "edgeR" |
method | 用于pipeline="edgeR"时：2个选项'glmLRT'： 用负二项广义对数线性模型拟合每个基因的read counts'exactTest'： 用genewise精确检验得出两组负二项分布计数的平均值的差异。 |
fdr.cut | 设置p值来筛选差异基因 |
logFC.cut | 设置logFC阈值来筛选差异基因 |
batch.factors | 批处理纠正选项："Plate", "TSS", "Year", "Portion", "Center", "Patients" |
ClinicalDF | GDCquery_clinic()返回的数据框，用于提取年份数据 |
paired | 成对的样本设置为TRUE |
log.trans | 当要转换为log-CPM值时，设置为TRUE |
voom | 当进行voom转化时，设置为TRUE（voom参数详细可见limma包） |
trend | 要进行经验贝叶斯先验趋势时，设置为TRUE |
MAT | 一个基因表达矩阵，若已提供mat1 和 mat2，此参数不用设置 |
contrast.formula | 自主设置系数和对比度 |
Condtypes | MAT中的样本分组 |
富集分析
02
# 设置logFC，挑选表达有差异的基因进行富集分析
Data_DEGs_high_expr <- Data_DEGs[Data_DEGs$logFC >=1,]
Genelist <- rownames(Data_DEGs_high_expr)# 富集分析
ansEA <- TCGAanalyze_EAcomplete(TFname="DEA genes Normal Vs Tumor",
Genelist)# 富集分析的结果
# 富集分析结果可视化
TCGAvisualize_EAbarplot(tf = rownames(ansEA$ResBP),
GOBPTab = ansEA$ResBP,
GOCCTab = ansEA$ResCC,
GOMFTab = ansEA$ResMF,
PathTab = ansEA$ResPat,
nRGTab = Genelist,
nBar = 10, #显示条形图的数量
filename = "TCGAvisualize_EAbarplot_Output.pdf")
TCGAanalyze_EAcomplete 函数的用法：
TCGAanalyze_EAcomplete(TFname, RegulonList)参数：
TFname | 设置基因列表的名称 |
|---|---|
RegulonList | 差异基因列表 |
TCGAvisualize_EAbarplot函数的用法：
TCGAvisualize_EAbarplot(tf, GOMFTab, GOBPTab, GOCCTab, PathTab, nBar,
nRGTab, filename = "TCGAvisualize_EAbarplot_Output.pdf",
text.size = 1, mfrow = c(2, 2), xlim = NULL, color = c("orange",
"cyan", "green", "yellow"))参数：
tf | 基因名列表 |
|---|---|
GOMFTab | TCGAanalyze_EAcomplete 函数的 Molecular Function (MF) 结果 |
GOBPTab | TCGAanalyze_EAcomplete 函数的 Biological Process (BP) 结果 |
GOCCTab | TCGAanalyze_EAcomplete 函数的 Cellular Component (CC) 结果 |
PathTab | TCGAanalyze_EAcomplete 函数的 Pathways EA 结果 |
nBar | 显示的条形图直条的数量，默认是10 |
nRGTab | 经过 logFC 筛选的差异基因列表 |
filename | pdf的名称，若 null，直接返回图片 |
text.size | 文本字号 |
mfrow | 图片的行数和列数，默认"c(2,2)"：2行2列 |
xlim | X轴的范围 |
color | 每个条形图的颜色，默认：c("orange", "cyan","green","yellow") |
聚类分析
02
res.hc <- TCGAanalyze_Clustering(Data_DEGs,
method = "hclust",
methodHC = "ward.D2")# 聚类分析结果：
TCGAanalyze_Clustering 函数的用法：
TCGAanalyze_Clustering(tabDF, method, methodHC = "ward.D2")参数：
tabDF | 一个数据框或数值型矩阵，行是基因，每列是一个样本（来自TCGAPrepare） |
|---|---|
method | 使用的聚类方法，如"hclust"（层次聚类） or "consensus"（一致性聚类） |
methodHC | 使用的层次聚类的方法，使用的方法有"ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (=UPGMC). |
今天又是认真学习的一天呢，我们下一期将要分享的是TCGAbiolinks包的一些可视化功能，如热图、火山图等。哈哈，不见不散哦~
腾讯云开发者
