Changeset 5833

Timestamp:

Feb 24, 2020, 11:05:56 AM (3 years ago)

Author:

Nicklas Nordborg

Message:

References #1218: Implement MIPs alignment

Started to re-design the MIPs alignment script. Moved a lot of functionality into an external script (mips_functions.sh). The Trimmomatic script has been fixed including logging and error handling.

File:

• extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/MipsAlignJobCreator.java (modified) (3 diffs)

extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/MipsAlignJobCreator.java

@@ -291 +291 @@
       script.cmd("ArchiveRoot="+archiveRoot);
       script.cmd("FastqFolder=${ArchiveRoot}"+fastQFolder);
-      script.cmd("FastqPrefix=( \\");
-      script.cmd("  " + Values.getString(fastqPrefix, " \\\n  ", true) + " \\");
-      script.cmd("  )");
+      script.cmd("FastqPrefix=( ");
+      script.cmd("  " + Values.getString(fastqPrefix, "\n  ", true));
+      script.cmd(")");
+      script.cmd("TrimOptionsAdapter=\""+align_trimmomaticOptions1+"\"");
+      script.cmd("TrimOptionsQual=\""+align_trimmomaticOptions2+"\"");
       
       script.newLine();
@@ -309 +311 @@
       script.cmd("cd ${TMPDIR}");
       script.cmd("mkdir fastq");
-      script.cmd("mkdir trim.adapter");
-      script.cmd("mkdir trim.quality");
-      script.cmd("mkdir bam");
-      script.cmd("mkdir bam.umi");
-      script.cmd("mkdir bam.aligned");
-      script.cmd("mkdir bam.sorted");
-      script.cmd("mkdir bam.merged");
-      script.cmd("mkdir bam.final");
+      script.cmd("mkdir tmp");
+      script.cmd("mkdir out");
       script.cmd("cp ${ScriptDir}/stdwrap.sh .");
       script.cmd("cp ${ScriptDir}/stderrwrap.sh .");
-      script.cmd("cp ${ScriptDir}/mips_adapters.sh .");
       script.cmd("cp ${ScriptDir}/picard2 .");
-      script.cmd("cp ${ScriptDir}/convert_picard_parameters.sh .");
-      script.newLine();
-
-      script.comment("Setup adapters in a separate script");
-      script.cmd(". ./mips_adapters.sh");
+      script.cmd("cp ${ScriptDir}/mips_functions.sh .");
+      script.newLine();
+
+      script.comment("Load script with helper functions");
+      script.cmd(". ./mips_functions.sh");
       script.newLine();
       
       script.comment("Copy FASTQ and .rg files to tmp folder");
       script.progress(10, "Copying FASTQ files");
-      script.cmd("if [ ! -d \"${FastqFolder}\" ]; then");
-      script.cmd("echo \"Can't find FASTQ data folder ${FastqFolder}\" 1>&2");
+      script.cmd("copy_fastq_and_rg \"${FastqFolder}\" \"fastq\"");
+      script.newLine();
+      
+      script.comment("Run Trimmomatic");
+      script.progress(20, "Running Trimmomatic");
+      script.comment("Create adapter.fa for use with trimmomatic");
+      script.cmd("create_adapter_fa \"adapter.fa\"");
+      script.newLine();
+      
+      script.cmd("n=0");
+      script.cmd("nAfterTrim=0");
+      script.cmd("for prefix in ${FastqPrefix[@]} ; do");
+      script.cmd("   n=$((n + 1))");
+      script.cmd("   run_trimmomatic \"${prefix}\" $n");
+      script.cmd("done");
+      script.cmd("trimlog_stats \"tmp/*.adaptrim.log\" \"out/adaptrim.hist.txt\"");
+      script.cmd("trimlog_stats \"tmp/*.qualtrim.log\" \"out/qualtrim.hist.txt\"");
+      script.cmd("if [ ${nAfterTrim} -eq 0 ]; then");
+      script.cmd("  echo \"No paired end reads surviving adapter and quality trimming for any fastq prefix\" 1>&2");
       script.cmd("exit 1");
       script.cmd("fi");
-      script.cmd("cp ${FastqFolder}/*.fastq.gz fastq");
-      script.cmd("cp ${FastqFolder}/*.rg fastq");
-      script.newLine();
-
-      script.comment("Run Trimmomatic");
-      script.cmd("for prefix in ${FastqPrefix[@]} ; do");
-      
-      String trimCmd1 = "./stdwrap.sh ${JAVA} -jar ${TrimmomaticJAR} PE";
-      trimCmd1 += " -threads ${NumThreads}";
-      trimCmd1 += " -phred33 -trimlog trim.adapter/${prefix}.log";
-      trimCmd1 += " fastq/${prefix}_R1.fastq.gz";
-      trimCmd1 += " fastq/${prefix}_R2.fastq.gz";
-      trimCmd1 += " trim.adapter/${prefix}_R1.fastq";
-      trimCmd1 += " trim.adapter/un_${prefix}_R1.fastq";
-      trimCmd1 += " trim.adapter/${prefix}_R2.fastq";
-      trimCmd1 += " trim.adapter/un_${prefix}_R2.fastq";
-      if (align_trimmomaticOptions1 != null)
-      {
-        trimCmd1 += " " + align_trimmomaticOptions1;
-      }
-      trimCmd1 += " >> trim.adapter/trimmomatic.out";
-      
-      String trimCmd2 = "./stdwrap.sh ${JAVA} -jar ${TrimmomaticJAR} PE";
-      trimCmd2 += " -threads ${NumThreads}";
-      trimCmd2 += " -phred33 -trimlog trim.quality/${prefix}.log";
-      trimCmd2 += " trim.adapter/${prefix}_R1.fastq";
-      trimCmd2 += " trim.adapter/${prefix}_R2.fastq";
-      trimCmd2 += " trim.quality/${prefix}_R1.fastq";
-      trimCmd2 += " trim.quality/un_${prefix}_R1.fastq";
-      trimCmd2 += " trim.quality/${prefix}_R2.fastq";
-      trimCmd2 += " trim.quality/un_${prefix}_R2.fastq";
-      if (align_trimmomaticOptions2 != null)
-      {
-        trimCmd2 += " " + align_trimmomaticOptions2;
-      }
-      trimCmd2 += " >> trim.quality/trimmomatic.out";
-      script.cmd("   " + trimCmd1);
-      script.cmd("   " + trimCmd2);
-      script.cmd("done");
-      script.newLine();
-      
+      
+      /*
       script.comment("Run Picard FastqToSam and annotat the BAM files with UMI");
       script.cmd("for prefix in ${FastqPrefix[@]} ; do");
@@ -467 +440 @@
       markDupCmd += " >> markduplicates.out";
       script.cmd(markDupCmd);
-      
+      */
       /*
       script.progress(95, "Copying result files to project archive");