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

Last change on this file since 68 was 68, checked in by robert@…, 8 years ago

Added extra stacktrace information

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