Differential expression gene analysis based on DEFB1 expression levels

Preparing dataset of Esophageal Carcinoma (ESCA)

library(DESeq2)
library(stringr)
load("./dataset/expr_TCGA-ESCA.rda")
counts <- assay(`expr_TCGA-ESCA`)

library(gdata)
sample_group <- read.xls("sample.xlsx",sheet=1,stringsAsFactors=F)
counts_subset  <- counts[,str_sub(colnames(counts),1,15) %in% sample_group$sample]
colnames(counts_subset) <- str_sub(colnames(counts_subset),1,15)

rownames(sample_group) <- sample_group$sample
sample_group <- sample_group[,-c(1,2)]
sample_group$group <- factor(sample_group$group,levels=c("L","H"))
counts_subset <- counts_subset[,rownames(sample_group)]
dds <- DESeqDataSetFromMatrix(countData=counts_subset,colData=sample_group,design=~group)
### dds <- DESeq(dds)
### save(dds,file="./output/dds.rda")
### save.image("./output/tmp190307.RData")
load("./output/dds.rda")
res <- results(dds,contrast=c("group","H","L"))
resOrdered <- res[order(res$pvalue),]
write.csv(as.data.frame(resOrdered), file="./output/Table 1 Group_HvsL_results.csv")

library(org.Hs.eg.db)
ensg2eg <- select(org.Hs.eg.db, keys=rownames(res),columns=c("ENSEMBL","ENTREZID"),keytype="ENSEMBL")
res_eg <- res[ensg2eg$ENSEMBL[!is.na(ensg2eg$ENTREZID)],]
rownames(res_eg) <-  ensg2eg$ENTREZID[match(rownames(res_eg),ensg2eg$ENSEMBL)]
dds_eg <- dds[ensg2eg$ENSEMBL[!is.na(ensg2eg$ENTREZID)],]
rownames(dds_eg) <-  ensg2eg$ENTREZID[match(rownames(dds_eg),ensg2eg$ENSEMBL)]

if(FALSE){
library(ReportingTools)
library(biomaRt)
desReport <- HTMLReport(shortName = 'Differential expression genes analysis', title ='Differential expression genes analysis result', reportDirectory = "./output")
publish(dds_eg, desReport, pvalueCutoff=0.05, annotation.db="org.Hs.eg.db", factor = colData(dds_eg)$group, reportDirectory ="./output")
finish(des2Report)
library(pcaExplorer)
pcaExplorer(dds = dds)
}

deg_sig <- subset(res_eg,padj<0.05 & abs(log2FoldChange) > 1,)
write.csv(as.data.frame(deg_sig), file="./output/Table 1b Group_HvsL_results_padj0.05fc2.csv")
deg_sig_list <- 2^(deg_sig$log2FoldChange)
names(deg_sig_list) <- rownames(deg_sig)
deg_sig_list <- sort(deg_sig_list,decreasing=TRUE)

library(clusterProfiler)
ego_bp <- enrichGO(gene=rownames(deg_sig),OrgDb=org.Hs.eg.db,keyType="ENTREZID",ont="BP",pAdjustMethod="BH",pvalueCutoff=0.01,qvalueCutoff=0.05)
ego_bp <- setReadable(ego_bp,OrgDb=org.Hs.eg.db)
write.csv(ego_bp,file="./output/Table 2a ego_bp.csv")
### ego_bp_simp <- simplify(ego_bp)
### write.csv(ego_bp_simp,file="./output/Table 2b ego_bp_simp.csv")

ego_mf <- enrichGO(gene=rownames(deg_sig),OrgDb=org.Hs.eg.db,keyType="ENTREZID",ont="MF",pAdjustMethod="BH",pvalueCutoff=0.01,qvalueCutoff=0.05)
ego_mf <- setReadable(ego_mf,OrgDb=org.Hs.eg.db)
write.csv(ego_mf,file="./output/Table 3a ego_mf.csv")
### ego_mf_simp <- simplify(ego_mf)
### write.csv(ego_mf_simp,file="./output/Table 3b ego_mf_simp.csv")

