# HG changeset patch # User fubar # Date 1370757604 14400 # Node ID f24f30eae81954472bf2573807eed7a075e409f3 # Parent aef0ecc0068312b0b69ab9191e65c40bf8e87323 Uploaded diff -r aef0ecc00683 -r f24f30eae819 gseaChip_3.1_symbols.loc.sample --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gseaChip_3.1_symbols.loc.sample Sun Jun 09 02:00:04 2013 -0400 @@ -0,0 +1,141 @@ +# generated by makeGSEAloc.py at 24/01/2013 10:43:16 from path = /data/genomes/gsea/3.1 +# loc file for the rgGSEA tool pointing to files from ftp://gseaftp.broadinstitute.org/pub/gsea/annotations +# name value path +AFFYMETRIX AFFYMETRIX /data/genomes/gsea/3.1/AFFYMETRIX.chip +AG AG /data/genomes/gsea/3.1/AG.chip +APPLERA_ABI1700 APPLERA_ABI1700 /data/genomes/gsea/3.1/APPLERA_ABI1700.chip +ATH1_121501 ATH1_121501 /data/genomes/gsea/3.1/ATH1_121501.chip +Agilent_Human1A Agilent_Human1A /data/genomes/gsea/3.1/Agilent_Human1A.chip +Agilent_Human1Av2 Agilent_Human1Av2 /data/genomes/gsea/3.1/Agilent_Human1Av2.chip +Agilent_Human1B Agilent_Human1B /data/genomes/gsea/3.1/Agilent_Human1B.chip +Agilent_Human1_cDNA Agilent_Human1_cDNA /data/genomes/gsea/3.1/Agilent_Human1_cDNA.chip +Agilent_HumanGenome Agilent_HumanGenome /data/genomes/gsea/3.1/Agilent_HumanGenome.chip +Agilent_MouseDev Agilent_MouseDev /data/genomes/gsea/3.1/Agilent_MouseDev.chip +Agilent_MouseGenome Agilent_MouseGenome /data/genomes/gsea/3.1/Agilent_MouseGenome.chip +Agilent_MouseOligo Agilent_MouseOligo /data/genomes/gsea/3.1/Agilent_MouseOligo.chip +Agilent_Mouse_cDNA Agilent_Mouse_cDNA /data/genomes/gsea/3.1/Agilent_Mouse_cDNA.chip +Agilent_RatGenome_G4131A Agilent_RatGenome_G4131A /data/genomes/gsea/3.1/Agilent_RatGenome_G4131A.chip +Agilent_RatOligo Agilent_RatOligo /data/genomes/gsea/3.1/Agilent_RatOligo.chip +AtlasMouse1.2 AtlasMouse1.2 /data/genomes/gsea/3.1/AtlasMouse1.2.chip +AtlasRat1.2 AtlasRat1.2 /data/genomes/gsea/3.1/AtlasRat1.2.chip +BIAO BIAO /data/genomes/gsea/3.1/BIAO.chip +Bovine Bovine /data/genomes/gsea/3.1/Bovine.chip +CNMCMuscleChip CNMCMuscleChip /data/genomes/gsea/3.1/CNMCMuscleChip.chip +Caltech16KMouse Caltech16KMouse /data/genomes/gsea/3.1/Caltech16KMouse.chip +Caltech16KOligoMouse Caltech16KOligoMouse /data/genomes/gsea/3.1/Caltech16KOligoMouse.chip +Canine_2 Canine_2 /data/genomes/gsea/3.1/Canine_2.chip +Clontech_Atlas_13K Clontech_Atlas_13K /data/genomes/gsea/3.1/Clontech_Atlas_13K.chip +Clontech_BD_Atlas Clontech_BD_Atlas /data/genomes/gsea/3.1/Clontech_BD_Atlas.chip +CodeLink_Human_Whole_Genome CodeLink_Human_Whole_Genome /data/genomes/gsea/3.1/CodeLink_Human_Whole_Genome.chip +CodeLink_UniSet_Human_20K_I_Bioarray CodeLink_UniSet_Human_20K_I_Bioarray /data/genomes/gsea/3.1/CodeLink_UniSet_Human_20K_I_Bioarray.chip +CodeLink_UniSet_Human_II_Bioarray CodeLink_UniSet_Human_II_Bioarray /data/genomes/gsea/3.1/CodeLink_UniSet_Human_II_Bioarray.chip +CodeLink_UniSet_Human_I_Bioarray CodeLink_UniSet_Human_I_Bioarray /data/genomes/gsea/3.1/CodeLink_UniSet_Human_I_Bioarray.chip +CodeLink_UniSet_Rat_I_Bioarray CodeLink_UniSet_Rat_I_Bioarray /data/genomes/gsea/3.1/CodeLink_UniSet_Rat_I_Bioarray.chip +DrosGenome1 DrosGenome1 /data/genomes/gsea/3.1/DrosGenome1.chip +Drosophila_2 Drosophila_2 /data/genomes/gsea/3.1/Drosophila_2.chip +G4110A G4110A /data/genomes/gsea/3.1/G4110A.chip +G4110Av2 G4110Av2 /data/genomes/gsea/3.1/G4110Av2.chip +GENE_SYMBOL GENE_SYMBOL /data/genomes/gsea/3.1/GENE_SYMBOL.chip +GenosysCytokineV2 GenosysCytokineV2 /data/genomes/gsea/3.1/GenosysCytokineV2.chip +HC_G110 HC_G110 /data/genomes/gsea/3.1/HC_G110.chip +HG_Focus HG_Focus /data/genomes/gsea/3.1/HG_Focus.chip +HG_U133A HG_U133A /data/genomes/gsea/3.1/HG_U133A.chip +HG_U133AAOFAV2 HG_U133AAOFAV2 /data/genomes/gsea/3.1/HG_U133AAOFAV2.chip +HG_U133A_2 HG_U133A_2 /data/genomes/gsea/3.1/HG_U133A_2.chip +HG_U133B HG_U133B /data/genomes/gsea/3.1/HG_U133B.chip +HG_U133_Plus_2 HG_U133_Plus_2 /data/genomes/gsea/3.1/HG_U133_Plus_2.chip +HG_U95Av2 HG_U95Av2 /data/genomes/gsea/3.1/HG_U95Av2.chip +HG_U95B HG_U95B /data/genomes/gsea/3.1/HG_U95B.chip +HG_U95C HG_U95C /data/genomes/gsea/3.1/HG_U95C.chip +HG_U95D HG_U95D /data/genomes/gsea/3.1/HG_U95D.chip +HG_U95E HG_U95E /data/genomes/gsea/3.1/HG_U95E.chip +HPCGGCompugenAnnotations HPCGGCompugenAnnotations /data/genomes/gsea/3.1/HPCGGCompugenAnnotations.chip +HT_HG_U133A HT_HG_U133A /data/genomes/gsea/3.1/HT_HG_U133A.chip +HT_HG_U133A_EA HT_HG_U133A_EA /data/genomes/gsea/3.1/HT_HG_U133A_EA.chip +HU6800 HU6800 /data/genomes/gsea/3.1/HU6800.chip +Hu35KsubA Hu35KsubA /data/genomes/gsea/3.1/Hu35KsubA.chip +Hu35KsubB Hu35KsubB /data/genomes/gsea/3.1/Hu35KsubB.chip +Hu35KsubC Hu35KsubC /data/genomes/gsea/3.1/Hu35KsubC.chip +Hu35KsubD Hu35KsubD /data/genomes/gsea/3.1/Hu35KsubD.chip +HuEx_1_0_STv2 HuEx_1_0_STv2 /data/genomes/gsea/3.1/HuEx_1_0_STv2.chip +HuGene_1_0_st HuGene_1_0_st /data/genomes/gsea/3.1/HuGene_1_0_st.chip +Illimina_Mus6_v1_1 Illimina_Mus6_v1_1 /data/genomes/gsea/3.1/Illimina_Mus6_v1_1.chip +Illimuna_Mus6_v1 Illimuna_Mus6_v1 /data/genomes/gsea/3.1/Illimuna_Mus6_v1.chip +Illumina_Hum6_v1 Illumina_Hum6_v1 /data/genomes/gsea/3.1/Illumina_Hum6_v1.chip +Illumina_Hum6_v2 Illumina_Hum6_v2 /data/genomes/gsea/3.1/Illumina_Hum6_v2.chip +Illumina_HumRef8_v1 Illumina_HumRef8_v1 /data/genomes/gsea/3.1/Illumina_HumRef8_v1.chip +Illumina_HumRef8_v2 Illumina_HumRef8_v2 /data/genomes/gsea/3.1/Illumina_HumRef8_v2.chip +Illumina_Human Illumina_Human /data/genomes/gsea/3.1/Illumina_Human.chip +Illumina_Mus6_v2 Illumina_Mus6_v2 /data/genomes/gsea/3.1/Illumina_Mus6_v2.chip +Illumina_MusRef8_v1 Illumina_MusRef8_v1 /data/genomes/gsea/3.1/Illumina_MusRef8_v1.chip +Illumina_MusRef8_v1_1 Illumina_MusRef8_v1_1 /data/genomes/gsea/3.1/Illumina_MusRef8_v1_1.chip +Illumina_RatRef12_v1 Illumina_RatRef12_v1 /data/genomes/gsea/3.1/Illumina_RatRef12_v1.chip +MG_U74Av2 MG_U74Av2 /data/genomes/gsea/3.1/MG_U74Av2.chip +MG_U74Bv2 MG_U74Bv2 /data/genomes/gsea/3.1/MG_U74Bv2.chip +MG_U74Cv2 MG_U74Cv2 /data/genomes/gsea/3.1/MG_U74Cv2.chip +MOE430A MOE430A /data/genomes/gsea/3.1/MOE430A.chip +MOE430B MOE430B /data/genomes/gsea/3.1/MOE430B.chip +MWG_Human_30K_A MWG_Human_30K_A /data/genomes/gsea/3.1/MWG_Human_30K_A.chip +MWG_Human_30K_B MWG_Human_30K_B /data/genomes/gsea/3.1/MWG_Human_30K_B.chip +MoEx_1_0_st MoEx_1_0_st /data/genomes/gsea/3.1/MoEx_1_0_st.chip +MoGene_1_0_st MoGene_1_0_st /data/genomes/gsea/3.1/MoGene_1_0_st.chip +MoGene_1_1_st MoGene_1_1_st /data/genomes/gsea/3.1/MoGene_1_1_st.chip +Mouse430A_2 Mouse430A_2 /data/genomes/gsea/3.1/Mouse430A_2.chip +Mouse430_2 Mouse430_2 /data/genomes/gsea/3.1/Mouse430_2.chip +MouseHsHomolSym MouseHsHomolSym /data/genomes/gsea/3.1/mm_hs_homologene.chip +MouseHsJaxSym MouseHsJaxSym /data/genomes/gsea/3.1/mm_hs_orthologues.chip +Mu11KsubA Mu11KsubA /data/genomes/gsea/3.1/Mu11KsubA.chip +Mu11KsubB Mu11KsubB /data/genomes/gsea/3.1/Mu11KsubB.chip +Mu19KsubA Mu19KsubA /data/genomes/gsea/3.1/Mu19KsubA.chip +Mu19KsubB Mu19KsubB /data/genomes/gsea/3.1/Mu19KsubB.chip +Mu19KsubC Mu19KsubC /data/genomes/gsea/3.1/Mu19KsubC.chip +NIA15k NIA15k /data/genomes/gsea/3.1/NIA15k.chip +Netherland_cancer_institute_operon_human_35k Netherland_cancer_institute_operon_human_35k /data/genomes/gsea/3.1/Netherland_cancer_institute_operon_human_35k.chip +Netherland_cancer_institute_operon_mouse_FOOk Netherland_cancer_institute_operon_mouse_FOOk /data/genomes/gsea/3.1/Netherland_cancer_institute_operon_mouse_FOOk.chip +OPERON_HUMANv2 OPERON_HUMANv2 /data/genomes/gsea/3.1/OPERON_HUMANv2.chip +OPERON_HUMANv3 OPERON_HUMANv3 /data/genomes/gsea/3.1/OPERON_HUMANv3.chip +Ortholog_SEQ_ACCESSION_MOUSE Ortholog_SEQ_ACCESSION_MOUSE /data/genomes/gsea/3.1/Ortholog_SEQ_ACCESSION_MOUSE.chip +RAE230A RAE230A /data/genomes/gsea/3.1/RAE230A.chip +RAE230B RAE230B /data/genomes/gsea/3.1/RAE230B.chip +RG_U34A RG_U34A /data/genomes/gsea/3.1/RG_U34A.chip +RG_U34B RG_U34B /data/genomes/gsea/3.1/RG_U34B.chip +RG_U34C RG_U34C /data/genomes/gsea/3.1/RG_U34C.chip +RN_U34 RN_U34 /data/genomes/gsea/3.1/RN_U34.chip +RT_U34 RT_U34 /data/genomes/gsea/3.1/RT_U34.chip +RZPD_Human_Ensembl1.1 RZPD_Human_Ensembl1.1 /data/genomes/gsea/3.1/RZPD_Human_Ensembl1.1.chip +RZPD_Human_ORF_Clones_Gateway RZPD_Human_ORF_Clones_Gateway /data/genomes/gsea/3.1/RZPD_Human_ORF_Clones_Gateway.chip +RZPD_Human_Unigene3.1 RZPD_Human_Unigene3.1 /data/genomes/gsea/3.1/RZPD_Human_Unigene3.1.chip +Rat230_2 Rat230_2 /data/genomes/gsea/3.1/Rat230_2.chip +RatHsHomolSym RatHsHomolSym /data/genomes/gsea/3.1/rn_hs_homologene.chip +RefSeq_NP_Human RefSeq_NP_Human /data/genomes/gsea/3.1/RefSeq_NP_Human.chip +RefSeq_NP_Mouse RefSeq_NP_Mouse /data/genomes/gsea/3.1/RefSeq_NP_Mouse.chip +RefSeq_NP_Rat RefSeq_NP_Rat /data/genomes/gsea/3.1/RefSeq_NP_Rat.chip +RefSeq_human RefSeq_human /data/genomes/gsea/3.1/RefSeq_human.chip +Research_Genetics Research_Genetics /data/genomes/gsea/3.1/Research_Genetics.chip +Rosetta Rosetta /data/genomes/gsea/3.1/Rosetta.chip +Rosetta50K Rosetta50K /data/genomes/gsea/3.1/Rosetta50K.chip +Seq_Accession Seq_Accession /data/genomes/gsea/3.1/Seq_Accession.chip +Stanford Stanford /data/genomes/gsea/3.1/Stanford.chip +Stanford_Source_Accessions Stanford_Source_Accessions /data/genomes/gsea/3.1/Stanford_Source_Accessions.chip +TIGR_31K_Human_Set TIGR_31K_Human_Set /data/genomes/gsea/3.1/TIGR_31K_Human_Set.chip +TIGR_40K_Human_Set TIGR_40K_Human_Set /data/genomes/gsea/3.1/TIGR_40K_Human_Set.chip +TRC TRC /data/genomes/gsea/3.1/TRC.chip +TRC_DB TRC_DB /data/genomes/gsea/3.1/TRC_DB.chip +TRC_DB_v1 TRC_DB_v1 /data/genomes/gsea/3.1/TRC_DB_v1.chip +U133_X3P U133_X3P /data/genomes/gsea/3.1/U133_X3P.chip +UCLA_NIH_33K UCLA_NIH_33K /data/genomes/gsea/3.1/UCLA_NIH_33K.chip +Zebrafish Zebrafish /data/genomes/gsea/3.1/Zebrafish.chip +ilmn_HUMANREF_8_V3_0_R1_11282963_A_WGDASL ilmn_HUMANREF_8_V3_0_R1_11282963_A_WGDASL /data/genomes/gsea/3.1/ilmn_HUMANREF_8_V3_0_R1_11282963_A_WGDASL.chip +ilmn_HumanHT_12_V3_0_R3_11283641_A ilmn_HumanHT_12_V3_0_R3_11283641_A /data/genomes/gsea/3.1/ilmn_HumanHT_12_V3_0_R3_11283641_A.chip +ilmn_HumanHT_12_V4_0_R1_15002873_B ilmn_HumanHT_12_V4_0_R1_15002873_B /data/genomes/gsea/3.1/ilmn_HumanHT_12_V4_0_R1_15002873_B.chip +ilmn_HumanRef_8_V2_0_R4_11223162_A ilmn_HumanRef_8_V2_0_R4_11223162_A /data/genomes/gsea/3.1/ilmn_HumanRef_8_V2_0_R4_11223162_A.chip +ilmn_HumanRef_8_V3_0_R3_11282963_A ilmn_HumanRef_8_V3_0_R3_11282963_A /data/genomes/gsea/3.1/ilmn_HumanRef_8_V3_0_R3_11282963_A.chip +ilmn_HumanWG_6_V2_0_R4_11223189_A ilmn_HumanWG_6_V2_0_R4_11223189_A /data/genomes/gsea/3.1/ilmn_HumanWG_6_V2_0_R4_11223189_A.chip +ilmn_HumanWG_6_V3_0_R3_11282955_A ilmn_HumanWG_6_V3_0_R3_11282955_A /data/genomes/gsea/3.1/ilmn_HumanWG_6_V3_0_R3_11282955_A.chip +ilmn_MouseRef_8_V1_1_R4_11234312_A ilmn_MouseRef_8_V1_1_R4_11234312_A /data/genomes/gsea/3.1/ilmn_MouseRef_8_V1_1_R4_11234312_A.chip +ilmn_MouseRef_8_V2_0_R3_11278551_A ilmn_MouseRef_8_V2_0_R3_11278551_A /data/genomes/gsea/3.1/ilmn_MouseRef_8_V2_0_R3_11278551_A.chip +ilmn_MouseWG_6_V1_1_R4_11234304_A ilmn_MouseWG_6_V1_1_R4_11234304_A /data/genomes/gsea/3.1/ilmn_MouseWG_6_V1_1_R4_11234304_A.chip +ilmn_MouseWG_6_V2_0_R3_11278593_A ilmn_MouseWG_6_V2_0_R3_11278593_A /data/genomes/gsea/3.1/ilmn_MouseWG_6_V2_0_R3_11278593_A.chip +ilmn_RatRef_12_V1_0_R5_11222119_A ilmn_RatRef_12_V1_0_R5_11222119_A /data/genomes/gsea/3.1/ilmn_RatRef_12_V1_0_R5_11222119_A.chip +labonweb_human labonweb_human /data/genomes/gsea/3.1/labonweb_human.chip +lymphochip lymphochip /data/genomes/gsea/3.1/lymphochip.chip diff -r aef0ecc00683 -r f24f30eae819 gseaGMT_3.1_symbols.loc.sample --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gseaGMT_3.1_symbols.loc.sample Sun Jun 09 02:00:04 2013 -0400 @@ -0,0 +1,23 @@ +# generated by makeGSEAloc.py at 24/01/2013 10:43:16 from path = /data/genomes/gsea/3.1 +# loc file for the rgGSEA tool pointing to files from ftp://gseaftp.broadinstitute.org/pub/gsea/annotations +# name value path +c1.all.v3.1.symbols c1.all.v3.1.symbols /data/genomes/gsea/3.1/c1.all.v3.1.symbols.gmt +c2.all.v3.1.symbols c2.all.v3.1.symbols /data/genomes/gsea/3.1/c2.all.v3.1.symbols.gmt +c2.cgp.v3.1.symbols c2.cgp.v3.1.symbols /data/genomes/gsea/3.1/c2.cgp.v3.1.symbols.gmt +c2.cp.biocarta.v3.1.symbols c2.cp.biocarta.v3.1.symbols /data/genomes/gsea/3.1/c2.cp.biocarta.v3.1.symbols.gmt +c2.cp.kegg.v3.1.symbols c2.cp.kegg.v3.1.symbols /data/genomes/gsea/3.1/c2.cp.kegg.v3.1.symbols.gmt +c2.cp.reactome.v3.1.symbols c2.cp.reactome.v3.1.symbols /data/genomes/gsea/3.1/c2.cp.reactome.v3.1.symbols.gmt +c2.cp.v3.1.symbols c2.cp.v3.1.symbols /data/genomes/gsea/3.1/c2.cp.v3.1.symbols.gmt +c3.all.v3.1.symbols c3.all.v3.1.symbols /data/genomes/gsea/3.1/c3.all.v3.1.symbols.gmt +c3.mir.v3.1.symbols c3.mir.v3.1.symbols /data/genomes/gsea/3.1/c3.mir.v3.1.symbols.gmt +c3.tft.v3.1.symbols c3.tft.v3.1.symbols /data/genomes/gsea/3.1/c3.tft.v3.1.symbols.gmt +c4.all.v3.1.symbols c4.all.v3.1.symbols /data/genomes/gsea/3.1/c4.all.v3.1.symbols.gmt +c4.cgn.v3.1.symbols c4.cgn.v3.1.symbols /data/genomes/gsea/3.1/c4.cgn.v3.1.symbols.gmt +c4.cm.v3.1.symbols c4.cm.v3.1.symbols /data/genomes/gsea/3.1/c4.cm.v3.1.symbols.gmt +c5.all.v3.1.symbols c5.all.v3.1.symbols /data/genomes/gsea/3.1/c5.all.v3.1.symbols.gmt +c5.bp.v3.1.symbols c5.bp.v3.1.symbols /data/genomes/gsea/3.1/c5.bp.v3.1.symbols.gmt +c5.cc.v3.1.symbols c5.cc.v3.1.symbols /data/genomes/gsea/3.1/c5.cc.v3.1.symbols.gmt +c5.mf.v3.1.symbols c5.mf.v3.1.symbols /data/genomes/gsea/3.1/c5.mf.v3.1.symbols.gmt +c6.all.v3.1.symbols c6.all.v3.1.symbols /data/genomes/gsea/3.1/c6.all.v3.1.symbols.gmt +msigdb.v3.1.symbols msigdb.v3.1.symbols /data/genomes/gsea/3.1/msigdb.v3.1.symbols.gmt + diff -r aef0ecc00683 -r f24f30eae819 rgGSEA.py --- a/rgGSEA.py Sun Jun 09 01:30:20 2013 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,495 +0,0 @@ -""" -April 2013 -eeesh GSEA does NOT respect the mode flag! - -Now realise that the creation of the input rank file for gsea needs to take the lowest p value for duplicate -feature names. To make Ish's life easier, remove duplicate gene ids from any gene set to stop GSEA from -barfing. - -October 14 2012 -Amazingly long time to figure out that GSEA fails with useless error message if any filename contains a dash "-" -eesh. - -Added history .gmt source - requires passing a faked name to gsea -Wrapper for GSEA http://www.broadinstitute.org/gsea/index.jsp -Started Feb 22 -Copyright 2012 Ross Lazarus -All rights reserved -Licensed under the LGPL - -called eg as - -#!/bin/sh -GALAXY_LIB="/data/extended/galaxy/lib" -if [ "$GALAXY_LIB" != "None" ]; then - if [ -n "$PYTHONPATH" ]; then - PYTHONPATH="$GALAXY_LIB:$PYTHONPATH" - else - PYTHONPATH="$GALAXY_LIB" - fi - export PYTHONPATH -fi - -cd /data/extended/galaxy/database/job_working_directory/027/27311 -python /data/extended/galaxy/tools/rgenetics/rgGSEA.py --input_tab "/data/extended/galaxy/database/files/033/dataset_33806.dat" --adjpvalcol "5" --signcol "2" ---idcol "1" --outhtml "/data/extended/galaxy/database/files/034/dataset_34455.dat" --input_name "actaearly-Controlearly-actalate-Controllate_topTable.xls" ---setMax "500" --setMin "15" --nPerm "1000" --plotTop "20" ---gsea_jar "/data/extended/galaxy/tool-data/shared/jars/gsea2-2.0.12.jar" ---output_dir "/data/extended/galaxy/database/job_working_directory/027/27311/dataset_34455_files" --mode "Max_probe" - --title " actaearly-Controlearly-actalate-Controllate_interpro_GSEA" --builtin_gmt "/data/genomes/gsea/3.1/IPR_DOMAIN.gmt" - - -""" -import optparse -import tempfile -import os -import sys -import subprocess -import time -import shutil -import glob -import math -import re - -KEEPSELECTION = False # detailed records for selection of multiple probes - -def timenow(): - """return current time as a string - """ - return time.strftime('%d/%m/%Y %H:%M:%S', time.localtime(time.time())) - - - -def fix_subdir(adir,destdir): - """ Galaxy wants everything in the same files_dir - if os.path.exists(adir): - for (d,dirs,files) in os.path.walk(adir): - for f in files: - sauce = os.path.join(d,f) - shutil.copy(sauce,destdir) - """ - - def fixAffycrap(apath=''): - """class='richTable'>RUNNING ESCORE ENRICHMENT1 - LBR1113 - 0.194-0.1065No2 - GGPS143090.014-0.4328 - No - """ - html = [] - try: - html = open(apath,'r').readlines() - except: - return html - for i,row in enumerate(html): - row = re.sub('https\:\/\/www.affymetrix.com\/LinkServlet\?probeset=',"http://www.genecards.org/index.php?path=/Search/keyword/",row) - html[i] = row - return html - - cleanup = False - if os.path.exists(adir): - flist = os.listdir(adir) # get all files created - for f in flist: - apath = os.path.join(adir,f) - dest = os.path.join(destdir,f) - if not os.path.isdir(apath): - if os.path.splitext(f)[1].lower() == '.html': - html = fixAffycrap(apath) - fixed = open(apath,'w') - fixed.write('\n'.join(html)) - fixed.write('\n') - fixed.close() - if not os.path.isfile(dest): - shutil.copy(apath,dest) - else: - fix_subdir(apath,destdir) - if cleanup: - try: - shutil.rmtree(path=adir,ignore_errors=True) - except: - pass - - - -def getFileString(fpath, outpath): - """ - format a nice file size string - """ - size = '' - fp = os.path.join(outpath, fpath) - s = fpath - if os.path.isfile(fp): - n = float(os.path.getsize(fp)) - if n > 2**20: - size = ' (%1.1f MB)' % (n/2**20) - elif n > 2**10: - size = ' (%1.1f KB)' % (n/2**10) - elif n > 0: - size = ' (%d B)' % (int(n)) - s = '%s %s' % (fpath, size) - return s - -class gsea_wrapper: - """ - GSEA java desktop client has a CL interface. CL can be gleaned by clicking the 'command line' button after setting up an analysis - We don't want gsea to do the analysis but it can read .rnk files containing rows of identifiers and an evidence weight such as the signed t statistic from limma for differential expression -(vgalaxy)rlazarus@iaas1:~/public_html/idle_illumina_analysis$ cat gseaHumanREFSEQ.sh -#!/bin/bash -for RNK in `ls *.rnk` -do -DIRNAME=${RNK%.*} -echo $DIRNAME -qsub -cwd -b y java -Xmx4096m -cp /data/app/bin/gsea2-2.07.jar xtools.gsea.GseaPreranked -gmx ../msigdb.v3.0.symbols.gmt -collapse true -mode Max_probe -norm meandiv --nperm 1000 -rnk $RNK -scoring_scheme weighted -rpt_label $RNK -chip ../RefSeq_human.chip -include_only_symbols true -make_sets true -plot_top_x 20 -rnd_seed timestamp --set_max 500 -set_min 15 -zip_report false -out gseaout/${DIRNAME} -gui false -done - """ - - def __init__(self,myName=None,opts=None): - """ setup cl for gsea - """ - self.idcol = 0 - self.signcol = 0 - self.adjpvalcol = 0 - self.progname=myName - self.opts = opts - remove_duplicates=True - if not os.path.isdir(opts.output_dir): - try: - os.makedirs(opts.output_dir) - except: - print >> sys.stderr,'##Error: GSEA wrapper unable to create or find output directory %s. Stopping' % (opts.output_dir) - sys.exit(1) - fakeGMT = re.sub('[^a-zA-Z0-9_]+', '', opts.input_name) # gives a more useful title for the GSEA report - fakeGMT = os.path.join(opts.output_dir,fakeGMT) - fakeRanks = '%s.rnk' % fakeGMT - if not fakeGMT.endswith('.gmt'): - fakeGMT = '%s.gmt' % fakeGMT - if opts.builtin_gmt and opts.history_gmt: - newfile = open(fakeGMT,'w') - subprocess.call(['cat',opts.builtin_gmt,opts.history_gmt],stdout=newfile) - newfile.close() - elif opts.history_gmt: - subprocess.call(['cp',opts.history_gmt,fakeGMT]) - else: - subprocess.call(['cp',opts.builtin_gmt,fakeGMT]) - # remove dupes from each gene set - gmt = open(fakeGMT,'r').readlines() - gmt = [x for x in gmt if len(x.split('\t')) > 3] - ugmt = [] - for i,row in enumerate(gmt): - rows = row.rstrip().split('\t') - gmtname = rows[0] - gmtcomment = rows[1] - glist = list(set(rows[2:])) - newgmt = [gmtname,gmtcomment] - newgmt += glist - ugmt.append('\t'.join(newgmt)) - gmt = open(fakeGMT,'w') - gmt.write('\n'.join(ugmt)) - gmt.write('\n') - gmt.close() - if opts.input_ranks: - infname = opts.input_ranks - rdat = open(opts.input_ranks,'r').readlines() # suck in and remove blank ids that cause gsea to barf rml april 10 2012 - rdat = [x.rstrip().split('\t') for x in rdat[1:]] # ignore head - dat = [[x[0],x[1],x[1]] for x in rdat] - # fake same structure as input tabular file - try: - pvals = [float(x[1]) for x in dat] - signs = [float(x[1]) for x in dat] - except: - print >> sys.stderr, '## error converting floating point - cannot process this input' - sys.exit(99) - else: # read tabular - self.idcol = int(opts.idcol) - 1 - self.signcol = int(opts.signcol) - 1 - self.adjpvalcol = int(opts.adjpvalcol) - 1 - maxcol = max(self.idcol,self.signcol,self.adjpvalcol) - infname = opts.input_tab - indat = open(opts.input_tab,'r').readlines() - dat = [x.rstrip().split('\t') for x in indat[1:]] - dat = [x for x in dat if len(x) > maxcol] - dat = [[x[self.idcol],x[self.adjpvalcol],x[self.signcol]] for x in dat] # reduce to rank form - pvals = [float(x[1]) for x in dat] - outofrange = [x for x in pvals if ((x < 0.0) or (x > 1.0))] - assert len(outofrange) == 0, '## p values outside 0-1 encountered - was that the right column for adjusted p value?' - signs = [float(x[2]) for x in dat] - outofrange = [i for i,x in enumerate(signs) if (not x) and (dat[i][self.signcol] <> '0')] - bad = [dat[x][2] for x in outofrange] - assert len(outofrange) == 0, '## null numeric values encountered for sign - was that the right column? %s' % bad - ids = [x[0] for x in dat] - res = [] - self.comments = [] - useme = [] - for i,row in enumerate(dat): - if row[1].upper() != 'NA' and row[2].upper() != 'NA' and row[0] != '' : - useme.append(i) - lost = len(dat) - len(useme) - if lost <> 0: - newdat = [dat[x] for x in useme] - del dat - dat = newdat - print >> sys.stdout, '## %d lost - NA values or null id' % lost - if remove_duplicates: - uids = list(set(ids)) # complex procedure to get min pval for each unique id - if len(uids) <> len(ids): # dupes - deal with mode - print >> sys.stdout,'## Dealing with %d uids in %d ids' % (len(uids),len(ids)) - ures = {} - for i,id in enumerate(ids): - p = pvals[i] - ures.setdefault(id,[]) - ures[id].append((p,signs[i])) - for id in uids: - tlist = ures[id] - tp = [x[0] for x in tlist] - ts = [x[1] for x in tlist] - if len(tp) == 1: - p = tp[0] - sign = ts[0] - else: - if opts.mode == "Max_probe": - p = min(tp) - sign = ts[tp.index(p)] - else: # guess median - too bad if even count - tp.sort() - ltp = len(tp) - ind = ltp/2 # yes, this is wrong for evens but what if sign wobbles? - if ltp % 2 == 1: # odd - ind += 1 # take the median - p = tp[ind] - sign = ts[ind] - if KEEPSELECTION: - self.comments.append('## for id=%s, got tp=%s, ts=%s, chose p=%f, sign =%f'\ - % (id,str(tp),str(ts),p,sign)) - if opts.input_ranks: # must be a rank file - res.append((id,'%f' % p)) - else: - if p == 0.0: - p = 1e-99 - try: - lp = -math.log10(p) # large positive if low p value - except ValueError: - lp = 0.0 - if sign < 0: - lp = -lp # if negative, swap p to negative - res.append((id,'%f' % lp)) - else: # no duplicates - for i,row in enumerate(dat): - (id,p,sign) = (row[0],float(row[1]),float(row[2])) - if opts.input_ranks: # must be a rank file - res.append((id,'%f' % p)) - else: - if p == 0.0: - p = 1e-99 - try: - lp = -math.log10(p) # large positive if low p value - except ValueError: - lp = 0.0 - if sign < 0: - lp = -lp # if negative, swap p to negative - res.append((id,'%f' % lp)) - else: - for i,row in enumerate(dat): - (id,p,sign) = (row[0],float(row[1]),float(row[2])) - if opts.input_ranks: # must be a rank file - res.append((id,'%f' % p)) - else: - if p == 0.0: - p = 1e-99 - try: - lp = -math.log10(p) # large positive if low p value - except ValueError: - lp = 0.0 - if sign < 0: - lp = -lp # if negative, swap p to negative - res.append((id,'%f' % lp)) - len1 = len(ids) - ranks = res # [x for x in res if len(x[0]) > 0 and len(x[1].strip()) > 0 and x[0].upper() <> 'NA'] - # rd = dict(zip([x[0] for x in res],res)) # this does something mysterious to dupes - probably keeps last one - # ranks = rd.values() - len2 = len(ranks) - delta = len1 - len2 - if delta <> 0: - print >> sys.stdout,'NOTE: %d of %d rank input file %s rows deleted - dup, null or NA IDs, pvals or logFCs' % (delta,len1,infname) - ranks = [(float(x[1]),x) for x in ranks] # decorate - ranks.sort() - ranks.reverse() - ranks = [x[1] for x in ranks] # undecorate - if opts.chip == '': # if mouse - need HUGO - ranks = [[x[0].upper(),x[1]] for x in ranks] - print >> sys.stdout, '## Fixed any lower case - now have',','.join([x[0] for x in ranks[:5]]) - ranks = ['\t'.join(x) for x in ranks] - if len(ranks) < 2: - print >> sys.stderr,'Input %s has 1 or less rows with two tab delimited fields - please check the tool documentation' % infname - sys.exit(1) - rclean = open(fakeRanks,'w') - rclean.write('contig\tscore\n') - rclean.write('\n'.join(ranks)) - rclean.write('\n') - rclean.close() - cl = [] - a = cl.append - a('java -Xmx6G -cp') - a(opts.gsea_jar) - a('xtools.gsea.GseaPreranked') - a('-gmx %s' % fakeGMT) # ensure .gmt extension as required by GSEA - gene sets to use - a('-gui false') # use preranked file mode and no gui - a('-make_sets true -rnd_seed timestamp') # more things from the GUI command line display - a('-norm meandiv -zip_report false -scoring_scheme weighted') # ? need to set these? - a('-rnk %s' % fakeRanks) # input ranks file symbol (the chip file is the crosswalk for ids in first column) - a('-out %s' % opts.output_dir) - a('-set_max %s' % opts.setMax) - a('-set_min %s' % opts.setMin) - a('-mode %s' % opts.mode) - if opts.chip > '': - #a('-chip %s -collapse true -include_only_symbols true' % opts.chip) - a('-chip %s -collapse true' % opts.chip) - else: - a("-collapse false") # needed if no chip - a('-nperm %s' % opts.nPerm) - a('-rpt_label %s' % opts.title) - a('-plot_top_x %s' % opts.plotTop) - self.cl = cl - self.comments.append('## GSEA command line:') - self.comments.append(' '.join(self.cl)) - self.fakeRanks = fakeRanks - self.fakeGMT = fakeGMT - - def grepIds(self): - """ - """ - found = [] - allids = open(self.opts.input_ranks,'r').readlines() - allids = [x.strip().split() for x in allids] - allids = [x[0] for x in allids] # list of ids - gmtpath = os.path.split(self.opts.fakeGMT)[0] # get path to all chip - - def run(self): - """ - - """ - tlog = os.path.join(self.opts.output_dir,"gsea_runner.log") - sto = open(tlog,'w') - x = subprocess.Popen(' '.join(self.cl),shell=True,stdout=sto,stderr=sto,cwd=self.opts.output_dir) - retval = x.wait() - sto.close() - d = glob.glob(os.path.join(self.opts.output_dir,'%s*' % self.opts.title)) - if len(d) > 0: - fix_subdir(d[0],self.opts.output_dir) - htmlfname = os.path.join(self.opts.output_dir,'index.html') - try: - html = open(htmlfname,'r').readlines() - html = [x.strip() for x in html if len(x.strip()) > 0] - if len(self.comments) > 0: - s = ['
']
-                s += self.comments
-                s.append('
') - try: - i = html.index('