Genome-Model-Tools-Music

 view release on metacpan or  search on metacpan

lib/Genome/Model/Tools/Music/Survival.pm  view on Meta::CPAN

        doc => "Table of samples (y) vs. categorical clinical data category (x)",
    },
    glm_clinical_data_file => {
        is => 'Text', is_optional => 1,
        doc => "Clinical traits, mutational profiles, other mixed clinical data (See DESCRIPTION).",
    },
    phenotypes_to_include => {
        is => 'Text', is_optional => 1,
        doc => "Include only these genes and/or phenotypes in the anlaysis. (COMMA-DELIMITED)",
    },
    legend_placement => {
        is => 'Text', is_optional => 1, default => 'bottomleft',
        doc => "Choose one of 'bottomleft', 'topleft', 'topright', or 'bottomright'.",
    },
    skip_non_coding => {
        is => 'Boolean', is_optional => 1, default => 1,
        doc => "Skip non-coding mutations from the provided MAF file",
    },
    skip_silent => {
        is => 'Boolean', is_optional => 1, default => 1,
        doc => "Skip silent mutations from the provided MAF file",

lib/Genome/Model/Tools/Music/Survival.pm  view on Meta::CPAN

 Qunyuan Zhang, Ph.D.
EOS
}

sub execute {

    # parse input arguments
    my $self = shift;
    my $bam_list = $self->bam_list;
    my $genetic_data_type = $self->genetic_data_type;
    my $legend_placement = $self->legend_placement;

    # handle phenotype inclusions
    my @phenotypes_to_include;
    my @clinical_phenotypes_to_include;
    my @mutated_genes_to_include;
    if ($self->phenotypes_to_include) { @phenotypes_to_include = split /,/,$self->phenotypes_to_include; }

    # check genetic data type
    unless ($genetic_data_type =~ /^gene|variant$/i) {
        $self->error_message("Please enter either \"gene\" or \"variant\" for the --genetic-data-type parameter.");

lib/Genome/Model/Tools/Music/Survival.pm  view on Meta::CPAN

    my $output_dir = $self->output_dir . "/";
    unless (-e $output_dir) {
        $self->status_message("Creating output directory: $output_dir...");
        unless(mkdir $output_dir) {
            $self->error_message("Failed to create output directory: $!");
            return;
        }
    }

    # set up R command
    my $R_cmd = "R --slave --args < " . __FILE__ . ".R " . join( " ", $survival_data_file, $mutation_matrix, $legend_placement, $output_dir );
    print "R_cmd:\n$R_cmd\n";

    #run R command
    WIFEXITED( system $R_cmd ) or croak "Couldn't run: $R_cmd ($?)";

    return(1);
}

sub create_sample_gene_matrix_gene {

lib/Genome/Model/Tools/Music/Survival.pm.R  view on Meta::CPAN

### Survival analysis for mutation data ###

### original location of code: /gscuser/qzhang/gstat/survival/survival.R
### example input file: /gscuser/qzhang/gstat/survival/tcga.tsv

### Run it on command line like below 
### for example,   R --no-save --args < survival.R vital_status.input mut_matrix.input legend.placement output_dir & 

### clinical data /vital status input file, first three columns are sample_ID, survival_time, vital_status (0=living, 1=deceased)

######################## read input arguments

clinical.survival.data=commandArgs()[4];
mut.data=commandArgs()[5];
legend.placement=commandArgs()[6];
out.dir=commandArgs()[7];

######################## read and prepare data 

vitals = read.table(clinical.survival.data,header=T);
mut_matrix = read.table(mut.data,header=T);
x = merge(vitals,mut_matrix,by.x=1,by.y=1);
write.table(x,file=paste(out.dir,"survival_analysis_data_matrix.csv",sep="/"),quote=F,append=F,row.names=F,sep="\t")
colnames(x)[-c(1:3)]->phenos
if (class(x[,phenos])=="integer" & length(unique(x[,phenos]))<6) x[,phenos] [x[,phenos]>1]=1

lib/Genome/Model/Tools/Music/Survival.pm.R  view on Meta::CPAN

    mfit.by <- survfit(Surv(time, status == 1) ~ x1, data = loopdata)
    ## file name for plot
    bitmap(file=paste(out.dir,"/",phenotype,"_survival_plot.png",sep=""))
    ## create survival plot
    plot(mfit.by,lty=1:10,ylab="Survival Probability",xlab="Time",col=c(1:10))
    if (dim(table(x1))>1) {
        title(paste(phenotype,", P=",signif(p,3),sep=""));
    } else {
        title(paste(phenotype));
    }
    legend(x=legend.placement, legend=names(table(x1)), lty = 1:10, col=c(1:10)) 
    dev.off()

}

########################## calculate fdr

logr=logr[,-5]; 
if (length(phenos) < 2) { logr=(t(logr)); }
fdr=p.adjust(as.numeric(logr[,"p"]),"fdr")
logr=cbind(logr,fdr)