Changeset 5841

Timestamp:

Feb 25, 2020, 11:47:30 AM (3 years ago)

Author:

Nicklas Nordborg

Message:

References #1218: Implement MIPs alignment

Finalizing the MIPs alignment script by implementing the last two steps for collecting metrics. Results files are copied back to the project archive and linked in BASE to the alignment item.

Still need to implement importing some numbers as annotations from the metrics data.

Location:

extensions/net.sf.basedb.reggie/trunk

Files:

• config/reggie-config.xml (modified) (1 diff)
• src/net/sf/basedb/reggie/grid/MipsAlignJobCreator.java (modified) (10 diffs)

extensions/net.sf.basedb.reggie/trunk/config/reggie-config.xml

@@ -364 +364 @@
         <genome-fasta>${ReferenceDir}/gatk_bundle/2.8/b37/human_g1k_v37_decoy.fasta</genome-fasta>
         
-        <mark-duplicate-options>-CREATE_INDEX true -CREATE_MD5_FILE true -MAX_FILE_HANDLES 20000 -ALLOW_MISSING_UMIS false -DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES -MAX_EDIT_DISTANCE_TO_JOIN 1 -REMOVE_DUPLICATES true</mark-duplicate-options>
-        
+        <mark-duplicates-options>-CREATE_INDEX true -CREATE_MD5_FILE true -MAX_FILE_HANDLES 20000 -ALLOW_MISSING_UMIS false -DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES -MAX_EDIT_DISTANCE_TO_JOIN 1 -REMOVE_DUPLICATES true</mark-duplicates-options>
+        <pcr-metrics-options>-METRIC_ACCUMULATION_LEVEL null -METRIC_ACCUMULATION_LEVEL ALL_READS -MINIMUM_MAPPING_QUALITY 1 -MINIMUM_BASE_QUALITY 20 -CLIP_OVERLAPPING_READS true -COVERAGE_CAP 500 -NEAR_DISTANCE 5</pcr-metrics-options>
       </align-mips>
       

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

@@ -9 +9 @@
 import org.slf4j.LoggerFactory;
 
+import net.sf.basedb.core.AnyToAny;
+import net.sf.basedb.core.DataFileType;
 import net.sf.basedb.core.DbControl;
 import net.sf.basedb.core.DerivedBioAssay;
+import net.sf.basedb.core.Directory;
 import net.sf.basedb.core.File;
+import net.sf.basedb.core.FileServer;
+import net.sf.basedb.core.FileSetMember;
 import net.sf.basedb.core.InvalidDataException;
 import net.sf.basedb.core.ItemList;
@@ -17 +22 @@
 import net.sf.basedb.core.ItemSubtype;
 import net.sf.basedb.core.Job;
+import net.sf.basedb.core.Path;
 import net.sf.basedb.core.Protocol;
 import net.sf.basedb.core.Sample;
@@ -32 +38 @@
 import net.sf.basedb.reggie.Reggie;
 import net.sf.basedb.reggie.XmlConfig;
+import net.sf.basedb.reggie.dao.AlignedSequences;
 import net.sf.basedb.reggie.dao.Annotationtype;
 import net.sf.basedb.reggie.dao.BiomaterialList;
@@ -37 +44 @@
 import net.sf.basedb.reggie.dao.Dna;
 import net.sf.basedb.reggie.dao.DoNotUse;
+import net.sf.basedb.reggie.dao.Fileserver;
 import net.sf.basedb.reggie.dao.Library;
 import net.sf.basedb.reggie.dao.MergedSequences;
@@ -159 +167 @@
     String genome_dict = cfg.getRequiredConfig("align-mips/genome-dict", alignParameterSet);
     String genome_fasta = cfg.getRequiredConfig("align-mips/genome-fasta", alignParameterSet);
-    String mark_duplicate_options = cfg.getRequiredConfig("align-mips/mark-duplicate-options", alignParameterSet);
+    String mark_duplicates_options = cfg.getRequiredConfig("align-mips/mark-duplicates-options", alignParameterSet);
+    String pcr_metrics_options = cfg.getRequiredConfig("align-mips/pcr-metrics-options", alignParameterSet);
 
     // Load common items
@@ -294 +303 @@
       script.cmd("  " + Values.getString(fastqPrefix, "\n  ", true));
       script.cmd(")");
+      script.cmd("AlignFolder=${ArchiveRoot}"+alignFolder);
       script.cmd("TrimOptionsAdapter=\""+align_trimmomaticOptions1+"\"");
       script.cmd("TrimOptionsQual=\""+align_trimmomaticOptions2+"\"");
       script.cmd("NovoAlignOptions=\""+novoalign_options+"\"");
-      script.cmd("MarkDuplicatesOptions=\""+mark_duplicate_options+"\"");
+      script.cmd("MarkDuplicatesOptions=\""+mark_duplicates_options+"\"");
+      script.cmd("PcrMetricsOptions=\""+pcr_metrics_options+"\"");
       
       script.newLine();
@@ -364 +375 @@
       script.comment("Run novoalign");
       script.progress(40, "Aligning with novoalign");
-      script.cmd("min_insert=$(awk 'NR == 1 || $3 - $2 < min {min = $3 - $2}END{print min - 1}' \"${AmpliconsBed}\")");
-      script.cmd("max_insert=$(awk 'NR == 1 || $3 - $2 > max {max = $3 - $2}END{print max + 1}' \"${AmpliconsBed}\")");
+      script.cmd("get_min_max_insert");
       script.cmd("n=0");
       script.cmd("AlignedBams=''");
