TCGA 差异基因分析—科研工具箱

############################################################
# Producer=Jin Wang (https://www.jingege.wang)
############################################################
# 安装edgeR 软件包
# source(\"https://bioconductor.org/biocLite.R\")
# biocLite(\"edgeR\")
# install.packages(\"gplots\")
# install.packages(\"ggplot2\")
############################################################
# 装载 edgeR 包
library(\'edgeR\')
library(\'gplots\')
library(\"ggplot2\")

#######################
#  加载数据到R        #
#######################

# 设置脚本输出参数
args <- commandArgs(TRUE)
# 工作目录
#working_dir <- args[1]
#rawdata_file <- args[2]
#sample_info_file <- args[3]

# 设置 FDR 和 fold change 阈值
fdr = 0.01
fold_change = 2

# 工作目录
working_dir <- \"/Users/zhangqiuxue/Documents/Train/TCGA/lab/DEG/Gene\"
# 基因表达矩阵
rawdata_file <- \"/Users/zhangqiuxue/Documents/Train/TCGA/lab/Download_Data/Gene_lncRNA/GDC/TCGA_CHOL_Gene_HTSeq_Counts.txt\"

# 设置工作目录
setwd(working_dir)

# 读取基因表达的矩阵
rawdata=read.table(rawdata_file, header=TRUE, stringsAsFactors=FALSE, row.names=1)

# 原始数据的基因数和样本量
dim(rawdata)

# 针对基因的表达量进行过滤，过滤标准设置为：至少有25% 的样本，基因的表达量大于 2
quant <- apply(rawdata,1,quantile,0.75)
keep <- which((quant >= 2) == 1)
rawdata <- rawdata[keep,]
dim(rawdata)


# 样分组信息，可以通过文件传入，也可以采用下面的简单方法，基于barcode 判断是癌症还是正常组织
# sample_info =read.table(sample_info_file, header=TRUE, stringsAsFactors=FALSE, row.names=1)

# 基于样本的编号，判定是癌症还是正常样本
############################################################
# Barcode 实例：TCGA-W5-AA34-11A-11R-A41I-07
############################################################

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)))

# 标识基因的上下调情况
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)


plotMD(et)
print(\'et:\')
et

# 筛选出差异表的基因
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\")

# 绘制火山图
# 设置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)


# 绘图聚类热图
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()

# 样品相关性绘图
pdf(file=\"gene_cor_cluster.pdf\",width=60,height=90)
par(oma=c(10,3,3,7))
# 基于pearson 相关性对样品进行聚类
heatmap(cor(normData))
dev.off()