extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/plugins/BarcodeFilesForDemuxExporter.java

Code
Comments
Other
Rev Date Author Line
2422 13 May 14 nicklas 1 package net.sf.basedb.reggie.plugins;
2422 13 May 14 nicklas 2
2422 13 May 14 nicklas 3 import java.io.Writer;
5586 29 Aug 19 nicklas 4 import java.util.ArrayList;
2422 13 May 14 nicklas 5 import java.util.Collection;
5586 29 Aug 19 nicklas 6 import java.util.HashMap;
2422 13 May 14 nicklas 7 import java.util.HashSet;
2422 13 May 14 nicklas 8 import java.util.List;
5586 29 Aug 19 nicklas 9 import java.util.Map;
2422 13 May 14 nicklas 10 import java.util.Set;
3753 17 Feb 16 nicklas 11 import java.util.regex.Pattern;
2422 13 May 14 nicklas 12
2422 13 May 14 nicklas 13 import net.sf.basedb.core.DbControl;
2422 13 May 14 nicklas 14 import net.sf.basedb.core.DerivedBioAssay;
2422 13 May 14 nicklas 15 import net.sf.basedb.core.Extract;
3753 17 Feb 16 nicklas 16 import net.sf.basedb.core.InvalidDataException;
2422 13 May 14 nicklas 17 import net.sf.basedb.core.Tag;
2422 13 May 14 nicklas 18 import net.sf.basedb.reggie.dao.Annotationtype;
4306 17 Jan 17 nicklas 19 import net.sf.basedb.reggie.grid.ScriptUtil;
2422 13 May 14 nicklas 20 import net.sf.basedb.util.export.TableWriter;
2422 13 May 14 nicklas 21
2422 13 May 14 nicklas 22 /**
2422 13 May 14 nicklas 23   Export that is used to export barcode files for the demux step using Picard. 
2422 13 May 14 nicklas 24   The demux is made in two steps, each one requiring a different barcode file:
2422 13 May 14 nicklas 25   
3716 21 Jan 16 nicklas 26   1. ExtractIlluminaBarcodes: {@link #exportBarcodesFile(DbControl, DerivedBioAssay, int, List, Writer)}
3716 21 Jan 16 nicklas 27   2. IlluminaBasecallsToFastq: {@link #exportMultiplexFile(DbControl, DerivedBioAssay, int, List, Writer, String)}
2422 13 May 14 nicklas 28   
2422 13 May 14 nicklas 29   In both cases one file per lane is needed, so after creating the exporter,
2422 13 May 14 nicklas 30   a loop over the lanes is typically made with calls to the above methods.
2422 13 May 14 nicklas 31   
2422 13 May 14 nicklas 32   @author nicklas
3716 21 Jan 16 nicklas 33   @since 2.16, 4.1
2422 13 May 14 nicklas 34 */
2422 13 May 14 nicklas 35 public class BarcodeFilesForDemuxExporter 
2422 13 May 14 nicklas 36 {
2422 13 May 14 nicklas 37
2422 13 May 14 nicklas 38   private final Collection<Tag> allBarcodes;
3753 17 Feb 16 nicklas 39   private final Pattern barcodePattern;
3753 17 Feb 16 nicklas 40   private final Pattern tagNamePattern;
5586 29 Aug 19 nicklas 41   // Merged name --> List of demuxed file prefixes
5587 30 Aug 19 nicklas 42   private final Map<String, List<OutputFileData>> outputFilePrefixes;
2422 13 May 14 nicklas 43   
2422 13 May 14 nicklas 44   /**
2422 13 May 14 nicklas 45     Create a new exporter.
2422 13 May 14 nicklas 46     @param allBarcodes Optional, if given, the exported file contain dummy lines for all
2422 13 May 14 nicklas 47       borcodes that are not used by libraries on the flow cell
2422 13 May 14 nicklas 48   */
5865 12 Mar 20 nicklas 49   public BarcodeFilesForDemuxExporter(DbControl dc, Collection<Tag> allBarcodes)
2422 13 May 14 nicklas 50   {
2422 13 May 14 nicklas 51     this.allBarcodes = allBarcodes;
3753 17 Feb 16 nicklas 52     this.barcodePattern = Pattern.compile("[ACGT]+");
5585 26 Aug 19 nicklas 53     this.tagNamePattern = Pattern.compile("[a-zA-Z0-9\\-\\.]+");
5586 29 Aug 19 nicklas 54     this.outputFilePrefixes = new HashMap<>();
2422 13 May 14 nicklas 55   }
2422 13 May 14 nicklas 56   
2422 13 May 14 nicklas 57   /**
2429 16 May 14 nicklas 58     Export a file hat can be used as the "BARCODE_FILE" parameter
2429 16 May 14 nicklas 59     when running "picard ExtractIlluminaBarcodes" for the given lane.
2422 13 May 14 nicklas 60   */
5865 12 Mar 20 nicklas 61   public void exportBarcodesFile(DbControl dc, String flowCellId, int laneNo, int numBarcodeReads, boolean useReverseComplementOnBarcode2, List<DerivedBioAssay> mergedInLane, List<Extract> libIgnore, Writer out)
2422 13 May 14 nicklas 62   {
2422 13 May 14 nicklas 63     TableWriter tw = createTableWriter(out);
2422 13 May 14 nicklas 64     // Write column headers
5489 12 Jun 19 nicklas 65     Object[] headers = numBarcodeReads == 2 ? 
5489 12 Jun 19 nicklas 66       new Object[] { "barcode_sequence_1", "barcode_sequence_2", "barcode_name", "library_name" } : 
5489 12 Jun 19 nicklas 67       new Object[] { "barcode_sequence_1", "barcode_name", "library_name" };
2422 13 May 14 nicklas 68     tw.tablePrintData(headers);
2422 13 May 14 nicklas 69     
2422 13 May 14 nicklas 70     Object[] info = new Object[headers.length];
2422 13 May 14 nicklas 71     Set<Tag> barcodes = new HashSet<Tag>(allBarcodes);
5865 12 Mar 20 nicklas 72     Set<String> barcodeSets = new HashSet<>();
3716 21 Jan 16 nicklas 73     for (DerivedBioAssay merged : mergedInLane)
2422 13 May 14 nicklas 74     {
3716 21 Jan 16 nicklas 75       Extract lib = merged.getExtract();
3716 21 Jan 16 nicklas 76       Tag barcode = lib.getTag();
5865 12 Mar 20 nicklas 77       String barcodeSet = (String)Annotationtype.BARCODE_SET.getAnnotationValue(dc, barcode);
5865 12 Mar 20 nicklas 78       if (barcodeSet != null) barcodeSets.add(barcodeSet);
5489 12 Jun 19 nicklas 79
5489 12 Jun 19 nicklas 80       int index = 0;
5489 12 Jun 19 nicklas 81       info[index++] = checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE);
5489 12 Jun 19 nicklas 82       if (numBarcodeReads == 2)
5489 12 Jun 19 nicklas 83       {
5592 06 Sep 19 nicklas 84         info[index++] = reverseComplement(checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE_2), useReverseComplementOnBarcode2);
5489 12 Jun 19 nicklas 85       }
5489 12 Jun 19 nicklas 86       info[index++] = checkValidBarcodeName(barcode);
5489 12 Jun 19 nicklas 87       info[index++] = merged.getName();
3716 21 Jan 16 nicklas 88       tw.tablePrintData(info);
3716 21 Jan 16 nicklas 89       barcodes.remove(barcode);
2422 13 May 14 nicklas 90     }
2422 13 May 14 nicklas 91     
4985 28 Sep 18 nicklas 92     // Libraries that should be ignored
4985 28 Sep 18 nicklas 93     // We can't put them under the "UNUSED" tag since the demux will warn if it finds
4985 28 Sep 18 nicklas 94     // a large number of unused barcodes!
4985 28 Sep 18 nicklas 95     for (Extract lib : libIgnore)
4985 28 Sep 18 nicklas 96     {
4985 28 Sep 18 nicklas 97       Tag barcode = lib.getTag();
4985 28 Sep 18 nicklas 98       if (barcodes.remove(barcode))
4985 28 Sep 18 nicklas 99       {
5489 12 Jun 19 nicklas 100         int index = 0;
5489 12 Jun 19 nicklas 101         info[index++] = checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE);
5489 12 Jun 19 nicklas 102         if (numBarcodeReads == 2)
5489 12 Jun 19 nicklas 103         {
5592 06 Sep 19 nicklas 104           info[index++] = reverseComplement(checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE_2), useReverseComplementOnBarcode2);
5489 12 Jun 19 nicklas 105         }
5489 12 Jun 19 nicklas 106         info[index++] = checkValidBarcodeName(barcode);
5489 12 Jun 19 nicklas 107         info[index++] = "IGNORED";
4985 28 Sep 18 nicklas 108         tw.tablePrintData(info);
4985 28 Sep 18 nicklas 109       }
4985 28 Sep 18 nicklas 110     }
4985 28 Sep 18 nicklas 111     
2422 13 May 14 nicklas 112     // Remaining barcodes not used by any of the libs
2422 13 May 14 nicklas 113     if (barcodes.size() > 0)
2422 13 May 14 nicklas 114     {
2422 13 May 14 nicklas 115       for (Tag barcode : barcodes)
2422 13 May 14 nicklas 116       {
5865 12 Mar 20 nicklas 117         if (barcodeSets.size() > 0)
5865 12 Mar 20 nicklas 118         {
5865 12 Mar 20 nicklas 119           // Check if the barocode is in the same "barcode set" as used by the libraries
5865 12 Mar 20 nicklas 120           // If not, we do not need to output "UNUSED" entry for it
5865 12 Mar 20 nicklas 121           String barcodeSet = (String)Annotationtype.BARCODE_SET.getAnnotationValue(dc, barcode);
5865 12 Mar 20 nicklas 122           if (!barcodeSets.contains(barcodeSet)) continue;
5865 12 Mar 20 nicklas 123         }
5489 12 Jun 19 nicklas 124         int index = 0;
5489 12 Jun 19 nicklas 125         info[index++] = checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE);
5489 12 Jun 19 nicklas 126         if (numBarcodeReads == 2)
5489 12 Jun 19 nicklas 127         {
5592 06 Sep 19 nicklas 128           info[index++] = reverseComplement(checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE_2), useReverseComplementOnBarcode2);
5489 12 Jun 19 nicklas 129         }
5489 12 Jun 19 nicklas 130         info[index++] = checkValidBarcodeName(barcode);
5489 12 Jun 19 nicklas 131         info[index++] = "UNUSED";
2422 13 May 14 nicklas 132         tw.tablePrintData(info);
2422 13 May 14 nicklas 133       }
2422 13 May 14 nicklas 134     }
2422 13 May 14 nicklas 135     
2422 13 May 14 nicklas 136     tw.flush();    
2422 13 May 14 nicklas 137   }
2422 13 May 14 nicklas 138   
2422 13 May 14 nicklas 139   /**
2429 16 May 14 nicklas 140     Export a file hat can be used as the "MULTIPLEX_PARAMS" parameter
2429 16 May 14 nicklas 141     when running "picard IlluminaBasecallsToFastq" for the given lane.
2422 13 May 14 nicklas 142   */
5865 12 Mar 20 nicklas 143   public void exportMultiplexFile(DbControl dc, String flowcellId, int laneNo, int numBarcodeReads, boolean useReverseComplementOnBarcode2, List<DerivedBioAssay> mergedInLane, List<Extract> libIgnore, Writer out, String outFolder)
2422 13 May 14 nicklas 144   {
2422 13 May 14 nicklas 145     TableWriter tw = createTableWriter(out);
2422 13 May 14 nicklas 146     // Write column headers
5489 12 Jun 19 nicklas 147     Object[] headers = numBarcodeReads == 2 ?
5489 12 Jun 19 nicklas 148       new Object[] { "OUTPUT_PREFIX", "BARCODE_1", "BARCODE_2" } :
5489 12 Jun 19 nicklas 149       new Object[] { "OUTPUT_PREFIX", "BARCODE_1" };
2422 13 May 14 nicklas 150     tw.tablePrintData(headers);
2422 13 May 14 nicklas 151     
2422 13 May 14 nicklas 152     Object[] info = new Object[headers.length];
3716 21 Jan 16 nicklas 153     for (DerivedBioAssay merged : mergedInLane)
2422 13 May 14 nicklas 154     {
3716 21 Jan 16 nicklas 155       Extract lib = merged.getExtract();
3716 21 Jan 16 nicklas 156       Tag barcode = lib.getTag();
2422 13 May 14 nicklas 157       
5586 29 Aug 19 nicklas 158       String prefix = ScriptUtil.checkValidFilename(merged.getName()+"_" + flowcellId + "." + laneNo);
5587 30 Aug 19 nicklas 159       List<OutputFileData> filePrefixes = outputFilePrefixes.get(merged.getName());
5586 29 Aug 19 nicklas 160       if (filePrefixes == null)
5586 29 Aug 19 nicklas 161       {
5586 29 Aug 19 nicklas 162         filePrefixes = new ArrayList<>();
5586 29 Aug 19 nicklas 163         outputFilePrefixes.put(merged.getName(), filePrefixes);
5586 29 Aug 19 nicklas 164       }
5586 29 Aug 19 nicklas 165       
5587 30 Aug 19 nicklas 166       OutputFileData data = new OutputFileData(merged.getName(), flowcellId, laneNo, prefix);
5587 30 Aug 19 nicklas 167       filePrefixes.add(data);
5586 29 Aug 19 nicklas 168       
5586 29 Aug 19 nicklas 169       info[0] = outFolder+"/"+prefix;
5587 30 Aug 19 nicklas 170       data.barcode1 = checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE);
5587 30 Aug 19 nicklas 171       info[1] = data.barcode1;
5489 12 Jun 19 nicklas 172       if (numBarcodeReads == 2)
5489 12 Jun 19 nicklas 173       {
5592 06 Sep 19 nicklas 174         data.barcode2 = reverseComplement(checkValidBarcodeSequence(dc, barcode, Annotationtype.BARCODE_SEQUENCE_2), useReverseComplementOnBarcode2);
5587 30 Aug 19 nicklas 175         info[2] = data.barcode2;
5489 12 Jun 19 nicklas 176       }
3716 21 Jan 16 nicklas 177       tw.tablePrintData(info);
2422 13 May 14 nicklas 178     }
2422 13 May 14 nicklas 179     
2422 13 May 14 nicklas 180     tw.flush();
2422 13 May 14 nicklas 181   }
5586 29 Aug 19 nicklas 182   
5586 29 Aug 19 nicklas 183   
5586 29 Aug 19 nicklas 184   /**
5586 29 Aug 19 nicklas 185     Get the output file prefixes that has been used for the given merged bioassay.
5586 29 Aug 19 nicklas 186     Typically there is one prefix for each flowCell/lane combination where the
5586 29 Aug 19 nicklas 187     merged bioassay is present.
5586 29 Aug 19 nicklas 188   */
5587 30 Aug 19 nicklas 189   public List<OutputFileData> getOutputFilePrefixes(DerivedBioAssay merged)
5586 29 Aug 19 nicklas 190   {
5586 29 Aug 19 nicklas 191     return outputFilePrefixes.get(merged.getName());
5586 29 Aug 19 nicklas 192   }
5586 29 Aug 19 nicklas 193   
2422 13 May 14 nicklas 194   private TableWriter createTableWriter(Writer out)
2422 13 May 14 nicklas 195   {
2422 13 May 14 nicklas 196     TableWriter tw = new TableWriter(out);
2422 13 May 14 nicklas 197     tw.setDataSeparator("\t"); // Tab-separated file
2422 13 May 14 nicklas 198     tw.setNullValue("");
2422 13 May 14 nicklas 199     // Get rid of "bad" characters (comma, tab, newline, etc.)
2422 13 May 14 nicklas 200     tw.setEncoder(new CsvEncoderDecoder());
2422 13 May 14 nicklas 201     return tw;
2422 13 May 14 nicklas 202   }
2422 13 May 14 nicklas 203   
3753 17 Feb 16 nicklas 204   /**
3753 17 Feb 16 nicklas 205     Checks that the barcode sequence contains only letters ACGT.
3753 17 Feb 16 nicklas 206   */
5489 12 Jun 19 nicklas 207   private String checkValidBarcodeSequence(DbControl dc, Tag barcode, Annotationtype barcodeType)
3753 17 Feb 16 nicklas 208   {
5489 12 Jun 19 nicklas 209     String sequence = (String)barcodeType.getAnnotationValue(dc, barcode);
3753 17 Feb 16 nicklas 210     if (sequence == null)
3753 17 Feb 16 nicklas 211     {
5489 12 Jun 19 nicklas 212       throw new InvalidDataException("No " + barcodeType.getName() + " annotation on tag '" + barcode.getName()+"'");
3753 17 Feb 16 nicklas 213     }
3753 17 Feb 16 nicklas 214     if (!barcodePattern.matcher(sequence).matches())
3753 17 Feb 16 nicklas 215     {
5489 12 Jun 19 nicklas 216       throw new InvalidDataException(barcodeType.getName() + " '"+sequence+"' is not valid on tag '" + barcode.getName()+
3753 17 Feb 16 nicklas 217           "': Only 'ACGT' is allowed.");
3753 17 Feb 16 nicklas 218     }
3753 17 Feb 16 nicklas 219     return sequence;
3753 17 Feb 16 nicklas 220   }
3753 17 Feb 16 nicklas 221   
3753 17 Feb 16 nicklas 222   /**
3753 17 Feb 16 nicklas 223      Checks that barcode name contains only letters (a-z) and numbers.
3753 17 Feb 16 nicklas 224   */
3753 17 Feb 16 nicklas 225   private String checkValidBarcodeName(Tag barcode)
3753 17 Feb 16 nicklas 226   {
3753 17 Feb 16 nicklas 227     String name = barcode.getName();
3753 17 Feb 16 nicklas 228     if (!tagNamePattern.matcher(name).matches())
3753 17 Feb 16 nicklas 229     {
5585 26 Aug 19 nicklas 230       throw new InvalidDataException("Barcode name '"+barcode.getName()+"' is not valid. Only 'a-zA-Z0-9', '-' and '.' is allowed.");
3753 17 Feb 16 nicklas 231     }
3753 17 Feb 16 nicklas 232     return name;
3753 17 Feb 16 nicklas 233   }
3753 17 Feb 16 nicklas 234   
5586 29 Aug 19 nicklas 235   /**
5592 06 Sep 19 nicklas 236     Change the given sequence to it's reverse complement. This is needed since the second barcode read is 
5592 06 Sep 19 nicklas 237     handled differently on different sequencers. 
5592 06 Sep 19 nicklas 238     
5592 06 Sep 19 nicklas 239     See https://support.illumina.com/content/dam/illumina-support/documents/documentation/system_documentation/miseq/indexed-sequencing-overview-guide-15057455-05.pdf
5592 06 Sep 19 nicklas 240     
5592 06 Sep 19 nicklas 241     Note that we have decided to store barcode sequences according to NextSeq specifications
5592 06 Sep 19 nicklas 242     so we need to reverse-complement when demuxing data from HiSeq 2000/2500, MiSeq, etc.
5592 06 Sep 19 nicklas 243   */
5592 06 Sep 19 nicklas 244   private String reverseComplement(String seq, boolean useReverseComplementOnBarcode2)
5592 06 Sep 19 nicklas 245   {
5592 06 Sep 19 nicklas 246     if (!useReverseComplementOnBarcode2) return seq;
5592 06 Sep 19 nicklas 247     StringBuilder rev = new StringBuilder();
5592 06 Sep 19 nicklas 248     for (int i = seq.length()-1; i >= 0; i--)
5592 06 Sep 19 nicklas 249     {
5592 06 Sep 19 nicklas 250       char c = seq.charAt(i);
5592 06 Sep 19 nicklas 251       switch (c)
5592 06 Sep 19 nicklas 252       {
5592 06 Sep 19 nicklas 253         case 'A': rev.append('T'); break;
5592 06 Sep 19 nicklas 254         case 'C': rev.append('G'); break;
5592 06 Sep 19 nicklas 255         case 'G': rev.append('C'); break;
5592 06 Sep 19 nicklas 256         case 'T': rev.append('A'); break;
5592 06 Sep 19 nicklas 257       }
5592 06 Sep 19 nicklas 258     }
5592 06 Sep 19 nicklas 259     return rev.toString();
5592 06 Sep 19 nicklas 260   }
5592 06 Sep 19 nicklas 261   
5592 06 Sep 19 nicklas 262   /**
5586 29 Aug 19 nicklas 263     Represents the output-file prefix for a demuxed data. The 'prefix' is the path to use for file
5586 29 Aug 19 nicklas 264     names. The other properties are useful metadata.
5586 29 Aug 19 nicklas 265   */
5587 30 Aug 19 nicklas 266   public static class OutputFileData
5586 29 Aug 19 nicklas 267   {
5586 29 Aug 19 nicklas 268     
5586 29 Aug 19 nicklas 269     public final String mergeName;
5586 29 Aug 19 nicklas 270     public final String flowCellId;
5586 29 Aug 19 nicklas 271     public final int lane;
5586 29 Aug 19 nicklas 272     public final String prefix;
5587 30 Aug 19 nicklas 273     public String barcode1;
5587 30 Aug 19 nicklas 274     public String barcode2;
5586 29 Aug 19 nicklas 275     
5587 30 Aug 19 nicklas 276     OutputFileData(String mergeName, String flowCellId, int lane, String prefix)
5586 29 Aug 19 nicklas 277     {
5586 29 Aug 19 nicklas 278       this.mergeName = mergeName;
5586 29 Aug 19 nicklas 279       this.flowCellId = flowCellId;
5586 29 Aug 19 nicklas 280       this.lane = lane;
5586 29 Aug 19 nicklas 281       this.prefix = prefix;
5586 29 Aug 19 nicklas 282     }
5586 29 Aug 19 nicklas 283     
5586 29 Aug 19 nicklas 284   }
2422 13 May 14 nicklas 285 }

Generated on Tue Nov 21 03:34:50 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)