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

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

Implemented importing of classifications

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