Mercurial > repos > pravs > protein_rna_correlation

comparison protein_rna_correlation.r @ 6:8e9428eca82c 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:09:54 -0400
parents 6bf0203ee17e
children 80892c607b1e
comparison
equal deleted inserted replaced
5:6bf0203ee17e 6:8e9428eca82c
534 #============================================================================================================= 534 #=============================================================================================================
535 # Distribution of proteome and transcriptome abundance (Box and Density plot) 535 # Distribution of proteome and transcriptome abundance (Box and Density plot)
536 #============================================================================================================= 536 #=============================================================================================================
537 cat("Generating Box and Density plot\n",file=logfile, append=T); 537 cat("Generating Box and Density plot\n",file=logfile, append=T);
538 outplot = paste(outdir,"/AbundancePlot.png",sep="",collapse=""); 538 outplot = paste(outdir,"/AbundancePlot.png",sep="",collapse="");
539 png(outplot); 539 #png(outplot);
540 bitmap(outplot, "png16m");
540 par(mfrow=c(2,2)); 541 par(mfrow=c(2,2));
541 boxplot(proteome_df[,2], ylab="Abundance", main="Proteome abundance", cex.lab=1.5); 542 boxplot(proteome_df[,2], ylab="Abundance", main="Proteome abundance", cex.lab=1.5);
542 plot(density(proteome_df[,2]), xlab="Protein Abundance", ylab="Density", main="Proteome abundance", cex.lab=1.5); 543 plot(density(proteome_df[,2]), xlab="Protein Abundance", ylab="Density", main="Proteome abundance", cex.lab=1.5);
543 boxplot(transcriptome_df[,2], ylab="Abundance", main="Transcriptome abundance", cex.lab=1.5); 544 boxplot(transcriptome_df[,2], ylab="Abundance", main="Transcriptome abundance", cex.lab=1.5);
544 plot(density(transcriptome_df[,2]), xlab="Transcript Abundance", ylab="Density", main="Transcriptome abundance", cex.lab=1.5); 545 plot(density(transcriptome_df[,2]), xlab="Transcript Abundance", ylab="Density", main="Transcriptome abundance", cex.lab=1.5);
552 #============================================================================================================= 553 #=============================================================================================================
553 # Scatter plot 554 # Scatter plot
554 #============================================================================================================= 555 #=============================================================================================================
555 cat("Generating scatter plot\n",file=logfile, append=T); 556 cat("Generating scatter plot\n",file=logfile, append=T);
556 outplot = paste(outdir,"/AbundancePlot_scatter.png",sep="",collapse=""); 557 outplot = paste(outdir,"/AbundancePlot_scatter.png",sep="",collapse="");
557 png(outplot); 558 #png(outplot);
559 bitmap(outplot,"png16m")
558 par(mfrow=c(1,1)); 560 par(mfrow=c(1,1));
559 scatter.smooth(transcriptome_df[,2], proteome_df[,2], xlab="Transcript Abundance", ylab="Protein Abundance", cex.lab=1.5); 561 scatter.smooth(transcriptome_df[,2], proteome_df[,2], xlab="Transcript Abundance", ylab="Protein Abundance", cex.lab=1.5);
560 562
561 cat( 563 cat(
562 "<font color='blue'><h3>Scatter plot between Proteome and Transcriptome Abundance</h3></font>\n", 564 "<font color='blue'><h3>Scatter plot between Proteome and Transcriptome Abundance</h3></font>\n",
634 cat( 636 cat(
635 "<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",
636 file = htmloutfile, append = TRUE); 638 file = htmloutfile, append = TRUE);
637 639
638 outplot = paste(outdir,"/PE_GE_lm_1.png",sep="",collapse=""); 640 outplot = paste(outdir,"/PE_GE_lm_1.png",sep="",collapse="");
639 png(outplot); 641 #png(outplot);
642 bitmap(outplot,"png16m");
640 par(mfrow=c(1,1)); 643 par(mfrow=c(1,1));
641 plot(regmodel, 1, cex.lab=1.5); 644 plot(regmodel, 1, cex.lab=1.5);
642 dev.off(); 645 dev.off();
643 646
644 outplot = paste(outdir,"/PE_GE_lm_2.png",sep="",collapse=""); 647 outplot = paste(outdir,"/PE_GE_lm_2.png",sep="",collapse="");
645 png(outplot); 648 #png(outplot);
649 bitmap(outplot,"png16m");
646 par(mfrow=c(1,1)); 650 par(mfrow=c(1,1));
647 plot(regmodel, 2, cex.lab=1.5); 651 plot(regmodel, 2, cex.lab=1.5);
648 dev.off(); 652 dev.off();
649 653
650 outplot = paste(outdir,"/PE_GE_lm_3.png",sep="",collapse=""); 654 outplot = paste(outdir,"/PE_GE_lm_3.png",sep="",collapse="");
651 png(outplot); 655 #png(outplot);
656 bitmap(outplot,"png16m");
652 par(mfrow=c(1,1)); 657 par(mfrow=c(1,1));
653 plot(regmodel, 3, cex.lab=1.5); 658 plot(regmodel, 3, cex.lab=1.5);
654 dev.off(); 659 dev.off();
655 660
656 outplot = paste(outdir,"/PE_GE_lm_5.png",sep="",collapse=""); 661 outplot = paste(outdir,"/PE_GE_lm_5.png",sep="",collapse="");
657 png(outplot); 662 #png(outplot);
663 bitmap(outplot,"png16m");
658 par(mfrow=c(1,1)); 664 par(mfrow=c(1,1));
659 plot(regmodel, 5, cex.lab=1.5); 665 plot(regmodel, 5, cex.lab=1.5);
660 dev.off(); 666 dev.off();
661 667
662 outplot = paste(outdir,"/PE_GE_lm.pdf",sep="",collapse=""); 668 outplot = paste(outdir,"/PE_GE_lm.pdf",sep="",collapse="");
703 709
704 cooksd <- cooks.distance(regmodel); 710 cooksd <- cooks.distance(regmodel);
705 711
706 cat("Generating cooksd plot\n",file=logfile, append=T); 712 cat("Generating cooksd plot\n",file=logfile, append=T);
707 outplot = paste(outdir,"/PE_GE_lm_cooksd.png",sep="",collapse=""); 713 outplot = paste(outdir,"/PE_GE_lm_cooksd.png",sep="",collapse="");
708 png(outplot); 714 #png(outplot);
715 bitmap(outplot,"png16m");
709 par(mfrow=c(1,1)); 716 par(mfrow=c(1,1));
710 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
711 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
712 #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
713 dev.off(); 720 dev.off();
763 770
764 #============================================================================================================= 771 #=============================================================================================================
765 # Scatter plot 772 # Scatter plot
766 #============================================================================================================= 773 #=============================================================================================================
767 outplot = paste(outdir,"/AbundancePlot_scatter_without_outliers.png",sep="",collapse=""); 774 outplot = paste(outdir,"/AbundancePlot_scatter_without_outliers.png",sep="",collapse="");
768 png(outplot); 775 #png(outplot);
776 bitmap(outplot,"png16m");
769 par(mfrow=c(1,1)); 777 par(mfrow=c(1,1));
770 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);
771 779
772 cat( 780 cat(
773 "<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",
805 "<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",
806 file = htmloutfile, append = TRUE); 814 file = htmloutfile, append = TRUE);
807 815
808 cat("Generating heatmap plot\n",file=logfile, append=T); 816 cat("Generating heatmap plot\n",file=logfile, append=T);
809 outplot = paste(outdir,"/PE_GE_heatmap.png",sep="",collapse=""); 817 outplot = paste(outdir,"/PE_GE_heatmap.png",sep="",collapse="");
810 png(outplot); 818 #png(outplot);
819 bitmap(outplot,"png16m");
811 par(mfrow=c(1,1)); 820 par(mfrow=c(1,1));
812 #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");
813 my_palette <- colorRampPalette(c("green", "white", "red"))(299); 822 my_palette <- colorRampPalette(c("green", "white", "red"))(299);
814 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");
815 dev.off(); 824 dev.off();
826 835
827 836
828 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);
829 cat("Generating kmeans plot\n",file=logfile, append=T); 838 cat("Generating kmeans plot\n",file=logfile, append=T);
830 outplot = paste(outdir,"/PE_GE_kmeans.png",sep="",collapse=""); 839 outplot = paste(outdir,"/PE_GE_kmeans.png",sep="",collapse="");
831 png(outplot); 840 #png(outplot);
841 bitmap(outplot,"png16m");
832 par(mfrow=c(1,1)); 842 par(mfrow=c(1,1));
833 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);
834 844
835 ind=which(k1$cluster==1); 845 ind=which(k1$cluster==1);
836 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);