ego_cc <- enrichGO(gene=rownames(deg_sig),OrgDb=org.Hs.eg.db,keyType="ENTREZID",ont="CC",pAdjustMethod="BH",pvalueCutoff=0.01,qvalueCutoff=0.05)
ego_cc <- setReadable(ego_cc,OrgDb=org.Hs.eg.db)
write.csv(ego_cc,file="./output/Table 4a ego_cc.csv")
### ego_cc_simp <- simplify(ego_cc)
### write.csv(ego_cc_simp,file="./output/Table 4b ego_cc_simp.csv")

ekegg <- enrichKEGG(gene=rownames(deg_sig),organism="hsa", pvalueCutoff=0.05)
ekegg <- setReadable(ekegg,OrgDb=org.Hs.eg.db,keytype="ENTREZID")
write.csv(as.data.frame(ekegg),file="./output/Table 5 ekegg.csv")


pdf(file="./output/enrichment_dotplot.pdf")
dotplot(ego_bp)
dotplot(ego_mf)
dotplot(ego_cc)
### dotplot(ego_bp_simp)
### dotplot(ego_mf_simp)
### dotplot(ego_cc_simp)
dotplot(ekegg)
dev.off()

pdf(file="./output/enrichment_netplot.pdf")
cnetplot(ego_bp, categorySize="pvalue",foldChange=deg_sig_list)
cnetplot(ego_mf, categorySize="pvalue",foldChange=deg_sig_list)
cnetplot(ego_cc, categorySize="pvalue",foldChange=deg_sig_list)
### cnetplot(ego_bp_simp,showCategory=10,categorySize="pvalue",foldChange=deg_sig_list)
### cnetplot(ego_mf_simp,showCategory=10,categorySize="pvalue",foldChange=deg_sig_list)
### cnetplot(ego_cc_simp,showCategory=10,categorySize="pvalue",foldChange=deg_sig_list)
cnetplot(ekegg, categorySize="pvalue",foldChange=deg_sig_list)
dev.off()

Differential expression genes analysis using padj<0.05, FC>2 or FC<1/2 (first column is gene id)

load("./output/enrichment_out.rd")
library(clusterProfiler)
library(DT)
deg_sig_df <- as.data.frame(deg_sig)
eg2symbol <- bitr(rownames(deg_sig_df),fromType="ENTREZID",toType="SYMBOL",OrgDb="org.Hs.eg.db")
deg_sig_df <- merge(deg_sig_df,eg2symbol,by.x='row.names',by.y="ENTREZID")
datatable(deg_sig_df[,c(1,8,2:7)],extensions="Buttons", options=list(dom="Bfrtip",buttons=c("csv","excel")))

GO BP enrichment analysis

datatable(as.data.frame(ego_bp)[,-8],extensions="Buttons", options=list(dom="Bfrtip",buttons=c("csv","excel")))

dotplot(ego_bp)

cnetplot(ego_bp, categorySize="pvalue",foldChange=deg_sig_list)

GO MF enrichment analysis

datatable(as.data.frame(ego_mf)[,-8],extensions="Buttons", options=list(dom="Bfrtip",buttons=c("csv","excel")))

dotplot(ego_mf)

cnetplot(ego_mf, categorySize="pvalue",foldChange=deg_sig_list)

GO MF enrichment analysis

datatable(as.data.frame(ego_cc)[,-8],extensions="Buttons", options=list(dom="Bfrtip",buttons=c("csv","excel")))

dotplot(ego_cc)

cnetplot(ego_cc, categorySize="pvalue",foldChange=deg_sig_list)

KEGG enrichment analysis

datatable(as.data.frame(ekegg)[,-8],extensions="Buttons", options=list(dom="Bfrtip",buttons=c("csv","excel")))

dotplot(ekegg)

cnetplot(ekegg, categorySize="pvalue",foldChange=deg_sig_list)