source: trunk/grails-app/services/nl/tno/massSequencing/FastaService.groovy @ 52

Last change on this file since 52 was 52, checked in by robert@…, 8 years ago
  • Updated a bug on production with org.apache.tools.zip
  • Added tooltips to all icons
File size: 33.6 KB
Line 
1package nl.tno.massSequencing
2
3import java.io.BufferedWriter;
4import java.io.File;
5import java.io.Writer;
6import java.util.ArrayList;
7import java.util.List;
8import org.codehaus.groovy.grails.commons.ConfigurationHolder
9import java.util.zip.*
10
11class FastaService {
12        def fileService
13        def fuzzySearchService
14        def sampleExcelService
15        def excelService
16
17        static transactional = true
18
19        /**
20         * Parses uploaded files and checks them for FASTA and QUAL files
21         * @param filenames             List of filenames currently existing in the upload directory
22         * @param onProgress    Closure to execute when progress indicators should be updated.
23         *                                              Has 4 parameters: numFilesProcessed, numBytesProcessed that indicate the number
24         *                                              of files and bytes that have been processed in total. totalFiles, totalBytes indicate the change
25         *                                              in total number of files and bytes (e.g. should take 1 if a new file is added to the list)
26         * @param directory             Directory to move the files to
27         * @return                              Structure with information about the parsed files. The 'success' files are
28         *                                              moved to the given directory
29         *
30         * [
31         *              success: [
32         *                      [filename: 'abc.fasta', type: FASTA, numSequences: 190]
33         *                      [filename: 'cde.fasta', type: FASTA, numSequences: 140]
34         *                      [filename: 'abc.qual', type: QUAL, numSequences: 190, avgQuality: 38]
35         *                      [filename: 'cde.qual', type: QUAL, numSequences: 140, avgQuality: 29]
36         *              ],
37         *              failure: [
38         *                      [filename: 'testing.xls', type: 'unknown', message: 'Type not recognized']
39         *              ]
40         * ]
41         *
42         */
43        def parseFiles( ArrayList filenames, Closure onProgress, File directory = null ) {
44                if( filenames.size() == 0 ) {
45                        return [ success: [], failure: [] ];
46                }
47
48                if( !directory ) {
49                        directory = fileService.absolutePath( ConfigurationHolder.config.massSequencing.fileDir )
50                }
51
52                def success = [];
53                def failure = [];
54
55                long filesProcessed = 0;
56                long bytesProcessed = 0;
57
58                // Loop through all filenames
59                for( int i = 0; i < filenames.size(); i++ ) {
60                        def filename = filenames[ i ];
61                       
62                        if( fileService.isZipFile( filename ) ) {
63                                // ZIP files are extracted and appended to the filenames list.
64                                def newfiles = fileService.extractZipFile( filename, { files, bytes, totalFiles, totalBytes ->
65                                        filesProcessed += files;
66                                        bytesProcessed += bytes;
67
68                                        onProgress( filesProcessed, bytesProcessed, totalFiles, totalBytes );
69                                } );
70                                if( newfiles ) {
71                                        newfiles.each {
72                                                filenames.add( it );
73                                        }
74                                       
75                                        // Delete the zip file itself
76                                        fileService.delete( filename );
77                                } else {
78                                        failure << [ filename: filename, originalfilename: fileService.originalFilename( filename ), type: "zip", message: 'Zip file could not be extracted' ];
79                                }
80                        } else {
81                                def file = fileService.get( filename );
82                                String filetype = fileService.determineFileType( file );
83
84                                if( !fileTypeValid( filetype ) ) {
85                                        // If files are not valid for parsing, delete them and return a message to the user
86                                        fileService.delete(filename);
87                                        failure << [ filename: filename, originalfilename: fileService.originalFilename( filename ), type: filetype, message: 'File type not accepted' ];
88                                } else {
89                                        try {
90                                                def contents = parseFile( file, filetype, { files, bytes ->
91                                                        filesProcessed += files;
92                                                        bytesProcessed += bytes;
93
94                                                        onProgress( filesProcessed, bytesProcessed, 0, 0 );
95                                                } );
96                                               
97                                                contents.filename = file.getName();
98                                                contents.originalfilename = fileService.originalFilename( contents.filename )
99
100                                                if( contents.success ) {
101                                                        success << contents;
102                                                } else {
103                                                        fileService.delete(filename);
104                                                        failure << contents;
105                                                }
106                                        } catch( Exception e ) {
107                                                // If anything fails during parsing, return an error message
108                                                fileService.delete(filename);
109                                                failure << [ filename: filename, originalfilename: fileService.originalFilename( filename ), type: filetype, message: e.getMessage() ];
110                                        }
111                                }
112                        }
113                }
114
115                return [ 'success': success, 'failure': failure ];
116        }
117
118        /**
119         * Matches uploaded fasta and qual files and combines them with the samples they probably belong to
120         * @param parsedFiles   Parsed files
121         * [
122         *              [filename: 'abc.fasta', type: FASTA, numSequences: 190]
123         *              [filename: 'cde.fasta', type: FASTA, numSequences: 140]
124         *              [filename: 'abc.qual', type: QUAL, numSequences: 190, avgQuality: 38]
125         *              [filename: 'cde.qual', type: QUAL, numSequences: 140, avgQuality: 29]
126         *              [filename: 'match.xls', type: EXCEL, matches: [ [ filename: 'abc.fasta', basename: 'abc', sample: 's1' ] ]
127         * ]
128         * @param samples               AssaySample objects to which the files should be matched.
129         * @return                              Structure with information about the matching.
130         * [
131         *              [
132         *                      fasta:  [filename: 'abc.fasta', type: FASTA, numSequences: 190],
133         *                      qual:   [filename: 'abc.qual', type: QUAL, numSequences: 190, avgQuality: 38],
134         *                      feasibleQuals: [
135         *                              [filename: 'abc.qual', type: QUAL, numSequences: 190, avgQuality: 38],
136         *                              [filename: 'def.qual', type: QUAL, numSequences: 190, avgQuality: 21]
137         *                      ]
138         *                      sample: AssaySample object     
139         */
140        def matchFiles( def parsedFiles, def samples ) {
141                def fastas = parsedFiles.findAll { it.type == "fasta" }
142                def quals = parsedFiles.findAll { it.type == "qual" }
143                def excels = parsedFiles.findAll { it.type == "excel" }
144               
145                samples = samples.toList()
146               
147                // Collect matches from all files
148                def matches = [];
149                excels.each { m ->
150                        if( m.matches ) {
151                                m.matches.each { matches << it }
152                        }
153                }
154
155                def files = [];
156
157                fastas.each { fastaFile ->
158                        // Remove extension
159                        def matchWith = fastaFile.originalfilename.substring( 0, fastaFile.originalfilename.lastIndexOf( '.' ) )
160
161                        // Determine feasible quals (based on number of sequences )
162                        def feasibleQuals = quals.findAll { it.numSequences == fastaFile.numSequences }
163
164                        // Best matching qual file
165                        def qualIdx = fuzzySearchService.mostSimilarWithIndex( matchWith + '.qual', feasibleQuals.originalfilename );
166
167                        def qual = null
168                        if( qualIdx != null )
169                                qual = feasibleQuals[ qualIdx ];
170
171                        // Best matching sample
172                        def assaySample = null
173                        def checked = true;
174                        if( matches ) {
175                                // Match with files from excelsheet
176                               
177                                // First find the best matching filename in the list of matches.
178                                def sampleNameIdx = fuzzySearchService.mostSimilarWithIndex( matchWith, matches*.basename, 0.8 );
179                               
180                                // If one is found, use the sample name associated with it to do the matching with samples
181                                if( sampleNameIdx != null ) {
182                                        matchWith = matches[ sampleNameIdx ].sample;
183                                } else {
184                                        // If no match is found in the excel sheet, this sample is unchecked by default
185                                        checked = false;
186                                }
187                        }
188                       
189                        // Match on filenames
190                        def sampleIdx = fuzzySearchService.mostSimilarWithIndex( matchWith, samples.sample.name );
191                        if( sampleIdx != null ) {
192                                assaySample = samples[ sampleIdx ];
193                        }
194
195                        files << [
196                                                fasta: fastaFile,
197                                                feasibleQuals: feasibleQuals,
198                                                qual: qual,
199                                                assaySample: assaySample,
200                                                checked: checked
201                                        ]
202                }
203
204                return files;
205
206
207        }
208
209        /**
210         * Determines whether a file can be processed.
211         * @param filetype      Filetype of the file
212         * @see determineFileType()
213         * @return
214         */
215        protected boolean fileTypeValid( String filetype ) {
216                switch( filetype ) {
217                        case "fasta":
218                        case "qual":
219                        case "excel":
220                                return true;
221                        default:
222                                return false;
223                }
224        }
225
226        /**
227         * Parses the given file
228         * @param file                  File to parse
229         * @param filetype              Type of the given file
230         * @param onProgress    Closure to execute when progress indicators should be updated.
231         *                                              Has 2 parameters: numFilesProcessed and numBytesProcessed that indicate the number
232         *                                              of files and bytes that have been processed in this file (so the first parameter should
233         *                                              only be 1 when the file is finished)
234         *
235         * @return                              List structure. Examples:
236         *
237         *   [ success: true, filename: 'abc.fasta', type: 'fasta', numSequences: 200 ]
238         *   [ success: true, filename: 'abc.qual', type: 'qual', numSequences: 200, avgQuality: 36 ]
239         *   [ success: false, filename: 'abc.txt', type: 'txt', message: 'Filetype could not be parsed.' ]
240         */
241        protected def parseFile( File file, String filetype, Closure onProgress ) {
242                switch( filetype ) {
243                        case "fasta":
244                                return parseFasta( file, onProgress );
245                        case "qual":
246                                return parseQual( file, onProgress );
247                        case "excel":
248                                return parseExcelMatch( file, onProgress );
249                        default:
250                                onProgress( 1, file.length() );
251                                return [ success: false, type: filetype, message: 'Filetype could not be parsed.' ]
252                }
253        }
254
255        /**
256         * Parses the given FASTA file
257         * @param file                  File to parse
258         * @param onProgress    Closure to execute when progress indicators should be updated.
259         *                                              Has 2 parameters: numFilesProcessed and numBytesProcessed that indicate the number
260         *                                              of files and bytes that have been processed in this file (so the first parameter should
261         *                                              only be 1 when the file is finished)
262         * @return                              List structure. Examples:
263         *
264         *   [ success: true, filename: 'abc.fasta', type: 'fasta', numSequences: 200 ]
265         *   [ success: false, filename: 'def.fasta', type: 'fasta', message: 'File is not a valid FASTA file' ]
266         */
267        protected def parseFasta( File file, Closure onProgress ) {
268
269                long startTime = System.nanoTime();
270                log.trace "Start parsing FASTA " + file.getName()
271
272                // Count the number of lines, starting with '>' (and where the following line contains a character other than '>')
273                long numSequences = 0;
274                long bytesProcessed = 0;
275                boolean lookingForSequence = false;
276
277                file.eachLine { line ->
278                        if( line ) {
279                                if( !lookingForSequence && line[0] == '>' ) {
280                                        lookingForSequence = true;
281                                } else if( lookingForSequence ) {
282                                        if( line[0] != '>' ) {
283                                                numSequences++;
284                                                lookingForSequence = false;
285                                        }
286                                }
287
288
289                                // Update progress every time 1MB is processed
290                                bytesProcessed += line.size();
291                                if( bytesProcessed > 1000000 ) {
292                                        onProgress( 0, bytesProcessed );
293                                        bytesProcessed = 0;
294                                }
295                        }
296
297                }
298
299                // Update progress and say we're finished
300                onProgress( 1, bytesProcessed );
301
302                log.trace "Finished parsing FASTA " + file.getName() + ": " + ( System.nanoTime() - startTime ) / 1000000L
303
304                return [ success: true, type: "fasta", filename: file.getName(), numSequences: numSequences ];
305        }
306
307        /**
308         * Parses the given FASTA file and returns a list with sequence lengths
309         * @param file                  File to parse
310         * @return List                 [ 400: 12, 410: 15, 'unknown': 14 ]
311         */
312        protected def parseFastaLengths( File file ) {
313
314                long startTime = System.nanoTime();
315                log.trace "Start parsing FASTA " + file.getName()
316
317                // Count the number of lines, starting with '>' (and where the following line contains a character other than '>')
318                def histogram = [:]
319                def lengthPattern = ~/length=(\d+)/
320                def length
321                def lengthMatches
322                try {
323                        file.eachLine { line ->
324                                if( line ) {
325                                        if( line[0] == '>' ) {
326                                                // Comments line: find length=###
327                                                lengthMatches = ( line =~ lengthPattern );
328                                                if( lengthMatches ) {
329                                                        length = Integer.valueOf( lengthMatches[0][1] );
330                                                } else {
331                                                        length = 'unknown'
332                                                }
333
334                                                histogram[ length ] = ( histogram[length] ?: 0 ) + 1;
335                                        }
336                                }
337                        }
338                } catch( Exception e ) {
339                        log.error( "Exception while reading fasta file " + file + ": " + e.getMessage() )
340                }
341
342                log.trace "Finished parsing FASTA " + file.getName() + ": " + ( System.nanoTime() - startTime ) / 1000000L
343
344                return histogram;
345        }
346
347        /**
348         * Parses the given QUAL file
349         * @param file                  File to parse
350         * @param onProgress    Closure to execute when progress indicators should be updated.
351         *                                              Has 2 parameters: numFilesProcessed and numBytesProcessed that indicate the number
352         *                                              of files and bytes that have been processed in this file (so the first parameter should
353         *                                              only be 1 when the file is finished)
354         * @return                              List structure. Examples:
355         *
356         *   [ success: true, filename: 'abc.qual', type: 'qual', numSequences: 200, avgQuality: 31 ]
357         *   [ success: false, filename: 'def.qual', type: 'qual', message: 'File is not a valid QUAL file' ]
358         */
359        protected def parseQual( File file, Closure onProgress ) {
360                long startTime = System.nanoTime();
361                log.trace "Start parsing QUAL " + file.getName()
362
363                // Count the number of lines, starting with '>'. After we've found such a character, we continue looking for
364                // quality scores
365                long numSequences = 0;
366                long bytesProcessed = 0;
367                def quality = [ quality: 0.0, number: 0L ]
368
369                boolean lookingForFirstQualityScores = false;
370                file.eachLine { line ->
371                        if( line ) {
372                                if( !lookingForFirstQualityScores && line[0] == '>' ) {
373                                        lookingForFirstQualityScores = true;
374                                } else if( lookingForFirstQualityScores ) {
375                                        if( line[0] != '>' ) {
376                                                numSequences++;
377                                                lookingForFirstQualityScores = false;
378
379                                                // Don't compute average quality because it takes too much time
380                                                //quality = updateQuality( quality, line );
381                                        }
382                                } else {
383                                        // Don't compute average quality because it takes too much time
384                                        //quality = updateQuality( quality, line );
385                                }
386
387                                // Update progress every time 1MB is processed
388                                bytesProcessed += line.size();
389                                if( bytesProcessed > 1000000 ) {
390                                        onProgress( 0, bytesProcessed );
391                                        bytesProcessed = 0;
392                                }
393                        }
394
395                }
396
397                // Update progress and say we're finished
398                onProgress( 1, bytesProcessed );
399
400                log.trace "Finished parsing QUAL " + file.getName() + ": " + ( System.nanoTime() - startTime ) / 1000000L
401
402                return [ success: true, type: "qual", filename: file.getName(), numSequences: numSequences, avgQuality: quality.quality ];
403        }
404
405        /**
406         * Parses a given excel file with a match between filenames and samples
407         * @param file                  File to parse
408         * @param onProgress    Closure to execute when progress indicators should be updated.
409         *                                              Has 2 parameters: numFilesProcessed and numBytesProcessed that indicate the number
410         *                                              of files and bytes that have been processed in this file (so the first parameter should
411         *                                              only be 1 when the file is finished)
412         * @return                              List structure. The matches array contains an array of matches between filenames and sample(name)s. 
413         *                                              The extension for all files are removed in the 'basename' parameter, in order to improve matching. 
414         *                                              Examples:
415         *
416         *   [ success: true, filename: 'abc.xls', type: 'excel', matches: [ [ filename: 's1.qual', basename: 's1', sample: 'sample a' ], [ filename: 's9.fna', basename: 's9', sample: 'sample b' ] ]
417         *   [ success: false, filename: 'def.xls', type: 'excel', message: 'File is not a valid XLS file' ]
418         */
419        protected def parseExcelMatch( File file, Closure onProgress ) {
420                long startTime = System.nanoTime();
421                log.trace "Start parsing XLS " + file.getName()
422
423                def matches = []
424
425                // Read excel file
426                def wb;
427                try {
428                        wb = excelService.open( file );
429                } catch( Exception e ) {
430                        // If an exception occurs, the file can't be opened. Return the error message
431                        return [ success: false, type: "excel", filename: file.getName(), message: "Excel file could not be opened or parsed."]
432                }
433
434                // Read all data into an array, and the header in a separate array
435                def header = excelService.readRow( wb, 0, 0 );
436                def data = excelService.readData( wb, 0, 1 );
437
438                // Check whether (at least) 2 columns are present
439                if( header.size() < 2 ) {
440                        return [ success: false, type: "excel", filename: file.getName(), message: "Excel file must contain at least 2 columns, one with filenames and one with sample names"]
441                }
442               
443                // Check the headers to see whether the default columns are switched
444                def filenameColumn = 0;
445                def sampleColumn = 1;
446               
447                header.eachWithIndex { cell, i ->
448                        switch( cell?.toLowerCase() ) {
449                                case "sample":
450                                        sampleColumn = i; break;
451                                case "file":
452                                case "filename":
453                                        filenameColumn = i; break;
454                        }
455                }
456               
457                // If both filenames and samples are found in the same column (happens if only one of the headers is given)
458                // an error is given
459                if( filenameColumn == sampleColumn )
460                        return [ success: false, type: "excel", filename: file.getName(), message: "Excel file must contain 'Sample' and 'Filename' headers."]
461               
462                // Loop through the data and create a match
463                def maxColumn = Math.max( filenameColumn, sampleColumn );
464                data.each { row ->
465                        if( row.size() >= maxColumn ) {
466                                def filename = row[ filenameColumn ];
467                                def sample = row[ sampleColumn ];
468                               
469                                if( sample && sample != "null"  && filename && filename != "null" ) {
470                                        // Remove extension from the filename, but only if it is
471                                        // .fna, .fasta, .qual, .fqa, .xls or .xlsx. Otherwise useful parts of the filename would be removed.
472                                        def basename = ( filename =~ /\.(fna|fasta|qual|fqa|xls|xlsx)$/ ).replaceAll( "" );
473                                       
474                                        matches << [ filename: filename, sample: sample, basename: basename ];
475                                }
476                        }
477                }
478
479                // Update progress and say we're finished
480                onProgress( 1, file.size() );
481
482                log.trace "Finished parsing XLS " + file.getName() + ": " + ( System.nanoTime() - startTime ) / 1000000L
483
484                return [ success: true, type: "excel", filename: file.getName(), matches: matches ];
485        }
486
487        /**
488         * Parses the given line and updates the average quality based on the scores
489         * @param quality       [quality: 0.0f, number: 0L]
490         * @param line          String of integer quality scores, separated by a whitespace
491         * @return                      [quality: 0.0f, number: 0L]
492         */
493        protected def updateQuality( def quality, String line ) {
494                // Determine current average
495                List tokens = line.tokenize();
496                Long total = 0;
497
498                tokens.each {
499                        total += Integer.parseInt( it );
500                }
501
502                int numTokens = tokens.size();
503
504                // Update the given average
505                if( numTokens > 0 ) {
506                        quality.number += numTokens;
507                        quality.quality = quality.quality + ( ( total / numTokens as double ) - quality.quality ) / quality.number * numTokens;
508                }
509
510                return quality
511        }
512
513        /**
514         * Moves a fasta and qual file to their permanent location, and returns information about these files
515         *
516         * @param fastaFile                     Filename of the fasta file in the temporary directory
517         * @param qualFile                      Filename of the fasta file in the temporary directory
518         * @param processedFiles        Structure with data about uploaded files and matches
519         * @return      [ fasta: <fasta filename>, qual: <qual filename>, numSequences: <number of sequences>, avgQuality: <average quality> ]
520         */
521        public def savePermanent( String fastaFile, String qualFile, def processedFiles ) {
522                File permanentDirectory = fileService.absolutePath( ConfigurationHolder.config.massSequencing.fileDir );
523                def returnStructure = [:];
524
525                if( fileService.fileExists( fastaFile ) ) {
526                        // Lookup the original filename
527                        def fastaData = processedFiles.parsed.success.find { it.filename == fastaFile };
528                        if( fastaData ) {
529                                returnStructure.fasta = fileService.moveFileToUploadDir( fileService.get( fastaFile ), fastaData.originalfilename, permanentDirectory );
530                                returnStructure.numSequences = fastaData.numSequences;
531                        } else {
532                                throw new Exception( "Fasta file wasn't uploaded the right way. Maybe the session has expired/" );
533                        }
534                } else {
535                        // If the file doesn't exist, we can't save anything to the database.
536                        throw new Exception( "Fasta file to save doesn't exist on disk" );
537                }
538
539                if( qualFile && fileService.fileExists( qualFile ) ) {
540                        // Lookup the original filename
541                        def qualData = processedFiles.parsed.success.find { it.filename == qualFile };
542
543                        if( qualData ) {
544                                returnStructure.qual = fileService.moveFileToUploadDir( fileService.get(qualFile ), qualData.originalfilename, permanentDirectory );
545                                returnStructure.avgQuality = qualData.avgQuality
546                        } else {
547                                // Delete the uploaded fasta file, since this is a serious error
548                                fileService.delete( returnStructure.fasta, permanentDirectory );
549                                throw new Exception( "Qual file wasn't uploaded the right way. Maybe the session has expired" );
550                        }
551                } else {
552                        // If the file doesn't exist, we don't save any quality information
553                        returnStructure.qual = null
554                        returnStructure.avgQuality = 0;
555                }
556
557                return returnStructure;
558        }
559
560        /*
561         * Creates a histogram of sequence lengths for the given assaysamples
562         * @param       assaySamples    List of assaySamples to create the histogram for
563         * @return      List                    List with length-number tuples. Example:
564         *                                                      [ [400, 12], [402, 9] ]
565         */
566        public List sequenceLengthHistogram( ArrayList assaySamples ) {
567                if( !assaySamples || assaySamples.size() == 0 )
568                        return [];
569                       
570                // Determine the directory the uploaded files are stored in
571                File permanentDirectory = fileService.absolutePath( ConfigurationHolder.config.massSequencing.fileDir );
572                       
573                // Initialize the histogram
574                def histogram = [:];   
575               
576                // Loop through all samples and find the histogram values for each of them
577                assaySamples.each { assaySample ->
578                        if( assaySample.numSequences() > 0 ) {
579                                assaySample.sequenceData.each { sequenceData ->
580                                        File sequenceFile = fileService.get( sequenceData.sequenceFile, permanentDirectory );
581                                       
582                                        parseFastaLengths( sequenceFile ).each { k, v ->
583                                                histogram[ k ] = ( histogram[ k ] ?: 0 ) + v; 
584                                        };
585                                }
586                        }
587                }
588               
589                def lengthList = [];
590                histogram.each { 
591                        if( it.key != 'unknown' )
592                                lengthList << [ it.key, it.value ]
593                }
594               
595                return lengthList;
596        }
597
598        /**
599         * Exports the fasta data of a list of assaysamples
600         * @param assaySamples  Assaysamples to export
601         * @param outStream             Outputstream to send the data to
602         * @return
603         */
604        public def export( List assaySamples, OutputStream outStream, String name = null ) {
605                if( !assaySamples || assaySamples.size() == 0 )
606                        return false;
607
608                // Retrieve the filename from configuration, if none is given
609                if( !name )
610                        name = ConfigurationHolder.config.massSequencing.exportFilename
611
612                // Determine the directory the uploaded files are stored in
613                File permanentDirectory = fileService.absolutePath( ConfigurationHolder.config.massSequencing.fileDir );
614
615                // First check whether qual files should be exported or not
616                // It is only exported if qual scores are available for all sequences
617                def exportQual = ( assaySamples*.numSequences().sum() == assaySamples*.numQualScores().sum() );
618
619                // First create tags for every sample
620                def tags = [];
621
622                // Determine new tag length. Since we can use 4 characters per
623                // tag position, we only have to use 4log( #samples)
624                // The minimum number of characters used is 10, to ensure the correct working of the other
625                // programs.
626                int tagLength = Math.max( 10, Math.ceil( Math.log( assaySamples.size() ) / Math.log( 4 ) ) );
627                int tagNumber = 0;
628
629                assaySamples.each { assaySample ->
630                        if( assaySample.numSequences() > 0 ) {
631                                // Create a new tag for this assaysample
632                                def tag = createTag( tagLength, tagNumber++);
633
634                                // Save the tag for exporting
635                                tags << [       assaySampleId: assaySample.id, sampleName: assaySample.sample.name,
636                                                        assayName: assaySample.assay.name, studyName: assaySample.assay.study.name,
637                                                        forwardPrimer: assaySample.fwPrimerSeq, reversePrimer: assaySample.revPrimerSeq,
638                                                        tag: tag
639                                                ];
640                        }
641                }
642
643                // Now create zip file for fasta and qual files
644                ZipOutputStream zipFile = new ZipOutputStream( new BufferedOutputStream( outStream ) );
645                BufferedWriter zipWriter = new BufferedWriter( new OutputStreamWriter( zipFile ) );
646
647                // We have to loop twice through the sequenceData, since we can't write part of the sequence
648                // file and part of the qual files mixed. We have to write the full sequence file first.
649                try {
650                        zipFile.putNextEntry( new ZipEntry( name + ".fna" ) );
651
652                        assaySamples.each { assaySample ->
653                                if( assaySample.numSequences() > 0 ) {
654                                        def currentTag = tags.find { it.assaySampleId == assaySample.id };
655
656                                        assaySample.sequenceData.each { sequenceData ->
657                                                copyFastaFileForExport( fileService.get( sequenceData.sequenceFile, permanentDirectory ), currentTag.tag, zipWriter)
658                                        }
659                                }
660                        }
661                        zipWriter.flush();
662                        zipFile.closeEntry();
663
664                        if( exportQual ) {
665                                zipFile.putNextEntry( new ZipEntry( name + ".qual" ) );
666
667                                assaySamples.each { assaySample ->
668                                        if( assaySample.numSequences() > 0 ) {
669                                                def currentTag = tags.find { it.assaySampleId == assaySample.id };
670
671                                                assaySample.sequenceData.each { sequenceData ->
672                                                        copyQualFileForExport( fileService.get( sequenceData.qualityFile, permanentDirectory ), currentTag.tag, zipWriter)
673                                                }
674                                        }
675                                }
676
677                                zipWriter.flush();
678                                zipFile.closeEntry();
679                        }
680
681                } catch( Exception e ) {
682                        log.error "Error while writing to fastafile or qualfile: " + e.getMessage();
683                } finally {
684                        // Always close zip entry
685                        try {
686                                zipFile.closeEntry();
687                        } catch( Exception e ) {
688                                log.error "Error while closing zip entry for fasta and qual: " + e.getMessage();
689                        }
690                }
691
692                // Export a tab delimited file with tags
693                zipFile.putNextEntry( new ZipEntry( name + ".tab" ) );
694                exportTabDelimitedSampleTagFile( tags, zipWriter );
695                zipWriter.flush();
696                zipFile.closeEntry();
697
698                // Export a mothur file with tags
699                zipFile.putNextEntry( new ZipEntry( name + ".oligos" ) );
700                exportMothurSampleTagFile( tags, zipWriter );
701                zipWriter.flush();
702                zipFile.closeEntry();
703
704                // Export an excel file with information about the samples
705                zipFile.putNextEntry( new ZipEntry( name + ".xls" ) );
706                sampleExcelService.exportExcelSampleData( assaySamples, tags, zipFile );
707                zipFile.flush();
708                zipFile.closeEntry();
709
710                zipFile.close();
711        }
712
713        /**
714         * Creates an oligos file for Mothur that represents the connection between samples
715         * and the artificial tags. 
716         * 
717         * @see http://www.mothur.org/wiki/Trim.seqs#allfiles
718         * @param tags          Map with newly created tags
719         * @param zipWriter     Writer to write the data to
720         */
721        protected void exportMothurSampleTagFile( List tags, Writer zipWriter ) {
722                // Add the forward and reverse primers, as found in the assaysamples
723                // The primers are already cut off, so they are not relevant anymore
724                // For that reason, a '#' is prepended to each line.
725                def fwPrimers = tags.collect { it.forwardPrimer }.findAll { it }.unique();
726                def revPrimers = tags.collect { it.reversePrimer }.findAll { it }.unique();
727
728                fwPrimers.each { zipWriter.write( "#forward\t" + it + "\n" )    }
729                revPrimers.each { zipWriter.write( "#reverse\t" + it + "\n" ) }
730
731                // Check whether the sample names are unique. If they aren't, the assay and study names
732                // are appended to the sample name
733                def sampleNames = tags*.sampleNames;
734                if( sampleNames.unique().size() < sampleNames.size() ) {
735                        tags.each {
736                                it.uniqueName = it.sampleName + " (" + it.assayName + " / " + it.studyName + ")";
737                        }
738                } else {
739                        tags.each {
740                                it.uniqueName = it.sampleName;
741                        }
742                }
743
744                tags.each {
745                        zipWriter.write( "barcode\t" + it.tag + "\t" + it.uniqueName + "\n" );
746                }
747        }
748
749        /**
750         * Creates a tab delimited file with two columns and column headers "Sequence" and "Samplename"
751         * @param tags          Map with newly created tags
752         * @param zipWriter     Writer to write the data to
753         */
754        protected void exportTabDelimitedSampleTagFile( List tags, Writer zipWriter ) {
755                zipWriter.write( "Sequence" + "\t" + "Samplename" + "\n" );
756
757                // Check whether the sample names are unique. If they aren't, the assay and study names
758                // are appended to the sample name
759                def sampleNames = tags*.sampleNames;
760                if( sampleNames.unique().size() < sampleNames.size() ) {
761                        tags.each {
762                                it.uniqueName = it.sampleName + " (" + it.assayName + " / " + it.studyName + ")";
763                        }
764                } else {
765                        tags.each {
766                                it.uniqueName = it.sampleName;
767                        }
768                }
769
770                tags.each {
771                        zipWriter.write( it.tag + "\t" + it.uniqueName + "\n" );
772                }
773        }
774       
775        /**
776         * Deletes all sequences from a list of assaysamples
777         * @param assaySample
778         * @return
779         */
780        public int deleteSequenceData( ArrayList assaySamples ) {
781                def numFiles = 0;
782                assaySamples?.each { assaySample ->
783                        numFiles += assaySample.deleteSequenceData();
784                }
785               
786                return numFiles;
787        }
788
789        /**
790         * Creates a unique tag for the given number
791         * @param length
792         * @param tagNumber
793         * @return
794         */
795        public String createTag( int length, int tagNumber ) {
796                def chars = ["C", "A", "G", "T"];
797                def numChars = chars.size();
798
799                if( tagNumber > numChars ** length )
800                        throw new Exception( "Given tag number (" + tagNumber + ") is too large for the specified length (" + length + ")")
801
802                String tag = "";
803
804                for( def i = 0; i < length; i++ ) {
805                        int currentChar = tagNumber % numChars
806
807                        // Append the new character to the end of the tag, to ensure that the first part of the tag is
808                        // the most volatile. This way it is easy to find the end of the tag and the beginning of the real
809                        // sequence on first sight.
810                        tag = tag + chars[ currentChar ];
811
812                        tagNumber = Math.floor( tagNumber / numChars );
813                }
814
815                return tag;
816        }
817
818        /**
819         * Copies the contents of the given sequence file to the output file and prepends the tag to every sequences
820         * @param inFile        Filename of the file to be read
821         * @param tag           
822         * @param outWriter
823         * @return
824         */
825        protected boolean copyFastaFileForExport( File inFile, String tag, BufferedWriter outWriter ) {
826                // Walk through the lines in the file, starting with '>'
827                // (and where the following line contains a character other than '>')
828
829                try {
830                        BufferedReader inReader = new BufferedReader( new FileReader( inFile ) );
831
832                        String line = null
833                        String newLine = null
834                        String sequence = "";
835
836                        def lengthPattern = ~/length=(\d+)/
837                        def lengthMatches
838                        int length = 0;
839                        int tagLength = tag.size();
840
841                        while( ( line = inReader.readLine()) != null) {
842                                if( line.size() == 0 ) {
843                                        // Print the sequence we collected, before writing the empty line
844                                        printSequence( outWriter, sequence, tag );
845                                        sequence = "";
846
847                                        // Empty line
848                                        outWriter.newLine();
849                                } else if( line[ 0 ] == '>' ) {
850                                        // Print the sequence we collected, before writing the new comments tag
851                                        printSequence( outWriter, sequence, tag );
852                                        sequence = "";
853
854                                        // Comments line: replace length=### with the
855                                        // updated length, and put the line in the
856                                        lengthMatches = ( line =~ lengthPattern );
857                                        if( lengthMatches ) {
858                                                length = Integer.valueOf( lengthMatches[0][1] ) + tagLength;
859                                                newLine = lengthMatches.replaceAll( "length=" + length );
860                                        }
861
862                                        outWriter.write(newLine);
863                                        outWriter.newLine();
864                                } else {
865                                        // This is part of the sequence. We collect the whole sequence and
866                                        // determine in the end how to write it to the file
867                                        sequence += line;
868                                }
869                        }
870
871                        // Print the sequence we collected, before ending the file
872                        printSequence( outWriter, sequence, tag );
873                        sequence = "";
874
875                } catch( Exception e ) {
876                        log.error( "An error occurred while copying contents from " + inFile.getName() + ": " + e.getMessage() );
877                        return false;
878                }
879        }
880
881        /**
882         * Prints a sequence to the output file
883         * @param outWriter
884         * @param sequence
885         * @param tag
886         */
887        private void printSequence( BufferedWriter outWriter, String sequence, String tag, int maxWidth = 60 ) {
888                // If no sequence is given, also don't prepend it with the tag
889                if( sequence.size() == 0 )
890                        return
891
892                // Prepend the tag to the sequence
893                sequence = tag + sequence;
894
895                // Write the sequence with a width of maxWidth characters per line
896                while( sequence ) {
897                        if( sequence.size() > maxWidth ) {
898                                outWriter.write( sequence[0..maxWidth-1] );
899                                sequence = sequence[maxWidth..-1]
900                        } else {
901                                outWriter.write( sequence );
902                                sequence = null;
903                        }
904                        outWriter.newLine();
905                }
906        }
907
908        /**
909         * Copies the contents of the given qual file to the output file and prepends the tag quality score to every sequence.
910         * For every tag character '40' is prepended to the qual scores
911         *
912         * @param inFile        Filename of the file to be read
913         * @param tag
914         * @param outWriter
915         * @return
916         */
917        protected boolean copyQualFileForExport( File inFile, String tag, BufferedWriter outWriter ) {
918                // Walk through the lines in the file, starting with '>'
919                // (and where the following line contains a character other than '>')
920                try {
921                        BufferedReader inReader = new BufferedReader( new FileReader( inFile ) );
922
923                        String line = null
924                        String newLine = null
925                        List<Integer> qualScores = []
926
927                        def lengthPattern = ~/length=(\d+)/
928                        def lengthMatches
929                        int length = 0;
930                        int tagLength = tag.size();
931
932                        while( ( line = inReader.readLine()) != null) {
933                                if( line.size() == 0 ) {
934                                        // Print the quality scores we collected, before writing the empty line
935                                        printQualScore( outWriter, qualScores, tagLength );
936                                        qualScores = [];
937
938                                        // Empty line
939                                        outWriter.newLine();
940                                } else if( line[ 0 ] == '>' ) {
941                                        // Print the quality scores we collected, before writing the empty line
942                                        printQualScore( outWriter, qualScores, tagLength );
943                                        qualScores = [];
944
945                                        // Comments line: replace length=### with the
946                                        // updated length, and put the line in the
947                                        lengthMatches = ( line =~ lengthPattern );
948                                        if( lengthMatches ) {
949                                                length = Integer.valueOf( lengthMatches[0][1] ) + tagLength;
950                                                newLine = lengthMatches.replaceAll( "length=" + length );
951                                        }
952
953                                        outWriter.write(newLine);
954                                        outWriter.newLine();
955                                } else {
956                                        // This is part of the quality score. We collect the whole set of quality
957                                        // scores and determine in the end how to write it to the file
958                                        qualScores += line.split( " " ).collect {
959                                                if( !it.isInteger() )
960                                                        return 0;
961                                                else
962                                                        return Integer.parseInt( it );
963                                        };
964                                }
965                        }
966
967                        // Print the quality scores we collected, before ending the file
968                        printQualScore( outWriter, qualScores, tagLength );
969                        qualScores = [];
970
971                } catch( Exception e ) {
972                        log.error( "An error occurred while copying contents from " + inFile.getName() + ": " + e.getMessage() );
973                        return false;
974                }
975        }
976
977        /**
978         * Prints a sequence to the output file
979         * @param outWriter
980         * @param sequence
981         * @param tag
982         */
983        private void printQualScore( BufferedWriter outWriter, List<Integer> qualScores, int tagLength, int maxWidth = 60 ) {
984                // If no qualScores are given, also don't prepend it with the tag
985                if( qualScores.size() == 0 )
986                        return
987
988                // Prepend the tag to the sequence
989                qualScores = Collections.nCopies( tagLength, 40 ) + qualScores;
990
991                // Write the sequence with a width of maxWidth characters per line
992                while( qualScores ) {
993                        if( qualScores.size() > maxWidth ) {
994                                outWriter.write( qualScores[0..maxWidth-1].join( " " ) );
995                                qualScores = qualScores[maxWidth..-1]
996                        } else {
997                                outWriter.write( qualScores.join( " " ) );
998                                qualScores = null;
999                        }
1000                        outWriter.newLine();
1001                }
1002        }
1003
1004}
Note: See TracBrowser for help on using the repository browser.