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

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