Created
March 18, 2015 22:38
-
-
Save cplaisier/47a10580ce6dd1a88e9e to your computer and use it in GitHub Desktop.
Using DESeq2
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Load counts | |
| d0 = read.csv('gexp_counts.csv',header=T,row.names=1) | |
| d0 = d0[which(apply(d0,1,sum)!=0),] | |
| library(DESeq2) | |
| conds = data.frame(subset=factor(c(rep('Cntl_CD44Low',4), rep('Cntl_Ag85B_Plus_PD1_Plus',4), rep('Cntl_ESAT6_Plus_PD1_Plus',4), rep('Stim_Ag85B_Plus_PD1_Plus',4), rep('Stim_ESAT6_Plus_PD1_Plus',4)))) | |
| rownames(conds) = colnames(d0) | |
| # Make a map of Entrez IDs to UCSC transcript IDs | |
| d1 = read.csv('entrez2ucsc.csv',header=F) | |
| entrez2ucsc = list() | |
| for(j in 1:length(d1[,1])) { | |
| tmp = intersect(strsplit(as.character(d1[j,2]),split=';')[[1]],rownames(d0)) | |
| if (length(tmp)>1) { | |
| entrez2ucsc[[as.character(j)]] = tmp | |
| } | |
| } | |
| # Get rid of redundant transcripts | |
| # 1. Correlate transcripts that map to a single Entrez ID | |
| # 2. Merge them if they are correlated R >= 0.8 | |
| c1 = c() | |
| ucscNames = rownames(d0) | |
| merge = 0 | |
| mergeGenes = 0 | |
| noMerge = 0 | |
| noMergeGenes = 0 | |
| for(entrez in names(entrez2ucsc)) { | |
| n1 = entrez2ucsc[[entrez]] | |
| tmp = cor(t(d0[n1,])) | |
| c1 = c(c1,tmp[upper.tri(tmp)]) | |
| if(min(tmp)>=0.8) { | |
| min1 = apply(d0[n1,],1,min) | |
| keeper = names(which(min1==min(min1))) | |
| n2 = n1[which(n1!=keeper)] | |
| ucscNames = ucscNames[which(!(ucscNames %in% n2))] | |
| merge = merge + 1 | |
| } else { | |
| noMerge = noMerge + 1 | |
| print(n1) | |
| print(tmp[upper.tri(tmp)]) | |
| # Iterate through matrix | |
| for(i in n1) { | |
| for(j in n1) { | |
| if(!(i==j) && tmp[i,j]>0.8) { | |
| min1 = apply(d0[n1,],1,min) | |
| keeper = names(which(min1==min(min1))) | |
| n2 = n1[which(n1!=keeper)] | |
| ucscNames = ucscNames[which(!(ucscNames %in% n2))] | |
| } | |
| } | |
| } | |
| } | |
| } | |
| # Get rid of transcripts that are redundant | |
| d0 = d0[ucscNames,] | |
| dds = DESeqDataSetFromMatrix(countData = d0, colData = conds, design= ~ subset) | |
| # Load conversion dictionaries | |
| d2 = read.table('knownToLocusLink.txt',sep='\t') | |
| ucsc2entrez = d2[,2] | |
| names(ucsc2entrez) = d2[,1] | |
| d3 = read.delim('Mus_musculus.gene_info',sep='\t',skip=1) | |
| entrez2symbol = d3[,3] | |
| names(entrez2symbol) = as.character(d3[,2]) | |
| # Cluster samples | |
| library(pheatmap) | |
| #rld = rlog(dds,blind=F) | |
| rld = rlog(dds,blind=T) | |
| sd1 = dist(t(assay(rld))) | |
| sdm1 = as.matrix(sd1) | |
| diag(sdm1) = NA | |
| pheatmap(sdm1,color=colorRampPalette(c('blue4','white','gold'))(32)) # Looks much better | |
| c1 = cor(assay(rld),method='spearman') | |
| diag(c1) = NA | |
| #pheatmap(c1,cluster_rows=F,cluster_cols=F,color=colorRampPalette(c('blue4','white','gold'))(32)) | |
| pdf('clusteringOfSamples_rlog_counts.pdf') | |
| pheatmap(c1,color=colorRampPalette(c('blue4','white','gold'))(32),clustering_method='single') # Looks much better | |
| dev.off() | |
| pdf('pca_rlog_counts.pdf') | |
| plotPCA(rld, intgroup=c('subset')) | |
| dev.off() | |
| # Do differential expression analyses | |
| #cds = newCountDataSet(d0, conds) | |
| #cds = estimateSizeFactors(cds) | |
| #cds = estimateDispersions(cds) | |
| #pdf('dispersions_DESeq_counts.pdf') | |
| #plotDispEsts(cds) | |
| #dev.off() | |
| # Run pipeline | |
| dds = DESeq(dds) | |
| # c(rep('Cntl_CD44Low',4), rep('Cntl_Ag85B_Plus_PD1_Plus',4), rep('Cntl_ESAT6_Plus_PD1_Plus',4), rep('Stim_Ag85B_Plus_PD1_Plus',4), rep('Stim_ESAT6_Plus_PD1_Plus',4)) | |
| # Cntl_Ag85B_Plus_PD1_Plus vs. Cntl_CD44Low | |
| res.cntl_ag85b_plus_pd1_plus.vs.cd44low = results(dds, contrast=c('subset', 'Cntl_Ag85B_Plus_PD1_Plus', 'Cntl_CD44Low')) | |
| write.csv(res.cntl_ag85b_plus_pd1_plus.vs.cd44low, 'res.cntl_ag85b_plus_pd1_plus.vs.cd44low_DESeq2.csv') | |
| # Cntl_ESAT6_Plus_PD1_Plus vs. Cntl_CD44Low | |
| res.cntl_esat6_plus_pd1_plus.vs.cd44low = results(dds, contrast=c('subset', 'Cntl_ESAT6_Plus_PD1_Plus', 'Cntl_CD44Low')) | |
| write.csv(res.cntl_esat6_plus_pd1_plus.vs.cd44low, 'res.cntl_esat6_plus_pd1_plus.vs.cd44low_DESeq2.csv') | |
| # Stim_Ag85B_Plus_PD1_Plus vs. Cntl_CD44Low | |
| res.stim_ag85b_plus_pd1_plus.vs.cd44low = results(dds, contrast=c('subset', 'Stim_Ag85B_Plus_PD1_Plus', 'Cntl_CD44Low')) | |
| write.csv(res.stim_ag85b_plus_pd1_plus.vs.cd44low, 'res.stim_ag85b_plus_pd1_plus.vs.cd44low_DESeq2.csv') | |
| # Stim_ESAT6_Plus_PD1_Plus vs. Cntl_CD44Low | |
| res.stim_esat6_plus_pd1_plus.vs.cd44low = results(dds, contrast=c('subset', 'Stim_ESAT6_Plus_PD1_Plus', 'Cntl_CD44Low')) | |
| write.csv(res.stim_esat6_plus_pd1_plus.vs.cd44low, 'res.stim_esat6_plus_pd1_plus.vs.cd44low_DESeq2.csv') | |
| # Stim_Ag85B_Plus_PD1_Plus vs. Cntl_Ag85B_Plus_PD1_Plus | |
| res.stim_ag85b_plus_pd1_plus.vs.cntl_ag85b_plus_pd1_plus = results(dds, contrast=c('subset', 'Stim_Ag85B_Plus_PD1_Plus', 'Cntl_Ag85B_Plus_PD1_Plus')) | |
| write.csv(res.stim_ag85b_plus_pd1_plus.vs.cntl_ag85b_plus_pd1_plus, 'res.stim_ag85b_plus_pd1_plus.vs.cntl_ag85b_plus_pd1_plus_DESeq2.csv') | |
| # Stim_ESAT6_Plus_PD1_Plus vs. Cntl_CD44Low | |
| res.stim_esat6_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus = results(dds, contrast=c('subset', 'Stim_ESAT6_Plus_PD1_Plus', 'Cntl_ESAT6_Plus_PD1_Plus')) | |
| write.csv(res.stim_esat6_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus, 'res.stim_esat6_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus_DESeq2.csv') | |
| # Cntl_Ag85B_Plus_PD1_Plus vs. Cntl_ESAT6_Plus_PD1_Plus | |
| res.cntl_ag85b_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus = results(dds, contrast=c('subset', 'Cntl_Ag85B_Plus_PD1_Plus', 'Cntl_ESAT6_Plus_PD1_Plus')) | |
| write.csv(res.cntl_ag85b_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus, 'res.cntl_ag85b_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus_DESeq2.csv') | |
| # Stim_Ag85B_Plus_PD1_Plus vs. Stim_ESAT6_Plus_PD1_Plus | |
| res.stim_ag85b_plus_pd1_plus.vs.stim_esat6_plus_pd1_plus = results(dds, contrast=c('subset', 'Stim_Ag85B_Plus_PD1_Plus', 'Stim_ESAT6_Plus_PD1_Plus')) | |
| write.csv(res.stim_ag85b_plus_pd1_plus.vs.stim_esat6_plus_pd1_plus, 'res.stim_ag85b_plus_pd1_plus.vs.stim_esat6_plus_pd1_plus_DESeq2.csv') | |
| # A list of all DE genes | |
| allDEUp = c() # 14,934 | |
| allDEDown = c() # 11,375 | |
| # Plot volcano plots | |
| pdf('volcano_plots_DESeq2_NR.pdf') | |
| tmp = res.cntl_ag85b_plus_pd1_plus.vs.cd44low | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.cntl_ag85b_plus_pd1_plus.vs.cd44low_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.cntl_ag85b_plus_pd1_plus.vs.cd44low_DOWN.csv') | |
| # Add title | |
| title(paste('Control-Ag85B+PD1+ vs. CD44Low\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| tmp = res.cntl_esat6_plus_pd1_plus.vs.cd44low | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.cntl_esat6_plus_pd1_plus.vs.cd44low_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.cntl_esat6_plus_pd1_plus.vs.cd44low_DOWN.csv') | |
| # Add title | |
| title(paste('Control-ESAT6+PD1+ vs. CD44Low\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| tmp = res.stim_ag85b_plus_pd1_plus.vs.cd44low | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.stim_ag85b_plus_pd1_plus.vs.cd44low_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.stim_ag85b_plus_pd1_plus.vs.cd44low_DOWN.csv') | |
| # Add title | |
| title(paste('Peptide-Simulated-Ag85B+PD1+ vs. CD44Low\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| tmp = res.stim_esat6_plus_pd1_plus.vs.cd44low | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.stim_esat6_plus_pd1_plus.vs.cd44low_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.stim_esat6_plus_pd1_plus.vs.cd44low_DOWN.csv') | |
| # Add title | |
| title(paste('Peptide-Stimulated-ESAT6+PD1+ vs. CD44Low\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| tmp = res.stim_ag85b_plus_pd1_plus.vs.cntl_ag85b_plus_pd1_plus | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| ag85b_up = up | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| ag85b_down = down | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.stim_ag85b_plus_pd1_plus.vs.cntl_ag85b_plus_pd1_plus_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.stim_ag85b_plus_pd1_plus.vs.cntl_ag85b_plus_pd1_plus_DOWN.csv') | |
| # Add title | |
| title(paste('Peptide-Simulated-Ag85B+PD1+ vs.\nControl-Ag85B+PD1+\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| tmp = res.stim_esat6_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| esat6_up = up | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| esat6_down = down | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.stim_esat6_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.stim_esat6_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus_DOWN.csv') | |
| # Add title | |
| title(paste('Peptide-Stimulated-ESAT6+PD1+ vs.\nControl-ESAT6+PD1+\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| tmp = res.cntl_ag85b_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| esat6_up = up | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| esat6_down = down | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.cntl_ag85b_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.cntl_ag85b_plus_pd1_plus.vs.cntl_esat6_plus_pd1_plus_DOWN.csv') | |
| # Add title | |
| title(paste('Control-Ag85B+PD1+ vs.\nControl-ESAT6+PD1+\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| tmp = res.stim_ag85b_plus_pd1_plus.vs.stim_esat6_plus_pd1_plus | |
| plot(tmp$log2FoldChange, -log10(tmp$padj),pch=20,col=rgb(0,0,1,0.25),ylab='-log10(Adjusted P-Value)',xlab='log2(Fold-Change)') | |
| # Add markers | |
| pdcd1 = which(rownames(tmp)=='uc007cev.1') | |
| ifng = which(rownames(tmp)=='uc007hdy.1') | |
| points(tmp[pdcd1,]$log2FoldChange, -log10(tmp[pdcd1,]$padj),col=rgb(0,1,0,1),pch=15) | |
| points(tmp[ifng,]$log2FoldChange, -log10(tmp[ifng,]$padj),col=rgb(0.6,0,0.6,1),pch=15) | |
| # Add lines | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| up = intersect(rownames(tmp)[which(tmp$log2FoldChange>=1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| esat6_up = up | |
| down = intersect(rownames(tmp)[which(tmp$log2FoldChange<=-1)],rownames(tmp)[which(tmp$padj<=0.05)]) | |
| esat6_down = down | |
| allDEUp = c(allDEUp, up) | |
| allDEDown = c(allDEDown, down) | |
| # Write out differentially expressed transcripts | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% up),],entrez=as.character(ucsc2entrez[up]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[up])])),'res.stim_ag85b_plus_pd1_plus.vs.stim_esat6_plus_pd1_plus_UP.csv') | |
| write.csv(cbind(data.frame(tmp)[which(rownames(tmp) %in% down),],entrez=as.character(ucsc2entrez[down]),symbol=as.character(entrez2symbol[as.character(ucsc2entrez[down])])),'res.stim_ag85b_plus_pd1_plus.vs.stim_esat6_plus_pd1_plus_DOWN.csv') | |
| # Add title | |
| title(paste('Peptide-Stimulated-Ag85b+PD1+ vs.\nPeptide-Stimulated-ESAT6+PD1+\nDown: ',length(down),', Up: ',length(up))) | |
| # Add legend | |
| legend('top',c('Transcript','Pdcd1','Ifng'),col=c(rgb(0,0,1,0.25),rgb(0,1,0,1),rgb(0.6,0,0.6,1)),bg='white',pch=c(20,15,15)) | |
| dev.off() | |
| library(limma) | |
| geneNames = rownames(d0) | |
| c1 = cbind(Ag85b_Up=geneNames %in% ag85b_up,ESAT6_Up=geneNames %in% esat6_up,Ag85b_Down=geneNames %in% ag85b_down,ESAT6_Down=geneNames %in% esat6_down) | |
| rownames(c1) = geneNames | |
| write.csv(c1,'geneGroupings.csv') | |
| a1 = vennCounts(c1) | |
| pdf('vennDiagram.pdf') | |
| vennDiagram(a1) | |
| dev.off() | |
| # Plot volcano plots | |
| pdf('volcano_plots_DESeq2.pdf') | |
| plot(res.cd44low.vs.klrg1_plus_pe_minus$log2FoldChange, -log10(res.cd44low.vs.klrg1_plus_pe_minus$padj),pch=20,col=rgb(0,0,1,0.25)) | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| plot(res.cd44low.vs.klrg1_plus_pe_plus$log2FoldChange, -log10(res.cd44low.vs.klrg1_plus_pe_plus$padj),pch=20,col=rgb(0,0,1,0.25)) | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| plot(res.cd44low.vs.pd1_plus_pe_minus$log2FoldChange, -log10(res.cd44low.vs.pd1_plus_pe_minus$padj),pch=20,col=rgb(0,0,1,0.25)) | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| plot(res.klrg1_plus_pe_minus.vs.pd1_plus_pe_minus$log2FoldChange, -log10(res.klrg1_plus_pe_minus.vs.pd1_plus_pe_minus$padj),pch=20,col=rgb(0,0,1,0.25)) | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| plot(res.klrg1_plus_pe_plus.vs.pd1_plus_pe_minus$log2FoldChange, -log10(res.klrg1_plus_pe_plus.vs.pd1_plus_pe_minus$padj),pch=20,col=rgb(0,0,1,0.25)) | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| plot(res.klrg1_plus_pe_minus.vs.klrg1_plus_pe_plus$log2FoldChange, -log10(res.klrg1_plus_pe_minus.vs.klrg1_plus_pe_plus$padj),pch=20,col=rgb(0,0,1,0.25)) | |
| abline(h=-log10(0.05),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| abline(v=c(-1,1),col=rgb(1,0,0,0.5),lwd=2,lty=2) | |
| dev.off() | |
| # Plot differerntial expression plots | |
| pdf('MA_plots_DESeq.pdf') | |
| plotMA(res.cd44low.vs.klrg1_plus_pe_minus) | |
| plotMA(res.cd44low.vs.klrg1_plus_pe_plus) | |
| plotMA(res.cd44low.vs.pd1_plus_pe_minus) | |
| plotMA(res.klrg1_plus_pe_minus.vs.pd1_plus_pe_minus) | |
| plotMA(res.klrg1_plus_pe_plus.vs.pd1_plus_pe_minus) | |
| dev.off() | |
| # Load coverage without normalization | |
| d0 = read.csv('../gexp_coverage_NA.csv',header=T,row.names=1) | |
| d0.full = na.omit(d0) | |
| library(pheatmap) | |
| pheatmap(cor(d0.full,method='spearman'),color=colorRampPalette(c('blue4','white','gold'))(32)) | |
| # Load FPKMs | |
| dFPKM = read.csv('../gexp_fpkm_NA.csv',header=T,row.names=1) | |
| dFPKM.full = na.omit(dFPKM) | |
| library(pheatmap) | |
| pheatmap(cor(dFPKM.full,method='spearman'),color=colorRampPalette(c('blue4','white','gold'))(32)) | |
| # Load RNA-seq expression data | |
| d1 = read.table('gexp_coverage.normalized.NetworkCentralityScaling.tsv',header=T,row.names=1,sep='\t') | |
| d1 = d1[,-4] | |
| d1.full = na.omit(d1) | |
| library(pheatmap) | |
| c1 = cor(log10(d1.full),method='spearman') | |
| diag(c1) = NA | |
| pheatmap(c1,cluster_rows=F,cluster_cols=F,color=colorRampPalette(c('blue4','white','gold'))(32)) | |
| min(c1[1:4,1:4],na.rm=T) | |
| min(c1[5:9,5:9],na.rm=T) | |
| min(c1[10:12,10:12],na.rm=T) | |
| min(c1[13:7,13:17],na.rm=T) | |
| max(c1[1:4,1:4],na.rm=T) | |
| max(c1[5:9,5:9],na.rm=T) | |
| max(c1[10:12,10:12],na.rm=T) | |
| max(c1[13:7,13:17],na.rm=T) | |
| pdf('pairwisePlots.pdf') | |
| pairs(log10(d1[,1:4]),pch=18,col=rgb(0,0,1,0.5)) | |
| pairs(log10(d1[,5:9]),pch=18,col=rgb(0,0,1,0.5)) | |
| pairs(log10(d1[,10:12]),pch=18,col=rgb(0,0,1,0.5)) | |
| pairs(log10(d1[,13:17]),pch=18,col=rgb(0,0,1,0.5)) | |
| dev.off() | |
| # PCA analysis | |
| pca1 = princomp(d1.full) | |
| # With pseudo-counts | |
| d1.NoNA = d1 | |
| d1.NoNA[is.na(d1.NoNA)] = 0 | |
| d1.NoNA = d1.NoNA + 1 | |
| c1 = cor(log10(d1.NoNA),method='spearman') | |
| diag(c1) = NA | |
| #pheatmap(rbind(c1,c(seq(0,1,length.out=17))),cluster_rows=F,cluster_cols=F,color=colorRampPalette(c('blue4','white','gold'))(32)) | |
| pheatmap(c1,cluster_rows=F,cluster_cols=F,color=colorRampPalette(c('blue4','white','gold'))(32)) | |
| min(c1[1:4,1:4],na.rm=T) | |
| min(c1[5:9,5:9],na.rm=T) | |
| min(c1[10:12,10:12],na.rm=T) | |
| min(c1[13:7,13:17],na.rm=T) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment