Genome-Model-Tools-Music

 view release on metacpan or  search on metacpan

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

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
    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

140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
 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

295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
    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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
### 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

52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
    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)