|
2
|
1 #!/usr/bin/env Rscript
|
|
|
2 ##
|
|
|
3 ## formats raw clinical data from TCGA to contain a single status & time colums
|
|
|
4 ##
|
|
|
5 ## Input (required):
|
|
|
6 ## - clinical data
|
|
|
7 ## Input (optional):
|
|
|
8 ## - status & time columns: (NOT USED IN THIS SCRIPT - see comment below)
|
|
|
9 ## ideally, a better design would allow a user to specify 1 or more columns
|
|
|
10 ## to check for the status & time columns - however, due to the necessities
|
|
|
11 ## required to pre-process the TCGA clinical data, the script would not be
|
|
|
12 ## generalizeable - and for this reason, the TCGA columns are hard-coded.
|
|
|
13 ##
|
|
|
14 ## Output: a re-formatted clinical file containing 3 columns: sample-ID, status & time
|
|
|
15 ##
|
|
|
16 ## Date: August 21, 2012
|
|
|
17 ## Author: Peter Waltman
|
|
|
18 ##
|
|
|
19
|
|
|
20 ##usage, options and doc goes here
|
|
|
21 argspec <- c("format.raw.TCGA.clinical.data.R takes a clustering from ConsensusClusterPlus and clinical survival data
|
|
|
22 and generates a KM-plot, along with the log-rank p-values
|
|
|
23
|
|
|
24 Usage:
|
|
|
25 format.raw.TCGA.clinical.data.R -c <clinical.file>
|
|
|
26 Options:
|
|
|
27 -o <output file> (tab-delimited (3 col: sample_id <tab> status <tab> time))
|
|
|
28 ")
|
|
|
29 args <- commandArgs(TRUE)
|
|
|
30 if ( length( args ) == 1 && args =="--help") {
|
|
|
31 write(argspec, stderr())
|
|
|
32 q();
|
|
|
33 }
|
|
|
34
|
|
|
35 lib.load.quiet <- function( package ) {
|
|
|
36 package <- as.character(substitute(package))
|
|
|
37 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
|
|
|
38 }
|
|
|
39 lib.load.quiet(getopt)
|
|
|
40
|
|
|
41 spec <- matrix( c( "clinical.fname", "d", 1, "character",
|
|
|
42 "output.fname", "o", 2, "character"
|
|
|
43 ),
|
|
|
44 ncol=4,
|
|
|
45 byrow=TRUE
|
|
|
46 )
|
|
|
47 opt <- getopt( spec=spec )
|
|
|
48 save.image( "/tmp/format.dbg.rda")
|
|
|
49
|
|
|
50 ##set some reasonable defaults for the options that are needed,
|
|
|
51 ##but were not specified.
|
|
|
52 if ( is.null(opt$output.fname ) ) { opt$output.fname <-file.path( getwd(), "formated.TCGA.clinical.data" ) }
|
|
|
53
|
|
|
54 ##orig.clinical.data <- read.delim( opt$clinical.fname, as.is=TRUE, row.names=1 )
|
|
|
55 orig.clinical.data <- read.delim( opt$clinical.fname, as.is=TRUE )
|
|
|
56 orig.clinical.data <- unique( orig.clinical.data )
|
|
|
57 rownames( orig.clinical.data ) <- orig.clinical.data[,1]
|
|
|
58 orig.clinical.data <- orig.clinical.data[, -1 ]
|
|
|
59
|
|
|
60 ## ugh, some TCGA data sets have all NAs in the "days_to_..." columns
|
|
|
61 if ( "days_to_last_known_alive" %in% colnames( orig.clinical.data ) ) {
|
|
|
62 time.cols <- c( "days_to_death", "days_to_last_followup", "days_to_last_known_alive" )
|
|
|
63 } else {
|
|
|
64 time.cols <- c( "days_to_death", "days_to_last_followup" )
|
|
|
65 }
|
|
|
66 good.samps <- ! apply( orig.clinical.data[, time.cols ], 1, function(x) all( is.na(x) ) | all( x <= 0, na.rm=T ) )
|
|
|
67
|
|
|
68 orig.clinical.data <- orig.clinical.data[ good.samps, ]
|
|
|
69
|
|
|
70 if ( is.null(opt$status.column ) ) {
|
|
|
71 status.colname <- "vital_status"
|
|
|
72 if ( status.colname %in% colnames( orig.clinical.data ) ) {
|
|
|
73 opt$status.column <- which( colnames( orig.clinical.data ) %in% status.colname )
|
|
|
74 clinical.data <- orig.clinical.data[ , opt$status.column ]
|
|
|
75 }
|
|
|
76 else {
|
|
|
77 status.colname <- "days_to_death"
|
|
|
78 if ( status.colname %in% colnames( orig.clinical.data ) ) {
|
|
|
79 opt$status.column <- which( colnames( orig.clinical.data ) %in% status.colname )
|
|
|
80 clinical.data <- orig.clinical.data[ , opt$status.column ]
|
|
|
81 }
|
|
|
82 else {
|
|
|
83 stop( "can't find a valid entry with status info - have tried vital_status & days_to_death\n" )
|
|
|
84 }
|
|
|
85 }
|
|
|
86 clinical.data <- as.numeric( ! grepl( "(LIVING|Not)", clinical.data ) )
|
|
|
87 }
|
|
|
88 if ( is.null(opt$time.column ) ) {
|
|
|
89 time.colname <- "CDE.clinical_time"
|
|
|
90
|
|
|
91 if ( time.colname %in% colnames( orig.clinical.data ) ) {
|
|
|
92 opt$time.column <- which( colnames( orig.clinical.data ) %in% time.colname )
|
|
|
93 clinical.data <- cbind( clinical.data,
|
|
|
94 as.numeric( orig.clinical.data[, opt$time.column ] ) )
|
|
|
95 }
|
|
|
96 else {
|
|
|
97 dec.mat <- matrix( NA,
|
|
|
98 nc=length( time.cols ),
|
|
|
99 nr=nrow( orig.clinical.data ),
|
|
|
100 dimnames=list( rownames( orig.clinical.data ),
|
|
|
101 time.cols )
|
|
|
102 )
|
|
|
103 for ( cname in colnames( dec.mat ) ) {
|
|
|
104 if ( cname %in% colnames( orig.clinical.data ) ) {
|
|
|
105 dec.mat[, cname ] <- as.numeric( orig.clinical.data[, cname ] )
|
|
|
106 }
|
|
|
107 }
|
|
|
108
|
|
|
109
|
|
|
110
|
|
|
111 if ( "days_to_last_known_alive" %in% colnames( orig.clinical.data ) ) {
|
|
|
112
|
|
|
113 opt$time.column <- sapply( 1:length( clinical.data ),
|
|
|
114 function(i) {
|
|
|
115 if ( clinical.data[i] ) {
|
|
|
116 ## this is a deceased sample
|
|
|
117 return( ifelse( ( !is.na( dec.mat[ i, "days_to_death" ] ) ),
|
|
|
118 dec.mat[ i, "days_to_death" ],
|
|
|
119 ifelse( ( !is.na( dec.mat[ i, "days_to_last_known_alive" ] ) ),
|
|
|
120 dec.mat[ i, "days_to_last_known_alive" ],
|
|
|
121 dec.mat[ i, "days_to_last_followup" ] ) ) )
|
|
|
122
|
|
|
123 }
|
|
|
124 else {
|
|
|
125 return( max( dec.mat[ i, c( "days_to_last_followup","days_to_last_known_alive") ], na.rm=T ) )
|
|
|
126 }
|
|
|
127 }
|
|
|
128 )
|
|
|
129 } else {
|
|
|
130 opt$time.column <- sapply( 1:length( clinical.data ),
|
|
|
131 function(i) {
|
|
|
132 if ( clinical.data[i] ) {
|
|
|
133 ## this is a deceased sample
|
|
|
134 return( ifelse( ( !is.na( dec.mat[ i, "days_to_death" ] ) ),
|
|
|
135 dec.mat[ i, "days_to_death" ],
|
|
|
136 dec.mat[ i, "days_to_last_followup" ] ) )
|
|
|
137
|
|
|
138 }
|
|
|
139 else {
|
|
|
140 return( max( dec.mat[ i, c( "days_to_last_followup") ], na.rm=T ) )
|
|
|
141 }
|
|
|
142 }
|
|
|
143 )
|
|
|
144 }
|
|
|
145
|
|
|
146
|
|
|
147 clinical.data <- cbind( clinical.data,
|
|
|
148 as.numeric( opt$time.column ) )
|
|
|
149 }
|
|
|
150 }
|
|
|
151
|
|
|
152 clinical.data <- as.data.frame( clinical.data )
|
|
|
153 colnames( clinical.data ) <- c( "status", "time" )
|
|
|
154 rownames( clinical.data ) <- rownames( orig.clinical.data )
|
|
|
155
|
|
|
156
|
|
|
157 ## check to make sure that the id's are sync'd correctly
|
|
|
158 ## the default format is to use hyphens to separate the elt's of the name
|
|
|
159 ## and to only use the 1st 3 elements of the name
|
|
|
160 ## so we check to see if they're using something else as separators and/or using more than 3 elts
|
|
|
161 reformat.ids <- function( ids ) {
|
|
|
162
|
|
|
163 if ( grepl( "TCGA\\.", ids[1] ) ) {
|
|
|
164 ids <- sapply( strsplit( ids, "\\." ), function(x) paste( x[1:3], collapse="-" ) )
|
|
|
165 } else {
|
|
|
166 ## do this just in case there's more than 3 elements to the names
|
|
|
167 if ( grepl( "TCGA-", ids[1] ) ) {
|
|
|
168 ids <- sapply( strsplit( ids, "-" ), function(x) paste( x[1:min( c(3,length(x) ) )], collapse="-" ) )
|
|
|
169 }
|
|
|
170 }
|
|
|
171 return( ids )
|
|
|
172 }
|
|
|
173
|
|
|
174
|
|
|
175 new.samp.ids <- reformat.ids( rownames( clinical.data ) )
|
|
|
176 if ( any( duplicated( new.samp.ids ) ) ) {
|
|
|
177 ## in some cases, we have duplicate sample ids in the raw data after we truncate to
|
|
|
178 ## the 1st 3 elts in the barcode, so just simplify the data
|
|
|
179 uniqs <- ! duplicated( new.samp.ids )
|
|
|
180 clinical.data <- clinical.data[ uniqs, ]
|
|
|
181 new.samp.ids <- new.samp.ids[ uniqs ]
|
|
|
182 }
|
|
|
183
|
|
|
184 rownames( clinical.data ) <- new.samp.ids
|
|
|
185 write.table( clinical.data, opt$output.fname, sep="\t", quote=FALSE, col.names=NA )
|
