Changeset 2497

Timestamp:

Jun 5, 2014, 3:41:07 PM (9 years ago)

Author:

Nicklas Nordborg

Message:

References #547: Start Demux and Merge

Added trimmomatic to the script. This is called after the merge. Since the merge and copy to project_archive happened at the same time before the needed some bigger changes.

After the demux all FASTQ files are in the FASTQ folder (we can no longer use the COMPRESS_OUTPUTS option since trimmomatic don't like concatenated gz files).

The second step is to merge and compress the FASTQ files into 'fastq.merged' directory.

The third step runs trimmomatic which save the modified FASTQ files (compressed) in the 'fastq.trimmomatic' directory. Options for trimmomatic are configured in the reggie-ogs-hosts.xml file. Information about what has happened is written to a log file 'trimmomatic.out' in the job folder. When the job has been completed this information is parsed and number of surviving reads are stored in annotation PT_READS.

Last step is to simply copy the generated files to the project_archive.

Location:

extensions/net.sf.basedb.reggie/trunk

Files:

• reggie-ogs-hosts.xml (modified) (1 diff)
• src/net/sf/basedb/reggie/dao/Annotationtype.java (modified) (1 diff)
• src/net/sf/basedb/reggie/servlet/DemuxMergeServlet.java (modified) (16 diffs)
• src/net/sf/basedb/reggie/servlet/InstallServlet.java (modified) (2 diffs)

extensions/net.sf.basedb.reggie/trunk/reggie-ogs-hosts.xml

@@ -59 +59 @@
     </demux>
 
+    <!-- Settings for the Trimmomatic step -->
+    <trimmomatic>
+      <!-- full path to the JAR file with the Trimmomatic program -->
+      <jar-path>/usr/local/packages/trimmomatic/0.32/trimmomatic-0.32.jar</jar-path>
+      <!-- full path to the file with Illumina adapter information -->
+      <adapter-file>/usr/local/packages/trimmomatic/0.32/adapters/TruSeq3-PE-2.fa</adapter-file>
+      <!-- static options for Trimmomatic -->
+      <options>ILLUMINACLIP:${AdapterFile}:2:30:12:1:true LEADING:3 TRAILING:3 MAXINFO:40:0.9 MINLEN:20</options>
+    </trimmomatic>
+
     <pe-filter>
       <!-- parallel environment option to the queue system -->

extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/dao/Annotationtype.java

@@ -867 +867 @@
 
   /**
+    The "PT_READS" annotation, used for derived bioassays (MergedSequences).
+    The annotation is the number of sequences that remains for a library after 
+    running trimmomatic.
+    @since 2.16
+  */
+  public static final Annotationtype PT_READS =
+    new Annotationtype("PT_READS", Type.LONG, false, Item.DERIVEDBIOASSAY);
+
+  /**
     The "PF_NNNN_PCT" annotation, used for derived bioassays (DemuxedSequences).
     The annotation is the percentage of sequences that has could not be assigned

extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/servlet/DemuxMergeServlet.java

@@ -13 +13 @@
 import java.util.Set;
 import java.util.TreeMap;
+import java.util.regex.Matcher;
 import java.util.regex.Pattern;
 
@@ -285 +286 @@
           mg.loadAnnotations(dc, "READS", Annotationtype.READS, null);
           mg.loadAnnotations(dc, "PF_READS", Annotationtype.PF_READS, null);
+          mg.loadAnnotations(dc, "PT_READS", Annotationtype.PT_READS, null);
           jsonMerged.add(mg.asJSONObject());
         }
@@ -376 +378 @@
         String debug_tileLimit = cluster.getConfig("demux/debug-tile-limit", "16");
         
+        String trim_jarPath = cluster.getRequiredConfig("trimmomatic/jar-path");
+        String trim_adapterFile = cluster.getRequiredConfig("trimmomatic/adapter-file");
+        String trim_options = cluster.getRequiredConfig("trimmomatic/options");
+        
         String projectRoot = cluster.getRequiredConfig("project-archive");
         String projectFolder = projectName+(debug ? ".debug" : "");
@@ -430 +436 @@
             script.cmd("export PicardDir="+demux_picardDir);
             script.cmd("export PicardMemory="+demux_picardMemory);
+            script.cmd("export TrimmomaticJAR="+trim_jarPath);
+            script.cmd("export AdapterFile="+trim_adapterFile);
             script.newLine();
             script.cmd("mkdir -p " + tmpFolder);
             script.cmd("cd " + tmpFolder);
             script.cmd("mkdir fastq");
+            script.cmd("mkdir fastq.merged");
+            script.cmd("mkdir fastq.trimmomatic");
             script.cmd("cp ${ScriptDir}/picard .");
+            script.cmd("cp ${ScriptDir}/trimmomatic .");
             script.newLine();
           }
@@ -555 +566 @@
           List<Tag> allTags = barcodeQuery.list(dc);
           String totalMetricsFile = jobFolder + "/demultiplex_metrics.txt";
+          String trimmomaticOut = jobFolder + "/trimmomatic.out";
 
           int currentLane = 0;
@@ -612 +624 @@
             
               script.comment(flowCellBarcode + "; lane " + lane);
-              // Progress between 5-75%
-              int percent = 5 + (currentLane * 70 / totalLanes);
+              // Progress between 5-50%
+              int percent = 5 + ((currentLane * 45) / totalLanes);
               currentLane++;
               script.progress(percent, "Extracting barcodes: " + flowCellBarcode + "; lane "+ lane);
@@ -632 +644 @@
               if (debug)
               {
-                //extractCmd += " TILE_LIMIT="+debug_tileLimit;
+                extractCmd += " TILE_LIMIT="+debug_tileLimit;
               }
               extractCmd += " >> " + demuxName + "_" + lane + ".out";
@@ -654 +666 @@
                 fastqCmd += " MACHINE_NAME="+sequencerName;
               fastqCmd += " NUM_PROCESSORS=0";
-              fastqCmd += " COMPRESS_OUTPUTS=true";
+              //fastqCmd += " COMPRESS_OUTPUTS=true";
               if (demux_fastqOptions != null)
               {
@@ -681 +693 @@
             String fastqFolder = (String)Annotationtype.DATA_FILES_FOLDER.getAnnotationValue(dc, merged);
             script.comment("Consolidating FASTQ: " + merged.getName());
-            // Progress between 75-95%
-            int percent = 75 + (currentToMove * 20 / totalToMove);
+            // Progress between 50-95%
+            int percent = 50 + ((currentToMove * 45) / totalToMove);
             currentToMove++;
+            
+            String R1_name = merged.getName()+"_R1.fastq.gz";
+            String R2_name = merged.getName()+"_R2.fastq.gz";
+            
             script.progress(percent, "Consolidating FASTQ: " + merged.getName());
+            script.bkgr("cat fastq/" + merged.getName() + "_*.1.fastq | gzip > fastq.merged/"+R1_name);
+            script.cmd("cat fastq/" + merged.getName() + "_*.2.fastq | gzip > fastq.merged/"+R2_name);
+            script.newLine();
+
+            script.progress(percent, "Trimmomatic: " + merged.getName());
+            script.cmd("echo [" + merged.getName() + "] >> " + trimmomaticOut);
+            String trimCmd = "./trimmomatic PE";
+            //trimCmd += " -trimlog fastq.trimmomatic/" + merged.getName()+".log";
+            trimCmd += " fastq.merged/"+R1_name;
+            trimCmd += " fastq.merged/"+R2_name;
+            trimCmd += " fastq.trimmomatic/"+R1_name;
+            trimCmd += " fastq.trimmomatic/un_"+R1_name;
+            trimCmd += " fastq.trimmomatic/"+R2_name;
+            trimCmd += " fastq.trimmomatic/un_"+R2_name;
+
+            if (trim_options != null)
+            {
+              trimCmd += " " + trim_options;
+            }
+            trimCmd += " >> "+trimmomaticOut+" 2>&1";
+            script.cmd(trimCmd);
+            script.newLine();
+            
+            script.progress(percent, "Archiving FASTQ: " + merged.getName());
             script.cmd("mkdir -p " + projectRoot + "/" + fastqFolder);
-            script.cmd("cat fastq/" + merged.getName() + "_*.1.fastq.gz > " + projectRoot + "/" + fastqFolder + "/"+merged.getName()+"_R1.fastq.gz");
-            script.cmd("cat fastq/" + merged.getName() + "_*.2.fastq.gz > " + projectRoot + "/" + fastqFolder + "/"+merged.getName()+"_R2.fastq.gz");
+            script.bkgr("cat fastq.trimmomatic/" + R1_name + " > " + projectRoot + "/" + fastqFolder + "/"+R1_name);
+            script.bkgr("cat fastq.trimmomatic/" + R2_name + " > " + projectRoot + "/" + fastqFolder + "/"+R2_name);
             script.newLine();
           }
@@ -977 +1017 @@
     {
       CmdResult metrics = cluster.executeCmd(ssh, "cat " + cluster.getJobFolder() + "/" + jobStatus.getJobName() + "/demultiplex_metrics.txt", 2);
-      if (metrics.getExitStatus() == 0)
-      {
-        Reads total = parseDemultiplexMetrics(sc, metrics.getStdout());
-        String msg = Values.formatNumber(total.reads/1000000f, 1) + "M reads; " + Values.formatNumber(total.passedFilter/1000000f, 1) + "M passed filter";
+      CmdResult trimmomatic = cluster.executeCmd(ssh, "cat " + cluster.getJobFolder() + "/" + jobStatus.getJobName() + "/trimmomatic.out", 2);
+      if (metrics.getExitStatus() != 0)
+      {
+        job.setDescription(metrics.getStderr());
+      }
+      else if (trimmomatic.getExitStatus() != 0)
+      {
+        job.setDescription(trimmomatic.getStderr());
+      }
+      else
+      {
+        Reads total = parseDemultiplexMetrics(sc, metrics.getStdout(), trimmomatic.getStdout());
+        String msg = Values.formatNumber(total.reads/1000000f, 1) + "M reads; ";
+        msg += Values.formatNumber(total.passedFilter/1000000f, 1) + "M passed filter; ";
+        msg += Values.formatNumber(total.passedTrimmomatic/1000000f, 1) + "M passed trimmomatic; ";
         if (total.warnings.size() > 0)
         {
@@ -987 +1038 @@
         return msg;
       }
-      else
-      {
-        job.setDescription(metrics.getStderr());
-      }
       return null;
     }
     
-    private Reads parseDemultiplexMetrics(SessionControl sc, String metrics)
-    {
+    private Reads parseDemultiplexMetrics(SessionControl sc, String metrics, String trimmomatic)
+    {
+      Map<String, Reads> sumReads = new HashMap<String, Reads>();
+
+      // Parse the demultiplex_metrics file
       FlatFileParser ffp = new FlatFileParser();
       ffp.setIgnoreRegexp(Pattern.compile("#.*"));
@@ -1005 +1055 @@
       ffp.setMinDataColumns(10);
 
-      Map<String, Reads> sumReads = new HashMap<String, Reads>();
       Set<String> demuxNames = new HashSet<String>();
       try
@@ -1094 +1143 @@
         throw new RuntimeException(ex);
       }
-            
+      
+      // Parse the trimmomatic.out file
+      Pattern libPattern = Pattern.compile("\\[(.*)\\]");
+      Pattern dataPattern = Pattern.compile(".*Both Surviving:\\s+(\\d+).*");
+      Reads currentLib = null;
+      int lineNo = 0;
+      for (String line : trimmomatic.split("\n"))
+      {
+        lineNo++;
+        Matcher m = libPattern.matcher(line);
+        if (m.matches())
+        {
+          String libName = m.group(1);
+          currentLib = sumReads.get(libName);
+          if (currentLib == null)
+          {
+            logger.error("At line " + lineNo + ": Found trimmomatic section for lib '" + libName + "' but not demultiplex metrics");
+          }
+          continue;
+        }
+        m = dataPattern.matcher(line);
+        if (m.matches())
+        {
+          if (currentLib == null)
+          {
+            logger.error("At line " + lineNo + ": Found trimmomatic data but has not found a library name");
+          }
+          else
+          {
+            currentLib.passedTrimmomatic = Values.getLong(m.group(1), -1);
+            if (logger.isDebugEnabled())
+            {
+              logger.debug("Trimmomatic: " + currentLib.libName + "; " + currentLib.passedTrimmomatic);
+            }
+            currentLib = null;
+          }
+        }
+      }
+      
       DbControl dc = null;
       Reads total = new Reads(null);
@@ -1144 +1231 @@
           if (logger.isDebugEnabled())
           {
-            logger.debug(r.libName + "; " + r.reads + "; " + r.passedFilter);
+            logger.debug(r.libName + "; " + r.reads + "; " + r.passedFilter + "; " + r.passedTrimmomatic);
           }
           
@@ -1150 +1237 @@
           if (merged != null)
           {
-            Annotationtype.READS.setAnnotationValue(dc, merged.getItem(), r.reads);
-            Annotationtype.PF_READS.setAnnotationValue(dc, merged.getItem(), r.passedFilter);
+            DerivedBioAssay m = merged.getItem();
+            Annotationtype.READS.setAnnotationValue(dc, m, r.reads);
+            Annotationtype.PF_READS.setAnnotationValue(dc, m, r.passedFilter);
+            Annotationtype.PT_READS.setAnnotationValue(dc, m, r.passedTrimmomatic);
             total.reads += r.reads;
             total.passedFilter += r.passedFilter;
+            total.passedTrimmomatic += r.passedTrimmomatic;
           }
         }
@@ -1174 +1264 @@
     long reads = 0;
     long passedFilter = 0;
+    long passedTrimmomatic = 0;
     
     Reads(String libName)

extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/servlet/InstallServlet.java

@@ -388 +388 @@
         jsonChecks.add(checkAnnotationType(dc, Annotationtype.READS, 1, null, effectivePermissionsUse, createIfMissing));
         jsonChecks.add(checkAnnotationType(dc, Annotationtype.PF_READS, 1, null, effectivePermissionsUse, createIfMissing));
+        jsonChecks.add(checkAnnotationType(dc, Annotationtype.PT_READS, 1, null, effectivePermissionsUse, createIfMissing));
         jsonChecks.add(checkAnnotationType(dc, Annotationtype.PF_NNNN_PCT, 1, null, effectivePermissionsUse, createIfMissing));
         jsonChecks.add(checkAnnotationType(dc, Annotationtype.PF_UNUSED_PCT, 1, null, effectivePermissionsUse, createIfMissing));
@@ -505 +506 @@
 
         jsonChecks.add(checkAnnotationTypeCategory(dc, Subtype.MERGED_SEQUENCES, createIfMissing, 
-            Annotationtype.ANALYSIS_RESULT, Annotationtype.READS, Annotationtype.PF_READS,
+            Annotationtype.ANALYSIS_RESULT, 
+            Annotationtype.READS, Annotationtype.PF_READS, Annotationtype.PT_READS,
             Annotationtype.DATA_FILES_FOLDER,
             Annotationtype.AUTO_PROCESSING