Genome-Model-Tools-Music
view release on metacpan or search on metacpan
lib/Genome/Model/Tools/Music/Survival.pm view on Meta::CPAN
373839404142434445464748495051525354555657
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
140141142143144145146147148149150151152153154155156157158159160Qunyuan 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
295296297298299300301302303304305306307308309310311312313314
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
);
"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
12345678910111213141516171819202122232425### 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
525354555657585960616263646566676869707172
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)
( run in 0.312 second using v1.01-cache-2.11-cpan-4e96b696675 )