@@ -179,12 +179,20 @@ if (length(SAMPLES) == 0) {
179179
180180# (B) ANALYSE
181181
182- # Initialize dataframe
182+ # Initialize dataframes
183183df <- data.frame (
184184 sample = character (),
185185 type = character (),
186186 value = character (),
187187 stringsAsFactors = FALSE )
188+ abundances <- data.frame (
189+ sample = character (),
190+ ID = character (),
191+ query = character (),
192+ target = character (),
193+ observed_abund = numeric (),
194+ expected_abund = numeric (),
195+ stringsAsFactors = FALSE )
188196
189197# for each sample
190198for (sample in SAMPLES ) {
@@ -238,10 +246,10 @@ for (sample in SAMPLES) {
238246 data $ observed_abund <- data $ observed_abund / sum(data $ observed_abund ,na.rm = TRUE )
239247 # Sort by expected abundance
240248 data <- data [order(data $ expected_abund , decreasing = TRUE ), ]
241- # write table
242- outfile <- paste0 (sample ," _abundances.tsv " )
243- print(paste( " write " , outfile ))
244- write.table( data , file = outfile , row.names = FALSE , col.names = TRUE , quote = FALSE , na = ' ' , sep = " \t " )
249+ # append sample data to table
250+ data $ sample <- rep (sample , nrow( data ) )
251+ data <- data [, c( " sample " , " ID " , " query " , " target " , " observed_abund " , " expected_abund " )]
252+ abundances <- rbind( abundances , data )
245253
246254 # (3) STATISTICS
247255
@@ -274,12 +282,13 @@ for (sample in SAMPLES) {
274282 invisible (dev.off())
275283
276284 # Scatter plot: Observed vs. Expected Abundance (log-log)
285+ data_matches <- data [data $ observed_abund > 0 & data $ expected_abund > 0 ,] # only on matched observed to expected
277286 outfile <- paste0(sample ," _scatter_loglog" )
278287 print(paste(" write" ,outfile ))
279288 svg(paste0(outfile ," .svg" ), width = 8 , height = 6 )
280289 plot(
281- log10(data $ expected_abund ),
282- log10(data $ observed_abund ),
290+ log10(data_matches $ expected_abund ),
291+ log10(data_matches $ observed_abund ),
283292 xlab = " log10(Expected Abundance)" ,
284293 ylab = " log10(Observed Abundance)" ,
285294 main = " Scatter Plot: Observed vs. Expected Abundance (log-log)" ,
@@ -321,6 +330,12 @@ for (sample in SAMPLES) {
321330
322331# (5) TABLES - metrics overall
323332
333+ # write table
334+ outfile <- " abundances_per-sample.tsv"
335+ abundances $ tag <- rep(tag , nrow(abundances ))
336+ print(paste(" write" ,outfile ))
337+ write.table(abundances , file = outfile , row.names = FALSE , col.names = TRUE , quote = FALSE , na = ' ' , sep = " \t " )
338+
324339# Write detailed output
325340outfile <- " performance_per-sample.tsv"
326341df $ tag <- rep(tag , nrow(df ))
0 commit comments