TCGA基因差异表达分析流程(2)
2019-05-07 本文已影响20人
ZJCLucasC22
3 基因差异表达分析
3.1 基因差异表达分析
# 安装edgeR 软件包
source("https://bioconductor.org/biocLite.R")
biocLite("edgeR")
install.packages("gplots")
install.packages("ggplot2")
# 装载 edgeR 包
library('edgeR')
library('gplots')
library("ggplot2")
#######################
# 加载数据到R #
#######################
# 设置 FDR 和 fold change 阈值，值不固定
fdr = 0.05
fold_change = 2
# 基因表达矩阵
rawdata_file<-"/train/TCGA/lab/Download_Data/Gene_lncRNA/GDC/TCGA_CHOL_Gene_HTSeq_Counts.txt"
image.png
setwd("D:/train/TCGA/lab/DEG/Gene")
# 读取基因表达的矩阵
rawdata=read.table(rawdata_file,header=TRUE,stringsAsFactors=FALSE,row.names=1)
# 原始数据的基因数和样本量
dim(rawdata)
image.png
45个样本，19981个基因
# 针对基因的表达量进行过滤，过滤标准设置为：至少有25% 的样本，基因的表达量大于 2
quant <- apply(rawdata,1,quantile,0.75)
keep <- which((quant >= 2) == 1)
rawdata <- rawdata[keep,]
dim(rawdata)
#quant>=2表达量大于二则保留。否则去掉
image.png
image.png
image.png
#TCGA.ZH.A8Y4.01A.11R.A41I.07,将第11位转化成因子进行判断
group <- factor(substr(colnames(rawdata),14,14))
summary(group)
image.png
#36癌症 9正常
# 获得基因名称列表
#rownames举例：ENSG00000000003
genes=rownames(rawdata)
#接下来是edgeR包的分析流程
# 制作DEGlist 对象
y <- DGEList(counts=rawdata, genes=genes, group=group)
nrow(y)
image.png
# TMM 标准化
y <- calcNormFactors(y)
# 分布估计，负二项分布
y <- estimateCommonDisp(y, verbose=TRUE)
y <- estimateTagwiseDisp(y)
# 差异分析检验，t检验
et <- exactTest(y)
# 针对FDR, fold change 对基因的显著性表达进行‘多重检验’
et$table$FDR <- p.adjust(et$table$PValue, method="BH")
summary(de<-decideTestsDGE(et,adjust.method="BH",p.value=fdr,lfc=log2(fold_change)))
image.png
# 标识基因的上下调情况
et$table$regulate=as.factor(ifelse(et$table$FDR<fdr&abs(et$table$logFC) >=log2(fold_change), ifelse(et$table$logFC>log2(fold_change),'Up','Down'),'Normal'))
summary(et$table$regulate)
head(et$table)
image.png
#接下来是数据可视化
plotMD(et)
#红色为上调，蓝色为下调
image.png
image.png
# 筛选出差异表的基因
diff_expr_out <- et$table[et$table$regulate !='Normal',]
#################
# 输出结果 #
#################
# 保存结果
write.table(diff_expr_out, file="DE_genes.txt", quote=FALSE, row.names=T, sep="\t")
image.png
# 绘制火山图
# 设置FDR, Fold_Change 的辅助线
fdr_line = -log10(fdr)
fc_line = log2(fold_change)
# 采用ggplot 进行绘图
gp1<-ggplot(et$table,aes(x=logFC,y=-log10(FDR),colour=regulate)) +xlab("log2Fold Change") + ylab("-log10 FDR") +ylim(0,30)
# 基于上下调关系，制定不同的颜色
gp2 <- gp1 + geom_point() +scale_color_manual(values =c("green","black", "red"))
# 增加两条辅助线
gp3 <- gp2 + geom_hline(yintercept=fdr_line,linetype=4)+geom_vline(xintercept=c(-fc_line,fc_line),linetype=4)
gp3
# 保存图片
ggsave("DE_gene_Vocano.pdf",width=16,height=9)
image.png
image.png
# 绘图聚类热图
normData=y$pseudo.counts
heatmapData <- normData[rownames(diff_expr_out),]
# 有些表达量为0， 无法进行log 转换，增加一个小背景: 0.001
heatmapExp=log2(heatmapData+0.001)
heatmapMat=as.matrix(heatmapExp)
#输出的热图的储存路径
pdf(file="DE_gene_heatmap.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
heatmap.2(heatmapMat,col='greenred',trace="none")
dev.off()
image.png
局部放大
image.png
# 样品相关性绘图
pdf(file="gene_cor_cluster.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
# 基于pearson 相关性对样品进行聚类
heatmap(cor(normData))
dev.off()
image.png
3.2 lncRNA基因差异表达分析
# 安装edgeR 软件包
# source("https://bioconductor.org/biocLite.R")
# biocLite("edgeR")
# install.packages("gplots")
# install.packages("ggplot2")
#################################################
# 装载 edgeR 包
library('edgeR')
library('gplots')
library("ggplot2")
#######################
# 加载数据到R #
#######################
# 设置 FDR 和 fold change 阈值
fdr = 0.01
fold_change = 2
# 基因表达矩阵
rawdata_file <- "/Users/LUCAS/Documents/Train/TCGA/lab/Download_Data/Gene_lncRNA/GDC/TCGA_CHOL_LncRNA_HTSeq_Counts.txt"
# 设置工作目录
setwd("D:/train/TCGA/lab/DEG/LncRNA")
# 读取基因表达的矩阵
rawdata=read.table(rawdata_file, header=TRUE, stringsAsFactors=FALSE, row.names=1)
# 原始数据的基因数和样本量
dim(rawdata)
image.png
# 针对基因的表达量进行过滤，过滤标准设置为：至少有25% 的样本，基因的表达量大于 2
quant <- apply(rawdata,1,quantile,0.75)
keep <- which((quant >= 2) == 1)
rawdata <- rawdata[keep,]
dim(rawdata)
image.png
表达量很低,一下少了一半有多
# 基于样本的编号，判定是癌症还是正常样本
group <- factor(substr(colnames(rawdata),14,14))
summary(group)
# 获得基因名称列表
genes=rownames(rawdata)
# 制作DEGlist 对象
y <- DGEList(counts=rawdata, genes=genes, group=group)
nrow(y)
# TMM 标准化
y <- calcNormFactors(y)
# 分布估计
y <- estimateCommonDisp(y, verbose=TRUE)
y <- estimateTagwiseDisp(y)
# 差异分析检验
et <- exactTest(y)
# 针对FDR, fold change 对基因的显著性表达进行多重检验
et$table$FDR <- p.adjust(et$table$PValue, method="BH")
summary(de <- decideTestsDGE(et, adjust.method="BH", p.value=fdr, lfc=log2(fold_change)))
image.png
# 标识基因的上下调情况
et$table$regulate=as.factor(ifelse(et$table$FDR<fdr&abs(et$table$logFC)>=log2(fold_change),ifelse(et$table$logFC>log2(fold_change),'Up','Down'),'Normal'))
summary(et$table$regulate)
image.png
plotMD(et)
image.png
image.png
# 筛选出差异表的基因
diff_expr_out <- et$table[et$table$regulate !='Normal',]
#################
# 输出结果 #
#################
# 保存结果
write.table(diff_expr_out, file="DE_lncRNA.txt", quote=FALSE, row.names=T, sep="\t")
# 绘制火山图
# 设置FDR, Fold_Change 的辅助线
fdr_line = -log10(fdr)
fc_line = log2(fold_change)
# 采用ggplot 进行绘图
gp1 <- ggplot(et$table, aes(x=logFC,y=-log10(FDR), colour=regulate)) +xlab("log2Fold Change") + ylab("-log10 FDR") +ylim(0,30)
# 基于上下调关系，制定不同的颜色
gp2 <- gp1 + geom_point() +scale_color_manual(values =c("green","black", "red"))
# 增加两条辅助线
gp3 <- gp2 + geom_hline(yintercept=fdr_line,linetype=4)+geom_vline(xintercept=c(-fc_line,fc_line),linetype=4)
gp3
# 保存图片
ggsave("DE_lncRNA_Vocano.pdf",width=16,height=9)
image.png
image.png
# 绘图聚类热图
normData=y$pseudo.counts
heatmapData <- normData[rownames(diff_expr_out),]
# 有些表达量为0， 无法进行log 转换，增加一个小背景: 0.001
heatmapExp=log2(heatmapData+0.001)
heatmapMat=as.matrix(heatmapExp)
#输出的热图的储存路径
pdf(file="DE_lncRNA_heatmap.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
heatmap.2(heatmapMat,col='greenred',trace="none")
dev.off()
image.png
# 样品相关性绘图
#画这种图的必要性以后再说……
pdf(file="lncRNA_cor_cluster.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
# 基于pearson 相关性对样品进行聚类
heatmap(cor(normData))
dev.off()
image.png
3.3 miRNA差异表达分析
# 装载 edgeR 包
library('edgeR')
library('gplots')
library("ggplot2")
#######################
# 加载数据到R #
#######################
# 设置脚本输出参数
args <- commandArgs(TRUE)
# 设置 FDR 和 fold change 阈值
fdr = 0.05
fold_change = 2
setwd("D:/train/TCGA/lab/DEG/miRNA")
# 基因表达矩阵
rawdata_file<-"/train/TCGA/lab/Download_Data/Gene_lncRNA/GDC/TCGA_CHOL_miRNA_Counts.txt"
# 读取基因表达的矩阵
rawdata=read.table(rawdata_file, header=TRUE, stringsAsFactors=FALSE, row.names=1)
# 原始数据的基因数和样本量
dim(rawdata)
image.png
# 针对基因的表达量进行过滤，过滤标准设置为：至少有25% 的样本，基因的表达量大于 2
quant <- apply(rawdata,1,quantile,0.75)
keep <- which((quant >= 2) == 1)
rawdata <- rawdata[keep,]
dim(rawdata)
image.png
# 基于样本的编号，判定是癌症还是正常样本
group <- factor(substr(colnames(rawdata),25,25))
#在TCGA_CHOL_miRNA_Counts.txt中，TCGA的Barcode变成了read_count_TCGA-W5-AA34-11A-11R-A41D-13
#所以不是第14位而是第25位
summary(group)
image.png
# 获得基因名称列表
genes=rownames(rawdata)
# 制作DEGlist 对象
y <- DGEList(counts=rawdata, genes=genes, group=group)
nrow(y)
# TMM 标准化
y <- calcNormFactors(y)
# 分布估计
y <- estimateCommonDisp(y, verbose=TRUE)
y <- estimateTagwiseDisp(y)
# 差异分析检验
et <- exactTest(y)
# 针对FDR, fold change 对基因的显著性表达进行多重检验
et$table$FDR <- p.adjust(et$table$PValue, method="BH")
summary(de <- decideTestsDGE(et, adjust.method="BH", p.value=fdr, lfc=log2(fold_change)))
# 标识基因的上下调情况
et$table$regulate = as.factor(ifelse(et$table$FDR < fdr & abs(et$table$logFC) >=log2(fold_change), ifelse(et$table$logFC>log2(fold_change),'Up','Down'),'Normal'))
image.png
summary(et$table$regulate)
image.png
plotMD(et)
print('et:')
et
image.png
# 筛选出差异表的基因
diff_expr_out <- et$table[et$table$regulate !='Normal',]
#################
# 输出结果 #
#################
# 保存结果
write.table(diff_expr_out, file="DE_miRNA.txt", quote=FALSE, row.names=T, sep="\t")
# 绘制火山图
# 设置FDR, Fold_Change 的辅助线
fdr_line = -log10(fdr)
fc_line = log2(fold_change)
# 采用ggplot 进行绘图
gp1 <- ggplot(et$table, aes(x=logFC,y=-log10(FDR), colour=regulate)) +xlab("log2Fold Change") + ylab("-log10 FDR") +ylim(0,30)
# 基于上下调关系，制定不同的颜色
gp2 <- gp1 + geom_point() +scale_color_manual(values =c("green","black", "red"))
# 增加两条辅助线
gp3 <- gp2 + geom_hline(yintercept=fdr_line,linetype=4)+geom_vline(xintercept=c(-fc_line,fc_line),linetype=4)
gp3
# 保存图片
ggsave("DE_miRNA_Vocano.pdf",width=16,height=9)
image.png
# 绘图聚类热图
normData=y$pseudo.counts
heatmapData <- normData[rownames(diff_expr_out),]
# 有些表达量为0， 无法进行log 转换，增加一个小背景: 0.001
heatmapExp=log2(heatmapData+0.001)
heatmapMat=as.matrix(heatmapExp)
#输出的热图的储存路径
pdf(file="DE_miRNA_heatmap.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
heatmap.2(heatmapMat,col='greenred',trace="none")
dev.off()
image.png
# 样品相关性绘图
#至于为什么要画着张图，以后再更新……
pdf(file="miRNA_cor_cluster.pdf",width=60,height=90)
par(oma=c(10,3,3,7))
# 基于pearson 相关性对样品进行聚类
heatmap(cor(normData))
dev.off()
image.png
