Mercurial > repos > eschen42 > w4mcorcov

comparison w4mcorcov_calc.R @ 11:ddcc33ff3205 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/HegemanLab/w4mcorcov_galaxy_wrapper/tree/master commit 4428e3252d54c8a8e0e5d85e8eaaeb13e9b21de7
author eschen42
date Wed, 05 Sep 2018 22:31:21 -0400
parents 066b1f409e9f
children ddaf84e15d06
comparison
equal deleted inserted replaced
10:9a52306991b3 11:ddcc33ff3205
69 my_cor_vs_cov <- cor_vs_cov_x 69 my_cor_vs_cov <- cor_vs_cov_x
70 } 70 }
71 # print("str(my_cor_vs_cov)") 71 # print("str(my_cor_vs_cov)")
72 # str(my_cor_vs_cov) 72 # str(my_cor_vs_cov)
73 if (is.null(my_cor_vs_cov) || sum(!is.na(my_cor_vs_cov$tsv1$covariance)) < 2) { 73 if (is.null(my_cor_vs_cov) || sum(!is.na(my_cor_vs_cov$tsv1$covariance)) < 2) {
74 x_progress("No cor_vs_cov data produced") 74 if (is.null(cor_vs_cov_x)) {
75 x_progress("No cor_vs_cov data produced")
76 }
75 plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n") 77 plot(x=1, y=1, xaxt="n", yaxt="n", xlab="", ylab="", type="n")
76 text(x=1, y=1, labels="too few covariance data") 78 text(x=1, y=1, labels="too few covariance data")
77 return(my_cor_vs_cov) 79 return(my_cor_vs_cov)
78 } 80 }
79 with( 81 with(
559 col_selector <- vrbl_metadata_names[ overall_significant ] 561 col_selector <- vrbl_metadata_names[ overall_significant ]
560 my_matrix <- tdm[ chosen_samples, col_selector, drop = FALSE ] 562 my_matrix <- tdm[ chosen_samples, col_selector, drop = FALSE ]
561 plot_action <- function(fctr_lvl_1, fctr_lvl_2) { 563 plot_action <- function(fctr_lvl_1, fctr_lvl_2) {
562 progress_action( 564 progress_action(
563 sprintf("calculating/plotting contrast of %s vs. %s" 565 sprintf("calculating/plotting contrast of %s vs. %s"
564 , fctr_lvl_1, fctr_lvl_2)) 566 , fctr_lvl_1, fctr_lvl_2)
567 )
565 predictor <- sapply( 568 predictor <- sapply(
566 X = chosen_facC 569 X = chosen_facC
567 , FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2 570 , FUN = function(fac) if ( fac == fctr_lvl_1 ) fctr_lvl_1 else fctr_lvl_2
568 ) 571 )
569 my_cor_cov <- do_detail_plot( 572 my_cor_cov <- do_detail_plot(
570 x_dataMatrix = my_matrix 573 x_dataMatrix = my_matrix
571 , x_predictor = predictor 574 , x_predictor = predictor
572 , x_is_match = is_match 575 , x_is_match = TRUE
573 , x_algorithm = algoC 576 , x_algorithm = algoC
574 , x_prefix = if (pairSigFeatOnly) { 577 , x_prefix = if (pairSigFeatOnly) {
575 "Significantly contrasting features" 578 "Significantly contrasting features"
576 } else { 579 } else {
577 "Significant features" 580 "Significant features"
580 , x_progress = progress_action 583 , x_progress = progress_action
581 , x_crossval_i = min(7, length(chosen_samples)) 584 , x_crossval_i = min(7, length(chosen_samples))
582 , x_env = calc_env 585 , x_env = calc_env
583 ) 586 )
584 if ( is.null(my_cor_cov) ) { 587 if ( is.null(my_cor_cov) ) {
585 progress_action("NOTHING TO PLOT.") 588 progress_action("NOTHING TO PLOT")
586 } else { 589 } else {
587 my_tsv <- my_cor_cov$tsv1 590 my_tsv <- my_cor_cov$tsv1
588 my_tsv$mz <- mz_lookup(my_tsv$featureID) 591 my_tsv$mz <- mz_lookup(my_tsv$featureID)
589 my_tsv$rt <- rt_lookup(my_tsv$featureID) 592 my_tsv$rt <- rt_lookup(my_tsv$featureID)
590 my_tsv["level1Level2Sig"] <- vrbl_metadata[ 593 my_tsv["level1Level2Sig"] <- vrbl_metadata[
617 vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name 620 vrbl_metadata_col <- col_match[1] # ^^^^^^^^^^^^^^^^^^^^^ # Column name
618 fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1 621 fctr_lvl_1 <- col_match[2] # ^^ # Factor-level 1
619 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2 622 fctr_lvl_2 <- col_match[3] # ^^ # Factor-level 2
620 # only process this column if both factors are members of lvlCSV 623 # only process this column if both factors are members of lvlCSV
621 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2) 624 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
622 progress_action( 625 if (is_match) {
623 sprintf("calculating/plotting contrast of %s vs. %s" 626 progress_action(
624 , fctr_lvl_1, fctr_lvl_2)) 627 sprintf("calculating/plotting contrast of %s vs. %s."
625 # choose only samples with one of the two factors for this column 628 , fctr_lvl_1, fctr_lvl_2
626 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) 629 )
627 predictor <- smpl_metadata_facC[chosen_samples]
628 # extract only the significantly-varying features and the chosen samples
629 fully_significant <- 1 == vrbl_metadata[,vrbl_metadata_col] *
630 ( if (intersample_sig_col %in% colnames(vrbl_metadata)) {
631 vrbl_metadata[,intersample_sig_col]
632 } else {
633 1
634 }
635 ) 630 )
636 col_selector <- vrbl_metadata_names[ 631 # choose only samples with one of the two factors for this column
637 if ( pairSigFeatOnly ) fully_significant else overall_significant 632 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
638 ] 633 predictor <- smpl_metadata_facC[chosen_samples]
639 my_matrix <- tdm[ chosen_samples, col_selector, drop = FALSE ] 634 # extract only the significantly-varying features and the chosen samples
640 my_cor_cov <- do_detail_plot( 635 fully_significant <- 1 == vrbl_metadata[,vrbl_metadata_col] *
641 x_dataMatrix = my_matrix 636 ( if (intersample_sig_col %in% colnames(vrbl_metadata)) {
642 , x_predictor = predictor 637 vrbl_metadata[,intersample_sig_col]
643 , x_is_match = is_match 638 } else {
644 , x_algorithm = algoC 639 1
645 , x_prefix = if (pairSigFeatOnly) { 640 }
646 "Significantly contrasting features" 641 )
647 } else { 642 col_selector <- vrbl_metadata_names[
648 "Significant features" 643 if ( pairSigFeatOnly ) fully_significant else overall_significant
649 } 644 ]
650 , x_show_labels = labelFeatures 645 my_matrix <- tdm[ chosen_samples, col_selector, drop = FALSE ]
651 , x_progress = progress_action 646 my_cor_cov <- do_detail_plot(
652 , x_crossval_i = min(7, length(chosen_samples)) 647 x_dataMatrix = my_matrix
653 , x_env = calc_env 648 , x_predictor = predictor
654 ) 649 , x_is_match = is_match
655 if ( is.null(my_cor_cov) ) { 650 , x_algorithm = algoC
656 progress_action("NOTHING TO PLOT.") 651 , x_prefix = if (pairSigFeatOnly) {
652 "Significantly contrasting features"
653 } else {
654 "Significant features"
655 }
656 , x_show_labels = labelFeatures
657 , x_progress = progress_action
658 , x_crossval_i = min(7, length(chosen_samples))
659 , x_env = calc_env
660 )
661 if ( is.null(my_cor_cov) ) {
662 progress_action("NOTHING TO PLOT.")
663 } else {
664 tsv <- my_cor_cov$tsv1
665 tsv$mz <- mz_lookup(tsv$featureID)
666 tsv$rt <- rt_lookup(tsv$featureID)
667 tsv["level1Level2Sig"] <- vrbl_metadata[
668 match(tsv$featureID, vrbl_metadata_names)
669 , vrbl_metadata_col
670 ]
671 corcov_tsv_action(tsv)
672 did_plot <- TRUE
673 }
657 } else { 674 } else {
658 tsv <- my_cor_cov$tsv1 675 progress_action(
659 tsv$mz <- mz_lookup(tsv$featureID) 676 sprintf("skipping contrast of %s vs. %s."
660 tsv$rt <- rt_lookup(tsv$featureID) 677 , fctr_lvl_1, fctr_lvl_2
661 tsv["level1Level2Sig"] <- vrbl_metadata[ 678 )
662 match(tsv$featureID, vrbl_metadata_names) 679 )
663 , vrbl_metadata_col
664 ]
665 corcov_tsv_action(tsv)
666 did_plot <- TRUE
667 } 680 }
668 } 681 }
669 } 682 }
670 } 683 }
671 } else { # tesC == "none" 684 } else { # tesC == "none"
685 # find all the levels for factor facC in sampleMetadata
672 level_union <- unique(sort(smpl_metadata_facC)) 686 level_union <- unique(sort(smpl_metadata_facC))
687 # identify the selected levels for factor facC from sampleMetadata
688 level_include <- sapply(X = level_union, FUN = isLevelSelected)
689 # discard the non-selected levels for factor facC
690 level_union <- level_union[level_include]
673 if ( length(level_union) > 1 ) { 691 if ( length(level_union) > 1 ) {
674 if ( length(level_union) > 2 ) { 692 if ( length(level_union) > 2 ) {
675 ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ## 693 ## pass 1 - contrast each selected level with all other levels combined into one "super-level" ##
676 completed <- c() 694 completed <- c()
677 lapply( 695 lapply(
685 completed <<- c(completed, fctr_lvl_1) 703 completed <<- c(completed, fctr_lvl_1)
686 setdiff(level_union, fctr_lvl_1) 704 setdiff(level_union, fctr_lvl_1)
687 } 705 }
688 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2) 706 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
689 fctr_lvl_2 <- "other" 707 fctr_lvl_2 <- "other"
690 progress_action(
691 sprintf("calculating/plotting contrast of %s vs. %s"
692 , fctr_lvl_1, fctr_lvl_2)
693 )
694 if (length(unique(chosen_samples)) < 1) { 708 if (length(unique(chosen_samples)) < 1) {
695 progress_action("NOTHING TO PLOT...") 709 progress_action(
710 sprintf("Skipping contrast of %s vs. %s; there are no chosen samples."
711 , fctr_lvl_1, fctr_lvl_2)
712 )
696 } else { 713 } else {
697 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples]) 714 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
698 predictor <- sapply( 715 predictor <- sapply(
699 X = chosen_facC 716 X = chosen_facC
700 , FUN = function(fac) { 717 , FUN = function(fac) {
702 } 719 }
703 ) 720 )
704 my_matrix <- tdm[ chosen_samples, , drop = FALSE ] 721 my_matrix <- tdm[ chosen_samples, , drop = FALSE ]
705 # only process this column if both factors are members of lvlCSV 722 # only process this column if both factors are members of lvlCSV
706 is_match <- isLevelSelected(fctr_lvl_1) 723 is_match <- isLevelSelected(fctr_lvl_1)
724 if (is_match) {
725 progress_action(
726 sprintf("Calculating/plotting contrast of %s vs. %s"
727 , fctr_lvl_1, fctr_lvl_2)
728 )
729 my_cor_cov <- do_detail_plot(
730 x_dataMatrix = my_matrix
731 , x_predictor = predictor
732 , x_is_match = is_match
733 , x_algorithm = algoC
734 , x_prefix = "Features"
735 , x_show_labels = labelFeatures
736 , x_progress = progress_action
737 , x_crossval_i = min(7, length(chosen_samples))
738 , x_env = calc_env
739 )
740 if ( is.null(my_cor_cov) ) {
741 progress_action("NOTHING TO PLOT...")
742 } else {
743 tsv <- my_cor_cov$tsv1
744 tsv$mz <- mz_lookup(tsv$featureID)
745 tsv$rt <- rt_lookup(tsv$featureID)
746 corcov_tsv_action(tsv)
747 did_plot <<- TRUE
748 }
749 } else {
750 }
751 }
752 "dummy" # need to return a value; otherwise combn fails with an error
753 }
754 )
755 }
756 ## pass 2 - contrast each selected level with each of the other levels individually ##
757 completed <- c()
758 utils::combn(
759 x = level_union
760 , m = 2
761 , FUN = function(x) {
762 fctr_lvl_1 <- x[1]
763 fctr_lvl_2 <- x[2]
764 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
765 if (length(unique(chosen_samples)) < 1) {
766 progress_action(
767 sprintf("Skipping contrast of %s vs. %s. - There are no chosen samples."
768 , fctr_lvl_1, fctr_lvl_2
769 )
770 )
771 } else {
772 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
773 predictor <- chosen_facC
774 my_matrix <- tdm[ chosen_samples, , drop = FALSE ]
775 # only process this column if both factors are members of lvlCSV
776 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
777 if (is_match) {
778 progress_action(
779 sprintf("Calculating/plotting contrast of %s vs. %s."
780 , fctr_lvl_1, fctr_lvl_2)
781 )
707 my_cor_cov <- do_detail_plot( 782 my_cor_cov <- do_detail_plot(
708 x_dataMatrix = my_matrix 783 x_dataMatrix = my_matrix
709 , x_predictor = predictor 784 , x_predictor = predictor
710 , x_is_match = is_match 785 , x_is_match = is_match
711 , x_algorithm = algoC 786 , x_algorithm = algoC
714 , x_progress = progress_action 789 , x_progress = progress_action
715 , x_crossval_i = min(7, length(chosen_samples)) 790 , x_crossval_i = min(7, length(chosen_samples))
716 , x_env = calc_env 791 , x_env = calc_env
717 ) 792 )
718 if ( is.null(my_cor_cov) ) { 793 if ( is.null(my_cor_cov) ) {
719 progress_action("NOTHING TO PLOT") 794 progress_action("NOTHING TO PLOT.....")
720 } else { 795 } else {
721 tsv <- my_cor_cov$tsv1 796 tsv <- my_cor_cov$tsv1
722 tsv$mz <- mz_lookup(tsv$featureID) 797 tsv$mz <- mz_lookup(tsv$featureID)
723 tsv$rt <- rt_lookup(tsv$featureID) 798 tsv$rt <- rt_lookup(tsv$featureID)
724 corcov_tsv_action(tsv) 799 corcov_tsv_action(tsv)
725 did_plot <<- TRUE 800 did_plot <<- TRUE
726 } 801 }
727 }
728 "dummy" # need to return a value; otherwise combn fails with an error
729 }
730 )
731 }
732 ## pass 2 - contrast each selected level with each of the other levels individually ##
733 completed <- c()
734 utils::combn(
735 x = level_union
736 , m = 2
737 , FUN = function(x) {
738 fctr_lvl_1 <- x[1]
739 fctr_lvl_2 <- x[2]
740 chosen_samples <- smpl_metadata_facC %in% c(fctr_lvl_1, fctr_lvl_2)
741 progress_action(
742 sprintf("calculating/plotting contrast of %s vs. %s"
743 , fctr_lvl_1, fctr_lvl_2))
744 if (length(unique(chosen_samples)) < 1) {
745 progress_action("NOTHING TO PLOT...")
746 } else {
747 chosen_facC <- as.character(smpl_metadata_facC[chosen_samples])
748 predictor <- chosen_facC
749 my_matrix <- tdm[ chosen_samples, , drop = FALSE ]
750 # only process this column if both factors are members of lvlCSV
751 is_match <- isLevelSelected(fctr_lvl_1) && isLevelSelected(fctr_lvl_2)
752 my_cor_cov <- do_detail_plot(
753 x_dataMatrix = my_matrix
754 , x_predictor = predictor
755 , x_is_match = is_match
756 , x_algorithm = algoC
757 , x_prefix = "Features"
758 , x_show_labels = labelFeatures
759 , x_progress = progress_action
760 , x_crossval_i = min(7, length(chosen_samples))
761 , x_env = calc_env
762 )
763 if ( is.null(my_cor_cov) ) {
764 progress_action("NOTHING TO PLOT")
765 } else { 802 } else {
766 tsv <- my_cor_cov$tsv1 803 progress_action(
767 tsv$mz <- mz_lookup(tsv$featureID) 804 sprintf("Skipping contrast of %s vs. %s."
768 tsv$rt <- rt_lookup(tsv$featureID) 805 , fctr_lvl_1, fctr_lvl_2
769 corcov_tsv_action(tsv) 806 )
770 did_plot <<- TRUE 807 )
771 } 808 }
772 } 809 }
773 "dummy" # need to return a value; otherwise combn fails with an error 810 "dummy" # need to return a value; otherwise combn fails with an error
774 } 811 }
775 ) 812 )
776 } else { 813 } else {
777 progress_action("NOTHING TO PLOT....") 814 progress_action("NOTHING TO PLOT......")
778 } 815 }
779 } 816 }
780 if (!did_plot) { 817 if (!did_plot) {
781 failure_action( 818 failure_action(
782 sprintf( 819 sprintf(