@@ -377 +387 @@
       script.newLine();
 
-      script.comment("Re-header BAM file");
+      script.comment("Re-header and split BAM file");
       script.progress(70, "Post-processing aligned BAM file");
       script.cmd("reheader_bam");
-      script.newLine();
-      
-      script.comment("Split BAM file into concordant, discordant and unmapped");
       script.cmd("split_bam");
       script.newLine();
@@ -392 +399 @@
       script.newLine();
 
-      /*
+      script.comment("Collect metrics");
+      script.progress(90, "Collecting PCR and alignment metrics");
+      script.cmd("collect_pcr_metrics");
+      script.cmd("collect_alignment_metrics");
+      script.newLine();
+      
       script.progress(95, "Copying result files to project archive");
-      script.cmd("mkdir -p ${HisatFolder}");
-      script.cmd("rm -rf ${HisatFolder}/*");
-      script.cmd("cp masked/*.out ${WD}");
-      script.cmd("cp aligned/*.out ${WD}");
-      script.cmd("cp aligned/alignment_picardmetrics.csv ${WD}");
-      script.cmd("cp aligned/* ${HisatFolder}");
+      script.cmd("mkdir -p ${AlignFolder}");
+      script.cmd("rm -rf ${AlignFolder}/*");
+      script.cmd("cp out/* ${AlignFolder}");
+      // TODO - copy some files to $WD so that we can update some annotations like ALIGNED_PAIRS, etc.
       if (externalGroup != null)
       {
-        script.cmd("chgrp -R " + externalGroup + " ${HisatFolder} 2>> ${WD}/chgrp.out || echo [" + alignedName +"] >> ${WD}/chgrp.out");
-      }
-      script.cmd("ls -1 ${HisatFolder}/* >> ${WD}/files.out");
-      script.newLine();
-      */
+        script.cmd("chgrp -R " + externalGroup + " ${AlignFolder} 2>> ${WD}/chgrp.out || echo [" + alignedName +"] >> ${WD}/chgrp.out");
+      }
+      script.cmd("ls -1 ${AlignFolder}/* >> ${WD}/files.out");
+      script.newLine();
       
       script.progress(99, "Cleaning up temporary folders");
@@ -435 +444 @@
     public String jobCompleted(SessionControl sc, OpenGridSession session, Job job, JobStatus status)
     {
+      String jobName = status.getName();
+      String files = session.getJobFileAsString(jobName, "files.out", "UTF-8");
+
+      Metrics metrics = parseAlignedOut(sc, job, files);
       return null;
     }
     
+    private Metrics parseAlignedOut(SessionControl sc, Job job, String filesOut)
+    {
+      Metrics metrics = new Metrics();
+      
+      DbControl dc = null;
+      try
+      {
+        dc = sc.newDbControl();
+        
+        AlignedSequences alignedSequences = AlignedSequences.getByJob(dc, job);
+        DerivedBioAssay aligned = alignedSequences.getItem();
+
+        // Create file links
+        boolean useExternalProjectArchive = Reggie.isExternalItem(aligned.getName());
+        FileServer fileArchive = useExternalProjectArchive ? Fileserver.EXTERNAL_ARCHIVE.load(dc) : Fileserver.PROJECT_ARCHIVE.load(dc);
+        String analysisDir = useExternalProjectArchive ? Reggie.EXTERNAL_ANALYSIS_DIR : Reggie.SECONDARY_ANALYSIS_DIR;
+
+        String dataFilesFolder = (String)Annotationtype.DATA_FILES_FOLDER.getAnnotationValue(dc, aligned);
+        String baseFolder = Reggie.convertDataFilesFolderToBaseFolder(dataFilesFolder);
+        Directory localDataDir = Directory.getNew(dc, new Path(analysisDir+baseFolder, Path.Type.DIRECTORY));
+        DataFileType bamData = Datafiletype.BAM.load(dc);
+        ItemSubtype bamType = bamData.getGenericType();
+  
+        int lineNo = 0;
+        for (String line : filesOut.split("\n"))
+        {
+          lineNo++;
+          
+          File f = File.getFile(dc, localDataDir, line.substring(line.lastIndexOf("/")+1), true);
+          f.setFileServer(fileArchive);
+          String fileUrl = "sftp://" + fileArchive.getHost() + dataFilesFolder + "/" + f.getName();
+          try
+          {
+            f.setUrl(fileUrl, true);
+          }
+          catch (RuntimeException ex)
+          {
+            f.setUrl(fileUrl, false);
+          }
+          
+          if (!f.isInDatabase())
+          {
+            dc.saveItem(f);
+          }
+          if (f.getName().equals("concordant.bam"))
+          {
+            //f.setDescription(metrics.numReadsAfterAlign + " ALIGNED PAIRS");
+            f.setItemSubtype(bamType);
+            FileSetMember member = aligned.getFileSet().addMember(f, bamData);
+          }
+          else
+          {
+            AnyToAny link = AnyToAny.getNewOrExisting(dc, aligned, f.getName(), f, true);
+            if (!link.isInDatabase()) dc.saveItem(link);
+          }
+        }
+        
+        dc.commit();
+      }
+      finally
+      {
+        if (dc != null) dc.close();
+      }
+  
+      return metrics;
+    }
+
+  }
+
+  static class Metrics
+  {
   }