Mercurial > repos > iuc > isoformswitchanalyzer

comparison IsoformSwitchAnalyzeR.R @ 1:5ae218cee629 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/isoformswitchanalyzer commit 7d278967f9c8d2c6ce0a8f83be2c444822746bbf
author iuc
date Fri, 19 May 2023 21:25:39 +0000
parents bb611fa3bc3b
children db30020fee5b
comparison
equal deleted inserted replaced
0:bb611fa3bc3b 1:5ae218cee629
2 library(IsoformSwitchAnalyzeR, 2 library(IsoformSwitchAnalyzeR,
3 quietly = TRUE, 3 quietly = TRUE,
4 warn.conflicts = FALSE) 4 warn.conflicts = FALSE)
5 library(argparse, quietly = TRUE, warn.conflicts = FALSE) 5 library(argparse, quietly = TRUE, warn.conflicts = FALSE)
6 library(dplyr, quietly = TRUE, warn.conflicts = FALSE) 6 library(dplyr, quietly = TRUE, warn.conflicts = FALSE)
7 library(ggplot2, quietly = TRUE, warn.conflicts = FALSE)
8
7 9
8 # setup R error handling to go to stderr 10 # setup R error handling to go to stderr
9 options( 11 options(
10 show.error.messages = FALSE, 12 show.error.messages = FALSE,
11 error = function() { 13 error = function() {
738 740
739 pdf( 741 pdf(
740 file = outputFile, 742 file = outputFile,
741 onefile = FALSE, 743 onefile = FALSE,
742 height = 6, 744 height = 6,
743 width = 9 745 width = 12
744 ) 746 )
745 747
746 consequenceSummary <- extractConsequenceSummary( 748 consequenceSummary <- extractConsequenceSummary(
747 SwitchList, 749 SwitchList,
748 consequencesToAnalyze = "all", 750 consequencesToAnalyze = "all",
833 outputFile <- file.path(getwd(), "splicingSummary.pdf") 835 outputFile <- file.path(getwd(), "splicingSummary.pdf")
834 pdf( 836 pdf(
835 file = outputFile, 837 file = outputFile,
836 onefile = FALSE, 838 onefile = FALSE,
837 height = 6, 839 height = 6,
838 width = 9 840 width = 12
839 ) 841 )
840 splicingSummary <- extractSplicingSummary( 842 splicingSummary <- extractSplicingSummary(
841 SwitchList, 843 SwitchList,
842 splicingToAnalyze = "all", 844 splicingToAnalyze = "all",
843 asFractionTotal = args$asFractionTotal, 845 asFractionTotal = args$asFractionTotal,
861 ) 863 )
862 864
863 865
864 ### Volcano like plot: 866 ### Volcano like plot:
865 outputFile <- file.path(getwd(), "volcanoPlot.pdf") 867 outputFile <- file.path(getwd(), "volcanoPlot.pdf")
868
866 pdf( 869 pdf(
867 file = outputFile, 870 file = outputFile,
868 onefile = FALSE, 871 onefile = FALSE,
869 height = 6, 872 height = 6,
870 width = 9 873 width = 9
871 ) 874 )
872 ggplot(data = SwitchList$isoformFeatures, aes(x = dIF, y = -log10(isoform_switch_q_value))) + 875
873 geom_point(aes(color = abs(dIF) > 0.1 & 876 p <- ggplot(data = SwitchList$isoformFeatures, aes(x = dIF, y = -log10(isoform_switch_q_value))) +
874 isoform_switch_q_value < 0.05), # default cutoff 877 geom_point(
875 size = 1) + 878 aes(color = abs(dIF) > 0.1 & isoform_switch_q_value < 0.05), # default cutoff
879 size = 1
880 ) +
876 geom_hline(yintercept = -log10(0.05), linetype = "dashed") + # default cutoff 881 geom_hline(yintercept = -log10(0.05), linetype = "dashed") + # default cutoff
877 geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") + # default cutoff 882 geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") + # default cutoff
878 facet_wrap(~ condition_2) + 883 facet_wrap(~ condition_2) +
884 #facet_grid(condition_1 ~ condition_2) + # alternative to facet_wrap if you have overlapping conditions
879 scale_color_manual("Signficant\nIsoform Switch", values = c("black", "red")) + 885 scale_color_manual("Signficant\nIsoform Switch", values = c("black", "red")) +
880 labs(x = "dIF", y = "-Log10 ( Isoform Switch Q Value )") + 886 labs(x = "dIF", y = "-Log10 ( Isoform Switch Q Value )") +
881 theme_bw() 887 theme_bw()
888 print(p)
882 dev.off() 889 dev.off()
883
884 890
885 ### Switch vs Gene changes: 891 ### Switch vs Gene changes:
886 outputFile <- file.path(getwd(), "switchGene.pdf") 892 outputFile <- file.path(getwd(), "switchGene.pdf")
887 pdf( 893 pdf(
888 file = outputFile, 894 file = outputFile,
889 onefile = FALSE,
890 height = 6, 895 height = 6,
891 width = 9 896 width = 9
892 ) 897 )
893 ggplot(data = SwitchList$isoformFeatures, 898 p <- ggplot(data = SwitchList$isoformFeatures,
894 aes(x = gene_log2_fold_change, y = dIF)) + 899 aes(x = gene_log2_fold_change, y = dIF)) +
895 geom_point(aes(color = abs(dIF) > 0.1 & 900 geom_point(aes(color = abs(dIF) > 0.1 &
896 isoform_switch_q_value < 0.05), 901 isoform_switch_q_value < 0.05),
897 size = 1) + 902 size = 1) +
898 facet_wrap(~ condition_2) + 903 facet_wrap(~ condition_2) +
899 geom_hline(yintercept = 0, linetype = "dashed") + 904 geom_hline(yintercept = 0, linetype = "dashed") +
900 geom_vline(xintercept = 0, linetype = "dashed") + 905 geom_vline(xintercept = 0, linetype = "dashed") +
901 scale_color_manual("Signficant\nIsoform Switch", values = c("black", "red")) + 906 scale_color_manual("Signficant\nIsoform Switch", values = c("black", "red")) +
902 labs(x = "Gene log2 fold change", y = "dIF") + 907 labs(x = "Gene log2 fold change", y = "dIF") +
903 theme_bw() 908 theme_bw()
909 print(p)
904 dev.off() 910 dev.off()
905 911
906 outputFile <- file.path(getwd(), "splicingGenomewide.pdf") 912 outputFile <- file.path(getwd(), "splicingGenomewide.pdf")
907 pdf( 913 pdf(
908 file = outputFile, 914 file = outputFile,
909 onefile = FALSE, 915 onefile = FALSE,
910 height = 6, 916 height = 6,
911 width = 9 917 width = 14
912 ) 918 )
913 splicingGenomeWide <- extractSplicingGenomeWide( 919 splicingGenomeWide <- extractSplicingGenomeWide(
914 SwitchList, 920 SwitchList,
915 featureToExtract = "all", 921 featureToExtract = "isoformUsage",
916 splicingToAnalyze = c("A3", "MES", "ATSS"), 922 splicingToAnalyze = "all",
917 plot = TRUE, 923 plot = TRUE,
918 returnResult = TRUE 924 returnResult = TRUE
919 ) 925 )
920 dev.off() 926 dev.off()
921 } 927 }