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

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

Renamed Metagenomics module to MassSequencing? module

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