Mercurial > repos > pravs > protein_rna_correlation

comparison protein_rna_correlation.r @ 9:e407b1a7a8de draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author pravs
date Wed, 20 Jun 2018 21:42:36 -0400
parents 80892c607b1e
children
comparison
equal deleted inserted replaced
8:80892c607b1e 9:e407b1a7a8de
636 cat( 636 cat(
637 "<font color='blue'><h3>Plotting various regression diagnostics plots</h3></font>\n", 637 "<font color='blue'><h3>Plotting various regression diagnostics plots</h3></font>\n",
638 file = htmloutfile, append = TRUE); 638 file = htmloutfile, append = TRUE);
639 639
640 outplot = paste(outdir,"/PE_GE_lm_1.png",sep="",collapse=""); 640 outplot = paste(outdir,"/PE_GE_lm_1.png",sep="",collapse="");
641 png(outplot); 641 #png(outplot);
642 #bitmap(outplot,"png16m"); 642 bitmap(outplot,"png16m");
643 par(mfrow=c(1,1)); 643 par(mfrow=c(1,1));
644 plot(regmodel, 1, cex.lab=1.5); 644 plot(regmodel, 1, cex.lab=1.5);
645 dev.off(); 645 dev.off();
646 646
647 outplot = paste(outdir,"/PE_GE_lm_2.png",sep="",collapse=""); 647 outplot = paste(outdir,"/PE_GE_lm_2.png",sep="",collapse="");
648 png(outplot); 648 #png(outplot);
649 #bitmap(outplot,"png16m"); 649 bitmap(outplot,"png16m");
650 par(mfrow=c(1,1)); 650 par(mfrow=c(1,1));
651 plot(regmodel, 2, cex.lab=1.5); 651 plot(regmodel, 2, cex.lab=1.5);
652 dev.off(); 652 dev.off();
653 653
654 outplot = paste(outdir,"/PE_GE_lm_3.png",sep="",collapse=""); 654 outplot = paste(outdir,"/PE_GE_lm_3.png",sep="",collapse="");
655 png(outplot); 655 #png(outplot);
656 #bitmap(outplot,"png16m"); 656 bitmap(outplot,"png16m");
657 par(mfrow=c(1,1)); 657 par(mfrow=c(1,1));
658 plot(regmodel, 3, cex.lab=1.5); 658 plot(regmodel, 3, cex.lab=1.5);
659 dev.off(); 659 dev.off();
660 660
661 outplot = paste(outdir,"/PE_GE_lm_5.png",sep="",collapse=""); 661 outplot = paste(outdir,"/PE_GE_lm_5.png",sep="",collapse="");
662 png(outplot); 662 #png(outplot);
663 #bitmap(outplot,"png16m"); 663 bitmap(outplot,"png16m");
664 par(mfrow=c(1,1)); 664 par(mfrow=c(1,1));
665 plot(regmodel, 5, cex.lab=1.5); 665 plot(regmodel, 5, cex.lab=1.5);
666 dev.off(); 666 dev.off();
667 667
668 outplot = paste(outdir,"/PE_GE_lm.pdf",sep="",collapse=""); 668 outplot = paste(outdir,"/PE_GE_lm.pdf",sep="",collapse="");
709 709
710 cooksd <- cooks.distance(regmodel); 710 cooksd <- cooks.distance(regmodel);
711 711
712 cat("Generating cooksd plot\n",file=logfile, append=T); 712 cat("Generating cooksd plot\n",file=logfile, append=T);
713 outplot = paste(outdir,"/PE_GE_lm_cooksd.png",sep="",collapse=""); 713 outplot = paste(outdir,"/PE_GE_lm_cooksd.png",sep="",collapse="");
714 png(outplot); 714 #png(outplot);
715 #bitmap(outplot,"png16m"); 715 bitmap(outplot,"png16m");
716 par(mfrow=c(1,1)); 716 par(mfrow=c(1,1));
717 plot(cooksd, pch="*", cex=2, cex.lab=1.5,main="Influential Obs. by Cooks distance", ylab="Cook\'s distance", xlab="Observations") # plot cooks distance 717 plot(cooksd, pch="*", cex=2, cex.lab=1.5,main="Influential Obs. by Cooks distance", ylab="Cook\'s distance", xlab="Observations") # plot cooks distance
718 abline(h = 4*mean(cooksd, na.rm=T), col="red") # add cutoff line 718 abline(h = 4*mean(cooksd, na.rm=T), col="red") # add cutoff line
719 #text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""), col="red", pos=2) # add labels 719 #text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""), col="red", pos=2) # add labels
720 dev.off(); 720 dev.off();
770 770
771 #============================================================================================================= 771 #=============================================================================================================
772 # Scatter plot 772 # Scatter plot
773 #============================================================================================================= 773 #=============================================================================================================
774 outplot = paste(outdir,"/AbundancePlot_scatter_without_outliers.png",sep="",collapse=""); 774 outplot = paste(outdir,"/AbundancePlot_scatter_without_outliers.png",sep="",collapse="");
775 png(outplot); 775 #png(outplot);
776 #bitmap(outplot,"png16m"); 776 bitmap(outplot,"png16m");
777 par(mfrow=c(1,1)); 777 par(mfrow=c(1,1));
778 scatter.smooth(PE_GE_data_no_outlier[,"GE_abundance"], PE_GE_data_no_outlier[,"PE_abundance"], xlab="Transcript Abundance", ylab="Protein Abundance", cex.lab=1.5); 778 scatter.smooth(PE_GE_data_no_outlier[,"GE_abundance"], PE_GE_data_no_outlier[,"PE_abundance"], xlab="Transcript Abundance", ylab="Protein Abundance", cex.lab=1.5);
779 779
780 cat( 780 cat(
781 "<font color='blue'><h3>Scatter plot between Proteome and Transcriptome Abundance, after removal of outliers/influential observations</h3></font>\n", 781 "<font color='blue'><h3>Scatter plot between Proteome and Transcriptome Abundance, after removal of outliers/influential observations</h3></font>\n",
813 "<font color='blue'><h3>Heatmap of PE and GE abundance values</h3></font>\n", 813 "<font color='blue'><h3>Heatmap of PE and GE abundance values</h3></font>\n",
814 file = htmloutfile, append = TRUE); 814 file = htmloutfile, append = TRUE);
815 815
816 cat("Generating heatmap plot\n",file=logfile, append=T); 816 cat("Generating heatmap plot\n",file=logfile, append=T);
817 outplot = paste(outdir,"/PE_GE_heatmap.png",sep="",collapse=""); 817 outplot = paste(outdir,"/PE_GE_heatmap.png",sep="",collapse="");
818 png(outplot); 818 #png(outplot);
819 #bitmap(outplot,"png16m"); 819 bitmap(outplot,"png16m");
820 par(mfrow=c(1,1)); 820 par(mfrow=c(1,1));
821 #heatmap.2(as.matrix(PE_GE_data[,c("PE_abundance","GE_abundance")]), trace="none", cexCol=1, col=greenred(100),Colv=F, labCol=c("PE","GE"), scale="col"); 821 #heatmap.2(as.matrix(PE_GE_data[,c("PE_abundance","GE_abundance")]), trace="none", cexCol=1, col=greenred(100),Colv=F, labCol=c("PE","GE"), scale="col");
822 my_palette <- colorRampPalette(c("green", "white", "red"))(299); 822 my_palette <- colorRampPalette(c("green", "white", "red"))(299);
823 heatmap.2(as.matrix(PE_GE_data[,c("PE_abundance","GE_abundance")]), trace="none", cexCol=1, col=my_palette ,Colv=F, labCol=c("PE","GE"), scale="col", dendrogram = "row"); 823 heatmap.2(as.matrix(PE_GE_data[,c("PE_abundance","GE_abundance")]), trace="none", cexCol=1, col=my_palette ,Colv=F, labCol=c("PE","GE"), scale="col", dendrogram = "row");
824 dev.off(); 824 dev.off();
835 835
836 836
837 k1 = kmeans(PE_GE_data_kdata[,c("PE_abundance","GE_abundance")], 5); 837 k1 = kmeans(PE_GE_data_kdata[,c("PE_abundance","GE_abundance")], 5);
838 cat("Generating kmeans plot\n",file=logfile, append=T); 838 cat("Generating kmeans plot\n",file=logfile, append=T);
839 outplot = paste(outdir,"/PE_GE_kmeans.png",sep="",collapse=""); 839 outplot = paste(outdir,"/PE_GE_kmeans.png",sep="",collapse="");
840 png(outplot); 840 #png(outplot);
841 #bitmap(outplot,"png16m"); 841 bitmap(outplot,"png16m");
842 par(mfrow=c(1,1)); 842 par(mfrow=c(1,1));
843 scatter.smooth(PE_GE_data_kdata[,"GE_abundance"], PE_GE_data_kdata[,"PE_abundance"], xlab="Transcript Abundance", ylab="Protein Abundance", cex.lab=1.5); 843 scatter.smooth(PE_GE_data_kdata[,"GE_abundance"], PE_GE_data_kdata[,"PE_abundance"], xlab="Transcript Abundance", ylab="Protein Abundance", cex.lab=1.5);
844 844
845 ind=which(k1$cluster==1); 845 ind=which(k1$cluster==1);
846 points(PE_GE_data_kdata[ind,"GE_abundance"], PE_GE_data_kdata[ind,"PE_abundance"], col="red", pch=16); 846 points(PE_GE_data_kdata[ind,"GE_abundance"], PE_GE_data_kdata[ind,"PE_abundance"], col="red", pch=16);