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

Code
Comments
Other
Rev Date Author Line
4590 25 Sep 17 nicklas 1 package net.sf.basedb.reggie.grid;
4590 25 Sep 17 nicklas 2
4622 16 Nov 17 nicklas 3 import java.io.InputStream;
4622 16 Nov 17 nicklas 4 import java.io.OutputStream;
4590 25 Sep 17 nicklas 5 import java.util.ArrayList;
4590 25 Sep 17 nicklas 6 import java.util.Arrays;
4590 25 Sep 17 nicklas 7 import java.util.List;
4590 25 Sep 17 nicklas 8 import java.util.regex.Matcher;
4590 25 Sep 17 nicklas 9 import java.util.regex.Pattern;
4590 25 Sep 17 nicklas 10
4590 25 Sep 17 nicklas 11 import org.slf4j.LoggerFactory;
4590 25 Sep 17 nicklas 12
4590 25 Sep 17 nicklas 13 import net.sf.basedb.core.AnyToAny;
5579 20 Aug 19 nicklas 14 import net.sf.basedb.core.BaseException;
4590 25 Sep 17 nicklas 15 import net.sf.basedb.core.DataFileType;
4590 25 Sep 17 nicklas 16 import net.sf.basedb.core.DbControl;
4590 25 Sep 17 nicklas 17 import net.sf.basedb.core.DerivedBioAssay;
4590 25 Sep 17 nicklas 18 import net.sf.basedb.core.Directory;
4590 25 Sep 17 nicklas 19 import net.sf.basedb.core.File;
4590 25 Sep 17 nicklas 20 import net.sf.basedb.core.FileServer;
4590 25 Sep 17 nicklas 21 import net.sf.basedb.core.FileSetMember;
4598 27 Sep 17 nicklas 22 import net.sf.basedb.core.ItemList;
4590 25 Sep 17 nicklas 23 import net.sf.basedb.core.ItemNotFoundException;
4590 25 Sep 17 nicklas 24 import net.sf.basedb.core.ItemSubtype;
4590 25 Sep 17 nicklas 25 import net.sf.basedb.core.Job;
4590 25 Sep 17 nicklas 26 import net.sf.basedb.core.Path;
4590 25 Sep 17 nicklas 27 import net.sf.basedb.core.Protocol;
4590 25 Sep 17 nicklas 28 import net.sf.basedb.core.Sample;
4590 25 Sep 17 nicklas 29 import net.sf.basedb.core.SessionControl;
4590 25 Sep 17 nicklas 30 import net.sf.basedb.core.Software;
5547 07 Aug 19 nicklas 31 import net.sf.basedb.core.StringParameterType;
4590 25 Sep 17 nicklas 32 import net.sf.basedb.opengrid.JobDefinition;
4590 25 Sep 17 nicklas 33 import net.sf.basedb.opengrid.JobStatus;
4590 25 Sep 17 nicklas 34 import net.sf.basedb.opengrid.OpenGridCluster;
4590 25 Sep 17 nicklas 35 import net.sf.basedb.opengrid.OpenGridSession;
4590 25 Sep 17 nicklas 36 import net.sf.basedb.opengrid.ScriptBuilder;
4590 25 Sep 17 nicklas 37 import net.sf.basedb.opengrid.config.ClusterConfig;
4590 25 Sep 17 nicklas 38 import net.sf.basedb.opengrid.config.JobConfig;
6640 11 Mar 22 nicklas 39 import net.sf.basedb.opengrid.filetransfer.FilePermission;
4590 25 Sep 17 nicklas 40 import net.sf.basedb.opengrid.service.JobCompletionHandler;
4590 25 Sep 17 nicklas 41 import net.sf.basedb.reggie.Reggie;
4590 25 Sep 17 nicklas 42 import net.sf.basedb.reggie.XmlConfig;
4590 25 Sep 17 nicklas 43 import net.sf.basedb.reggie.dao.AlignedSequences;
4590 25 Sep 17 nicklas 44 import net.sf.basedb.reggie.dao.Annotationtype;
4598 27 Sep 17 nicklas 45 import net.sf.basedb.reggie.dao.BiomaterialList;
4590 25 Sep 17 nicklas 46 import net.sf.basedb.reggie.dao.Datafiletype;
5791 16 Dec 19 nicklas 47 import net.sf.basedb.reggie.dao.DoNotUse;
4590 25 Sep 17 nicklas 48 import net.sf.basedb.reggie.dao.Fileserver;
4592 26 Sep 17 nicklas 49 import net.sf.basedb.reggie.dao.FlowCell;
4590 25 Sep 17 nicklas 50 import net.sf.basedb.reggie.dao.Library;
4656 24 Jan 18 nicklas 51 import net.sf.basedb.reggie.dao.MaskedSequences;
4590 25 Sep 17 nicklas 52 import net.sf.basedb.reggie.dao.MergedSequences;
5547 07 Aug 19 nicklas 53 import net.sf.basedb.reggie.dao.Pipeline;
4590 25 Sep 17 nicklas 54 import net.sf.basedb.reggie.dao.Subtype;
4622 16 Nov 17 nicklas 55 import net.sf.basedb.reggie.vcf.VcfData;
4622 16 Nov 17 nicklas 56 import net.sf.basedb.reggie.vcf.VcfParser;
4622 16 Nov 17 nicklas 57 import net.sf.basedb.util.FileUtil;
4622 16 Nov 17 nicklas 58 import net.sf.basedb.util.InputStreamSplitter;
4590 25 Sep 17 nicklas 59 import net.sf.basedb.util.Values;
7079 27 Mar 23 nicklas 60 import net.sf.basedb.util.extensions.logging.ExtensionsLog;
7079 27 Mar 23 nicklas 61 import net.sf.basedb.util.extensions.logging.ExtensionsLogger;
4590 25 Sep 17 nicklas 62
4590 25 Sep 17 nicklas 63 /**
4590 25 Sep 17 nicklas 64   Helper class for creating items needed for executing Hisat as
4590 25 Sep 17 nicklas 65   well as generating the script and send it to the cluster for
4590 25 Sep 17 nicklas 66   execution.
4590 25 Sep 17 nicklas 67   
4590 25 Sep 17 nicklas 68   @author nicklas
4590 25 Sep 17 nicklas 69   @since 4.12
4590 25 Sep 17 nicklas 70 */
4590 25 Sep 17 nicklas 71 public class HisatAlignJobCreator 
6674 11 Apr 22 nicklas 72   extends AbstractJobCreator
4590 25 Sep 17 nicklas 73 {
4656 24 Jan 18 nicklas 74   private Software maskingSoftware;
4656 24 Jan 18 nicklas 75   private Protocol maskingProtocol;
4656 24 Jan 18 nicklas 76   
4590 25 Sep 17 nicklas 77   private Software alignSoftware;
4590 25 Sep 17 nicklas 78   private Protocol alignProtocol;
4590 25 Sep 17 nicklas 79
4590 25 Sep 17 nicklas 80   public HisatAlignJobCreator()
4590 25 Sep 17 nicklas 81   {}
4590 25 Sep 17 nicklas 82
4656 24 Jan 18 nicklas 83   /**
4656 24 Jan 18 nicklas 84     Set the software item to set on created MaskedSequences.
4656 24 Jan 18 nicklas 85     @see DerivedBioAssay#setSoftware(Software)
4656 24 Jan 18 nicklas 86     @since 4.15
4656 24 Jan 18 nicklas 87   */
4656 24 Jan 18 nicklas 88   public void setMaskingSoftware(Software software)
4656 24 Jan 18 nicklas 89   {
4656 24 Jan 18 nicklas 90     this.maskingSoftware = software;
4656 24 Jan 18 nicklas 91   }
4590 25 Sep 17 nicklas 92   
4590 25 Sep 17 nicklas 93   /**
4656 24 Jan 18 nicklas 94     Set the protocol item to set on created MaskedSequences.
4656 24 Jan 18 nicklas 95     @see DerivedBioAssay#setProtocol(Protocol)
4656 24 Jan 18 nicklas 96     @since 4.15
4656 24 Jan 18 nicklas 97   */
4656 24 Jan 18 nicklas 98   public void setMaskingProtocol(Protocol protocol)
4656 24 Jan 18 nicklas 99   {
4656 24 Jan 18 nicklas 100     this.maskingProtocol = protocol;
4656 24 Jan 18 nicklas 101   }
4656 24 Jan 18 nicklas 102
4656 24 Jan 18 nicklas 103   /**
4590 25 Sep 17 nicklas 104     Set the software item to set on created AlignedSequences.
4590 25 Sep 17 nicklas 105     @see DerivedBioAssay#setSoftware(Software)
4590 25 Sep 17 nicklas 106   */
4590 25 Sep 17 nicklas 107   public void setAlignSoftware(Software software)
4590 25 Sep 17 nicklas 108   {
4590 25 Sep 17 nicklas 109     this.alignSoftware = software;
4590 25 Sep 17 nicklas 110   }
4590 25 Sep 17 nicklas 111   
4590 25 Sep 17 nicklas 112   /**
4590 25 Sep 17 nicklas 113     Set the protocol item to set on created AlignedSequences.
4590 25 Sep 17 nicklas 114     @see DerivedBioAssay#setProtocol(Protocol)
4590 25 Sep 17 nicklas 115   */
4590 25 Sep 17 nicklas 116   public void setAlignProtocol(Protocol protocol)
4590 25 Sep 17 nicklas 117   {
4590 25 Sep 17 nicklas 118     this.alignProtocol = protocol;
4590 25 Sep 17 nicklas 119   }
4590 25 Sep 17 nicklas 120     
4590 25 Sep 17 nicklas 121   /**
4590 25 Sep 17 nicklas 122     Create a child bioassays for all given merged sequences and schedule
4590 25 Sep 17 nicklas 123     jobs on the given cluster for running Hisat alignment.
4590 25 Sep 17 nicklas 124     @return A list with the corresponding jobs in BASE
4590 25 Sep 17 nicklas 125   */
4599 28 Sep 17 nicklas 126   public List<JobDefinition> createHisatAlignJobs(DbControl dc, OpenGridCluster cluster, List<MergedSequences> mergedSequences)
4590 25 Sep 17 nicklas 127   {
4590 25 Sep 17 nicklas 128     SessionControl sc = dc.getSessionControl();
4590 25 Sep 17 nicklas 129
4590 25 Sep 17 nicklas 130     ClusterConfig clusterCfg = cluster.getConfig();
4590 25 Sep 17 nicklas 131     XmlConfig cfg = Reggie.getConfig(cluster.getId());
4590 25 Sep 17 nicklas 132     if (cfg == null)
4590 25 Sep 17 nicklas 133     {
4590 25 Sep 17 nicklas 134       throw new ItemNotFoundException("No configuration in reggie-config.xml for cluster: " + cluster.getId());
4590 25 Sep 17 nicklas 135     }
4590 25 Sep 17 nicklas 136     
4656 24 Jan 18 nicklas 137     String maskParameterSet = (String)Annotationtype.PARAMETER_SET.getAnnotationValue(dc, maskingSoftware);
4590 25 Sep 17 nicklas 138     String alignParameterSet = (String)Annotationtype.PARAMETER_SET.getAnnotationValue(dc, alignSoftware);
4590 25 Sep 17 nicklas 139     
4590 25 Sep 17 nicklas 140     // Get global options
6693 22 Apr 22 nicklas 141     String global_env = ScriptUtil.multilineIndent(cfg.getConfig("global-env"));
6654 23 Mar 22 nicklas 142     String projectArchive = cfg.getRequiredConfig("project-archive", null);
6654 23 Mar 22 nicklas 143     String externalArchive = cfg.getConfig("external-archive", null, projectArchive);
4590 25 Sep 17 nicklas 144
4590 25 Sep 17 nicklas 145     // Options for the programs
7372 06 Oct 23 nicklas 146     String align_submit = cfg.getConfig("align-hisat/submit", alignParameterSet, null);
7372 06 Oct 23 nicklas 147     String align_submit_debug = cfg.getConfig("align-hisat/submit-debug", alignParameterSet, null);
6654 23 Mar 22 nicklas 148     String hisat_env = ScriptUtil.multilineIndent(cfg.getRequiredConfig("align-hisat/env", alignParameterSet));
6665 05 Apr 22 nicklas 149     String hisat_envdebug = ScriptUtil.multilineIndent(cfg.getConfig("align-hisat/env-debug", alignParameterSet, null));
6654 23 Mar 22 nicklas 150     String hisat_execute = ScriptUtil.multilineIndent(cfg.getConfig("align-hisat/execute", alignParameterSet, "./hisat.sh"));
4590 25 Sep 17 nicklas 151
4590 25 Sep 17 nicklas 152     // Load common items
4656 24 Jan 18 nicklas 153     ItemSubtype maskedType = Subtype.MASKED_SEQUENCES.get(dc);
4590 25 Sep 17 nicklas 154     ItemSubtype alignType = Subtype.ALIGNED_SEQUENCES.get(dc);
4590 25 Sep 17 nicklas 155     ItemSubtype hisatAlignJobType = Subtype.HISAT_ALIGN_JOB.get(dc);
4598 27 Sep 17 nicklas 156   
4598 27 Sep 17 nicklas 157     // Selected items must be removed from this list
4598 27 Sep 17 nicklas 158     ItemList hisatPipeline = BiomaterialList.HISAT_PIPELINE.load(dc);
4590 25 Sep 17 nicklas 159
4590 25 Sep 17 nicklas 160     // Options common for all jobs
4590 25 Sep 17 nicklas 161     JobConfig jobConfig = new JobConfig();
4590 25 Sep 17 nicklas 162     if (priority != null) jobConfig.setPriority(priority);
7372 06 Oct 23 nicklas 163     if (partition != null) jobConfig.setSbatchOption("partition", ScriptUtil.checkValidScriptParameter(partition));
7372 06 Oct 23 nicklas 164     jobConfig.convertOptionsTo(clusterCfg.getType());
7372 06 Oct 23 nicklas 165     if (submitOptionsOverride != null)
7372 06 Oct 23 nicklas 166     {
7372 06 Oct 23 nicklas 167       ScriptUtil.addSubmitOptions(jobConfig, submitOptionsOverride, clusterCfg.getType());
7372 06 Oct 23 nicklas 168     }
7372 06 Oct 23 nicklas 169     else
7372 06 Oct 23 nicklas 170     {
7372 06 Oct 23 nicklas 171       ScriptUtil.addSubmitOptions(jobConfig, align_submit, clusterCfg.getType());
7372 06 Oct 23 nicklas 172       if (debug) ScriptUtil.addSubmitOptions(jobConfig, align_submit_debug, clusterCfg.getType());
7372 06 Oct 23 nicklas 173     }
4590 25 Sep 17 nicklas 174     
4590 25 Sep 17 nicklas 175     // We submit one job for each raw bioassay to the cluster
4590 25 Sep 17 nicklas 176     List<JobDefinition> jobDefs = new ArrayList<JobDefinition>(mergedSequences.size());
4590 25 Sep 17 nicklas 177
4599 28 Sep 17 nicklas 178     for (MergedSequences ms : mergedSequences)
4590 25 Sep 17 nicklas 179     {
4623 17 Nov 17 nicklas 180       ms = MergedSequences.getById(dc, ms.getId()); // Ensure item is loaded in this transaction
4599 28 Sep 17 nicklas 181       // Get some information about the merged data that we need
4599 28 Sep 17 nicklas 182       DerivedBioAssay merged = ms.getDerivedBioAssay();
5364 16 Apr 19 nicklas 183       hisatPipeline.removeItem(merged);
4590 25 Sep 17 nicklas 184       
4599 28 Sep 17 nicklas 185       Library lib = Library.get(merged.getExtract());
5596 11 Sep 19 nicklas 186       boolean isExternal = Reggie.isExternalItem(merged.getName());
6654 23 Mar 22 nicklas 187       String archiveFolder = isExternal ? externalArchive : projectArchive;
4599 28 Sep 17 nicklas 188       
4599 28 Sep 17 nicklas 189       Sample specimen =  (Sample)lib.findSingleParent(dc, Subtype.SPECIMEN);
4599 28 Sep 17 nicklas 190       String fastQFolder = (String)Annotationtype.DATA_FILES_FOLDER.getAnnotationValue(dc, merged);
4599 28 Sep 17 nicklas 191       
4599 28 Sep 17 nicklas 192       List<FlowCell> flowCells = ms.getFlowCells(dc);
4599 28 Sep 17 nicklas 193       List<String> flowCellIds = new ArrayList<>();
4599 28 Sep 17 nicklas 194       for (FlowCell fc : flowCells)
4590 25 Sep 17 nicklas 195       {
4599 28 Sep 17 nicklas 196         String fcId = (String)Annotationtype.FLOWCELL_ID.getAnnotationValue(dc, fc.getItem());
4599 28 Sep 17 nicklas 197         if (fcId != null) flowCellIds.add(ScriptUtil.checkValidScriptParameter(fcId));
4599 28 Sep 17 nicklas 198       }
4599 28 Sep 17 nicklas 199       
4599 28 Sep 17 nicklas 200       // Create job 
4599 28 Sep 17 nicklas 201       Job alignJob = Job.getNew(dc, null, null, null);
4599 28 Sep 17 nicklas 202       alignJob.setItemSubtype(hisatAlignJobType);
4599 28 Sep 17 nicklas 203       alignJob.setPluginVersion("reggie-"+Reggie.VERSION);
4599 28 Sep 17 nicklas 204       alignJob.setSendMessage(Values.getBoolean(sc.getUserClientSetting("plugins.sendmessage"), false));
4599 28 Sep 17 nicklas 205       alignJob.setName("HisatAlign " + merged.getName());
5547 07 Aug 19 nicklas 206       alignJob.setParameterValue("pipeline", new StringParameterType(), Pipeline.RNASEQ_HISAT_STRINGTIE.getId());
4599 28 Sep 17 nicklas 207       if (debug) alignJob.setName(alignJob.getName() + " (debug)");
6981 17 Jan 23 nicklas 208       if (partition != null) alignJob.setParameterValue("partition", new StringParameterType(), partition);
7372 06 Oct 23 nicklas 209       if (submitOptionsOverride != null) alignJob.setParameterValue("jobOptions", new StringParameterType(), submitOptionsOverride);
4599 28 Sep 17 nicklas 210       dc.saveItem(alignJob);
4590 25 Sep 17 nicklas 211
4656 24 Jan 18 nicklas 212       // Created MASKED derived bioassay set
4656 24 Jan 18 nicklas 213       String maskedName = ms.getNextMaskedSequencesName(dc);
4656 24 Jan 18 nicklas 214       DerivedBioAssay masked = DerivedBioAssay.getNew(dc, false, alignJob);
4656 24 Jan 18 nicklas 215       masked.setItemSubtype(maskedType);
5547 07 Aug 19 nicklas 216       Pipeline.RNASEQ_HISAT_STRINGTIE.setAnnotation(dc, masked);
4881 03 Jul 18 nicklas 217       masked.setName(maskedName);
4656 24 Jan 18 nicklas 218       masked.setExtract(merged.getExtract());
4656 24 Jan 18 nicklas 219       masked.setSoftware(maskingSoftware);
4656 24 Jan 18 nicklas 220       masked.setProtocol(maskingProtocol);
4656 24 Jan 18 nicklas 221       masked.addParent(merged);
5791 16 Dec 19 nicklas 222       DoNotUse.copyDoNotUseAnnotations(dc, merged, masked, false);
4656 24 Jan 18 nicklas 223       dc.saveItem(masked);
4656 24 Jan 18 nicklas 224       
4599 28 Sep 17 nicklas 225       // Create ALIGNED derived bioassay set
4881 03 Jul 18 nicklas 226       String alignedName = masked.getName()+".a"; // This is safe since the masked item is a new item
4599 28 Sep 17 nicklas 227       DerivedBioAssay aligned = DerivedBioAssay.getNew(dc, false, alignJob);
4599 28 Sep 17 nicklas 228       aligned.setItemSubtype(alignType);
5547 07 Aug 19 nicklas 229       Pipeline.RNASEQ_HISAT_STRINGTIE.setAnnotation(dc, aligned);
4881 03 Jul 18 nicklas 230       aligned.setName(alignedName);
4656 24 Jan 18 nicklas 231       aligned.setExtract(masked.getExtract());
4599 28 Sep 17 nicklas 232       aligned.setSoftware(alignSoftware);
4599 28 Sep 17 nicklas 233       aligned.setProtocol(alignProtocol);
4656 24 Jan 18 nicklas 234       aligned.addParent(masked);
5791 16 Dec 19 nicklas 235       DoNotUse.copyDoNotUseAnnotations(dc, merged, aligned, false);
4599 28 Sep 17 nicklas 236       dc.saveItem(aligned);
4599 28 Sep 17 nicklas 237       
4599 28 Sep 17 nicklas 238       String hisatFolder = fastQFolder + "/"+aligned.getName().substring(merged.getName().length()+1);
4599 28 Sep 17 nicklas 239       if (debug && !hisatFolder.startsWith("/debug"))
4599 28 Sep 17 nicklas 240       {
4599 28 Sep 17 nicklas 241         hisatFolder = "/debug" + hisatFolder;
4599 28 Sep 17 nicklas 242       }
4599 28 Sep 17 nicklas 243       Annotationtype.DATA_FILES_FOLDER.setAnnotationValue(dc, aligned, hisatFolder);
4599 28 Sep 17 nicklas 244       if (autoConfirm)
4599 28 Sep 17 nicklas 245       {
4599 28 Sep 17 nicklas 246         Annotationtype.AUTO_PROCESSING.setAnnotationValue(dc, aligned, "AutoConfirm");
4599 28 Sep 17 nicklas 247       }
4590 25 Sep 17 nicklas 248
4599 28 Sep 17 nicklas 249       // Checks to make sure no bad things are included in script file
4599 28 Sep 17 nicklas 250       ScriptUtil.checkValidPath(hisatFolder, true, true);
4599 28 Sep 17 nicklas 251       String mergedName = ScriptUtil.checkValidScriptParameter(merged.getName());
4599 28 Sep 17 nicklas 252       String libName = ScriptUtil.checkValidScriptParameter(lib.getName());
4599 28 Sep 17 nicklas 253       // Remove all after the first dot not followed by a number
4599 28 Sep 17 nicklas 254       String sampleName = libName.replaceFirst("\\.(?!\\d).*", "");
4599 28 Sep 17 nicklas 255       if (specimen != null && specimen.getExternalId() != null)
4599 28 Sep 17 nicklas 256       {
4599 28 Sep 17 nicklas 257         // Replace SCANB-ID in libName/mergedName/sampleName with Sample.externalId
4599 28 Sep 17 nicklas 258         sampleName = ScriptUtil.checkValidScriptParameter(specimen.getExternalId());
4599 28 Sep 17 nicklas 259         libName = libName.replace(specimen.getName(), specimen.getExternalId());
4599 28 Sep 17 nicklas 260         mergedName = mergedName.replace(specimen.getName(), specimen.getExternalId());
4599 28 Sep 17 nicklas 261       }
6640 11 Mar 22 nicklas 262       
6654 23 Mar 22 nicklas 263       String hisatRgOptions = "--rg-id " + mergedName;
6640 11 Mar 22 nicklas 264       hisatRgOptions += " --rg PU:" + Values.getString(flowCellIds, "", true);
6640 11 Mar 22 nicklas 265       hisatRgOptions += " --rg SM:" + sampleName;
6640 11 Mar 22 nicklas 266       hisatRgOptions += " --rg LB:" + libName;
4590 25 Sep 17 nicklas 267     
4599 28 Sep 17 nicklas 268       ScriptBuilder script = new ScriptBuilder();
6665 05 Apr 22 nicklas 269       script.cmd(debug ? "set -ex" : "set -e");
5596 11 Sep 19 nicklas 270       // Set file permissions based on consent or external group!
5596 11 Sep 19 nicklas 271       String externalGroup = isExternal ? Reggie.getExternalGroup(merged.getName()) : null;
5596 11 Sep 19 nicklas 272       ScriptUtil.setUmaskForItem(dc, lib, externalGroup, script);
4599 28 Sep 17 nicklas 273       script.newLine();
6693 22 Apr 22 nicklas 274       script.cmd(global_env);
6654 23 Mar 22 nicklas 275       script.export("ArchiveFolder", archiveFolder);
6654 23 Mar 22 nicklas 276       script.export("FastqFolder", "${ArchiveFolder}"+fastQFolder);
6654 23 Mar 22 nicklas 277       script.export("HisatFolder", "${ArchiveFolder}"+hisatFolder);
4599 28 Sep 17 nicklas 278       script.newLine();
6654 23 Mar 22 nicklas 279       script.cmd(hisat_env);
6665 05 Apr 22 nicklas 280       if (debug) script.cmd(hisat_envdebug);
6654 23 Mar 22 nicklas 281       script.export("HisatOptions", "${HisatOptions} "+hisatRgOptions);
6654 23 Mar 22 nicklas 282       script.newLine();
6640 11 Mar 22 nicklas 283       script.cmd(hisat_execute);
6654 23 Mar 22 nicklas 284       script.newLine();
5596 11 Sep 19 nicklas 285       if (externalGroup != null)
5596 11 Sep 19 nicklas 286       {
5930 06 May 20 nicklas 287         ScriptUtil.addChgrp(externalGroup, "${HisatFolder}", alignedName, null, script);
5596 11 Sep 19 nicklas 288       }
4599 28 Sep 17 nicklas 289       
6674 11 Apr 22 nicklas 290       JobDefinition jobDef = new JobDefinition("HisatAlign", jobConfig, batchConfig, alignJob);
6640 11 Mar 22 nicklas 291       jobDef.addFile(ScriptUtil.upload("hisat.sh"));
6640 11 Mar 22 nicklas 292       jobDef.addFile(ScriptUtil.upload("reggie-utils.sh"));
6640 11 Mar 22 nicklas 293       jobDef.addFile(ScriptUtil.upload("stdwrap.sh"));
6640 11 Mar 22 nicklas 294       jobDef.addFile(ScriptUtil.upload("alignment_statistics.sh"));
6640 11 Mar 22 nicklas 295       jobDef.addFile(ScriptUtil.upload("singlecolumnaverager.awk"), FilePermission.USER_RWX);
4599 28 Sep 17 nicklas 296       jobDef.setDebug(debug);
4599 28 Sep 17 nicklas 297       jobDef.setCmd(script.toString());
4599 28 Sep 17 nicklas 298       jobDefs.add(jobDef);
4590 25 Sep 17 nicklas 299     }
4590 25 Sep 17 nicklas 300     
4599 28 Sep 17 nicklas 301     return jobDefs;
4590 25 Sep 17 nicklas 302   }
4590 25 Sep 17 nicklas 303   
4590 25 Sep 17 nicklas 304   /**
4590 25 Sep 17 nicklas 305     Job completion handler for mask/align jobs. The handler downloads the
4590 25 Sep 17 nicklas 306     'filter.out' file from the job folder and parses out number of
4590 25 Sep 17 nicklas 307     reads remaining. The information is stored as {@link Annotationtype#READS}
4590 25 Sep 17 nicklas 308     annotation on FilteredSequences items.
4590 25 Sep 17 nicklas 309   */
4590 25 Sep 17 nicklas 310   public static class AlignJobCompletionHandler
4590 25 Sep 17 nicklas 311     implements JobCompletionHandler
4590 25 Sep 17 nicklas 312   {
7079 27 Mar 23 nicklas 313     private static final ExtensionsLogger logger = 
7079 27 Mar 23 nicklas 314       ExtensionsLog.getLogger(JobCompletionHandlerFactory.ID, true).wrap(LoggerFactory.getLogger(AlignJobCompletionHandler.class));
4590 25 Sep 17 nicklas 315     
4590 25 Sep 17 nicklas 316     public AlignJobCompletionHandler()
4590 25 Sep 17 nicklas 317     {}
4590 25 Sep 17 nicklas 318   
4590 25 Sep 17 nicklas 319     @Override
4590 25 Sep 17 nicklas 320     public String jobCompleted(SessionControl sc, OpenGridSession session, Job job, JobStatus status)
4590 25 Sep 17 nicklas 321     {
4590 25 Sep 17 nicklas 322       String jobName = status.getName();
4656 24 Jan 18 nicklas 323       String masked = session.getJobFileAsString(jobName, "masked.out", "UTF-8");
4590 25 Sep 17 nicklas 324       String alignStatistics = session.getJobFileAsString(jobName, "alignment_statistics.out", "UTF-8");
4590 25 Sep 17 nicklas 325       String picardMetrics = session.getJobFileAsString(jobName, "alignment_picardmetrics.csv", "UTF-8");
4590 25 Sep 17 nicklas 326       String fragments = session.getJobFileAsString(jobName, "fragments.out", "UTF-8");
4590 25 Sep 17 nicklas 327       String files = session.getJobFileAsString(jobName, "files.out", "UTF-8");
4590 25 Sep 17 nicklas 328
4656 24 Jan 18 nicklas 329       Metrics metrics = parseMaskedAndAlignedOut(sc, job, masked, alignStatistics, picardMetrics, fragments, files);
4676 08 Feb 18 nicklas 330       String msg = Values.formatNumber(metrics.numReadsAfterMask/1000000f, 1) + "M reads after mask; ";
4676 08 Feb 18 nicklas 331       msg += Values.formatNumber(metrics.numReadsAfterAlign/1000000f, 1) + "M reads after alignment; ";
4676 08 Feb 18 nicklas 332       msg += Values.formatNumber(metrics.fractionDuplication * 100,  1) + "% duplicates; ";
4676 08 Feb 18 nicklas 333       msg += Values.formatNumber(metrics.hetPercentage, 0) + "% HET";
4590 25 Sep 17 nicklas 334       return msg;
4590 25 Sep 17 nicklas 335     }
4590 25 Sep 17 nicklas 336     
4656 24 Jan 18 nicklas 337     private Metrics parseMaskedAndAlignedOut(SessionControl sc, Job job, String maskedOut, String alignOut, String picardMetrics, String fragments, String filesOut)
4590 25 Sep 17 nicklas 338     {
4590 25 Sep 17 nicklas 339       Metrics metrics = new Metrics();
4590 25 Sep 17 nicklas 340       
4656 24 Jan 18 nicklas 341       Pattern p = Pattern.compile("\\s+(\\d+).*aligned concordantly 0 times.*");
4656 24 Jan 18 nicklas 342       for (String line : maskedOut.split("\n"))
4656 24 Jan 18 nicklas 343       {
4656 24 Jan 18 nicklas 344         Matcher m = p.matcher(line);
4656 24 Jan 18 nicklas 345         if (m.matches())
4656 24 Jan 18 nicklas 346         {
4656 24 Jan 18 nicklas 347           metrics.numReadsAfterMask = Values.getLong(m.group(1), null);
4656 24 Jan 18 nicklas 348           if (logger.isDebugEnabled())
4656 24 Jan 18 nicklas 349           {
4656 24 Jan 18 nicklas 350             logger.debug("Found match: " + line + "; numReadsAfterMask="+metrics.numReadsAfterMask);
4656 24 Jan 18 nicklas 351           }
4656 24 Jan 18 nicklas 352           break;
4656 24 Jan 18 nicklas 353         }
4656 24 Jan 18 nicklas 354       }
4656 24 Jan 18 nicklas 355       
4656 24 Jan 18 nicklas 356       p = Pattern.compile("(\\d+)");
4590 25 Sep 17 nicklas 357       for (String line : alignOut.split("\n"))
4590 25 Sep 17 nicklas 358       {
4590 25 Sep 17 nicklas 359         Matcher m = p.matcher(line);
4590 25 Sep 17 nicklas 360         if (m.matches())
4590 25 Sep 17 nicklas 361         {
4590 25 Sep 17 nicklas 362           metrics.numReadsAfterAlign = Values.getLong(m.group(1), null);
4590 25 Sep 17 nicklas 363           if (logger.isDebugEnabled())
4590 25 Sep 17 nicklas 364           {
4590 25 Sep 17 nicklas 365             logger.debug("Found match: " + line + "; numReadsAfterAlign="+metrics.numReadsAfterAlign);
4590 25 Sep 17 nicklas 366           }
4590 25 Sep 17 nicklas 367           break;
4590 25 Sep 17 nicklas 368         }        
4590 25 Sep 17 nicklas 369       }
4590 25 Sep 17 nicklas 370       
4590 25 Sep 17 nicklas 371       int readPairsExaminedIndex = -1;
4590 25 Sep 17 nicklas 372       int readPairDuplicatesIndex = -1;
4590 25 Sep 17 nicklas 373       int percentDuplicationIndex = -1;
4590 25 Sep 17 nicklas 374       
4590 25 Sep 17 nicklas 375       for (String line : picardMetrics.split("\n"))
4590 25 Sep 17 nicklas 376       {
4590 25 Sep 17 nicklas 377         String[] cols = line.split("\t");
4590 25 Sep 17 nicklas 378         if (cols.length >= 9)
4590 25 Sep 17 nicklas 379         {
4590 25 Sep 17 nicklas 380           if (readPairsExaminedIndex == -1)
4590 25 Sep 17 nicklas 381           {
4590 25 Sep 17 nicklas 382             List<String> colsA = Arrays.asList(cols);
4590 25 Sep 17 nicklas 383             readPairsExaminedIndex = colsA.indexOf("READ_PAIRS_EXAMINED");
4590 25 Sep 17 nicklas 384             readPairDuplicatesIndex = colsA.indexOf("READ_PAIR_DUPLICATES");
4590 25 Sep 17 nicklas 385             percentDuplicationIndex = colsA.indexOf("PERCENT_DUPLICATION");
4590 25 Sep 17 nicklas 386           }
4590 25 Sep 17 nicklas 387           else
4590 25 Sep 17 nicklas 388           {
4590 25 Sep 17 nicklas 389             metrics.readPairsExamined = Values.getLong(cols[readPairsExaminedIndex], null);
4590 25 Sep 17 nicklas 390             metrics.readPairDuplicates = Values.getLong(cols[readPairDuplicatesIndex], null);
4590 25 Sep 17 nicklas 391             metrics.fractionDuplication = Values.getFloat(cols[percentDuplicationIndex], null);
4590 25 Sep 17 nicklas 392           }
4590 25 Sep 17 nicklas 393         }
4590 25 Sep 17 nicklas 394       }
4590 25 Sep 17 nicklas 395       
4590 25 Sep 17 nicklas 396       // Fragments
4590 25 Sep 17 nicklas 397       p = Pattern.compile("(\\d+)\\t(\\d+\\.?\\d*)\\t(\\d+\\.?\\d*)");
4590 25 Sep 17 nicklas 398       for (String line : fragments.split("\n"))
4590 25 Sep 17 nicklas 399       {
4590 25 Sep 17 nicklas 400         Matcher m = p.matcher(line);
4590 25 Sep 17 nicklas 401         if (m.matches())
4590 25 Sep 17 nicklas 402         {
4590 25 Sep 17 nicklas 403           metrics.fragmentSizeCount = Values.getInt(m.group(1), -1);
4590 25 Sep 17 nicklas 404           metrics.fragmentSizeAvg = Values.getInt(m.group(2), -1);
4590 25 Sep 17 nicklas 405           metrics.fragmentSizeStd = Values.getInt(m.group(3), -1);
4590 25 Sep 17 nicklas 406         }
4590 25 Sep 17 nicklas 407       }
4590 25 Sep 17 nicklas 408       
4590 25 Sep 17 nicklas 409       DbControl dc = null;
4590 25 Sep 17 nicklas 410       try
4590 25 Sep 17 nicklas 411       {
6599 22 Feb 22 nicklas 412         dc = sc.newDbControl("Reggie: Hisat alignment completed handler");
4590 25 Sep 17 nicklas 413         
4590 25 Sep 17 nicklas 414         AlignedSequences alignedSequences = AlignedSequences.getByJob(dc, job);
4656 24 Jan 18 nicklas 415         MaskedSequences maskedSequences = alignedSequences.getMaskedSequences(dc);
4656 24 Jan 18 nicklas 416
4590 25 Sep 17 nicklas 417         DerivedBioAssay aligned = alignedSequences.getItem();
4656 24 Jan 18 nicklas 418         DerivedBioAssay masked = maskedSequences.getItem();
4590 25 Sep 17 nicklas 419
4656 24 Jan 18 nicklas 420         Annotationtype.PM_READS.setAnnotationValue(dc, masked, metrics.numReadsAfterMask);
4590 25 Sep 17 nicklas 421         Annotationtype.ALIGNED_PAIRS.setAnnotationValue(dc, aligned, metrics.numReadsAfterAlign);
4590 25 Sep 17 nicklas 422         Annotationtype.READ_PAIRS_EXAMINED.setAnnotationValue(dc, aligned, metrics.readPairsExamined);
4590 25 Sep 17 nicklas 423         Annotationtype.READ_PAIR_DUPLICATES.setAnnotationValue(dc, aligned, metrics.readPairDuplicates);
4590 25 Sep 17 nicklas 424         Annotationtype.FRACTION_DUPLICATION.setAnnotationValue(dc, aligned, metrics.fractionDuplication);
4590 25 Sep 17 nicklas 425         Annotationtype.FRAGMENT_SIZE_AVG.setAnnotationValue(dc, aligned, metrics.fragmentSizeAvg);
4590 25 Sep 17 nicklas 426         Annotationtype.FRAGMENT_SIZE_STDEV.setAnnotationValue(dc, aligned, metrics.fragmentSizeStd);
4590 25 Sep 17 nicklas 427   
4590 25 Sep 17 nicklas 428         // Create file links
5553 12 Aug 19 nicklas 429         boolean useExternalProjectArchive = Reggie.isExternalItem(aligned.getName());
4590 25 Sep 17 nicklas 430         FileServer fileArchive = useExternalProjectArchive ? Fileserver.EXTERNAL_ARCHIVE.load(dc) : Fileserver.PROJECT_ARCHIVE.load(dc);
4590 25 Sep 17 nicklas 431         String analysisDir = useExternalProjectArchive ? Reggie.EXTERNAL_ANALYSIS_DIR : Reggie.SECONDARY_ANALYSIS_DIR;
4590 25 Sep 17 nicklas 432
4590 25 Sep 17 nicklas 433         String dataFilesFolder = (String)Annotationtype.DATA_FILES_FOLDER.getAnnotationValue(dc, aligned);
4590 25 Sep 17 nicklas 434         String baseFolder = Reggie.convertDataFilesFolderToBaseFolder(dataFilesFolder);
4590 25 Sep 17 nicklas 435         Directory localDataDir = Directory.getNew(dc, new Path(analysisDir+baseFolder, Path.Type.DIRECTORY));
4590 25 Sep 17 nicklas 436         DataFileType bamData = Datafiletype.BAM.load(dc);
4590 25 Sep 17 nicklas 437         ItemSubtype bamType = bamData.getGenericType();
4619 13 Nov 17 nicklas 438         ItemSubtype vcfType = Subtype.VARIANT_CALL_FORMAT.load(dc);
4590 25 Sep 17 nicklas 439   
4590 25 Sep 17 nicklas 440         int lineNo = 0;
4590 25 Sep 17 nicklas 441         for (String line : filesOut.split("\n"))
4590 25 Sep 17 nicklas 442         {
4590 25 Sep 17 nicklas 443           lineNo++;
4590 25 Sep 17 nicklas 444           
4590 25 Sep 17 nicklas 445           File f = File.getFile(dc, localDataDir, line.substring(line.lastIndexOf("/")+1), true);
4590 25 Sep 17 nicklas 446           f.setFileServer(fileArchive);
4590 25 Sep 17 nicklas 447           String fileUrl = "sftp://" + fileArchive.getHost() + dataFilesFolder + "/" + f.getName();
4590 25 Sep 17 nicklas 448           try
4590 25 Sep 17 nicklas 449           {
4590 25 Sep 17 nicklas 450             f.setUrl(fileUrl, true);
4590 25 Sep 17 nicklas 451           }
4590 25 Sep 17 nicklas 452           catch (RuntimeException ex)
4590 25 Sep 17 nicklas 453           {
4590 25 Sep 17 nicklas 454             f.setUrl(fileUrl, false);
4590 25 Sep 17 nicklas 455           }
4619 13 Nov 17 nicklas 456           
4590 25 Sep 17 nicklas 457           if (!f.isInDatabase())
4590 25 Sep 17 nicklas 458           {
4590 25 Sep 17 nicklas 459             dc.saveItem(f);
4590 25 Sep 17 nicklas 460           }
4590 25 Sep 17 nicklas 461           if (f.getName().equals("alignment.bam"))
4590 25 Sep 17 nicklas 462           {
4590 25 Sep 17 nicklas 463             f.setDescription(metrics.numReadsAfterAlign + " ALIGNED PAIRS");
4590 25 Sep 17 nicklas 464             f.setItemSubtype(bamType);
4590 25 Sep 17 nicklas 465             FileSetMember member = aligned.getFileSet().addMember(f, bamData);
4590 25 Sep 17 nicklas 466           }
4590 25 Sep 17 nicklas 467           else
4590 25 Sep 17 nicklas 468           {
4590 25 Sep 17 nicklas 469             AnyToAny link = AnyToAny.getNewOrExisting(dc, aligned, f.getName(), f, true);
4590 25 Sep 17 nicklas 470             if (!link.isInDatabase()) dc.saveItem(link);
4619 13 Nov 17 nicklas 471             
4619 13 Nov 17 nicklas 472             // We copy the VCF file to the BASE server
4622 16 Nov 17 nicklas 473             if ("qc_genotype.vcf".equals(f.getName()))
4619 13 Nov 17 nicklas 474             {
4619 13 Nov 17 nicklas 475               f.setMimeTypeAuto("text/plain", vcfType);
5579 20 Aug 19 nicklas 476               VcfData vcfData = copyToBaseAndParse(f, aligned);
4622 16 Nov 17 nicklas 477               if (vcfData != null)
4622 16 Nov 17 nicklas 478               {
4622 16 Nov 17 nicklas 479                 metrics.genotypeCount = vcfData.getGtCount();
4622 16 Nov 17 nicklas 480                 metrics.hetCount = vcfData.getHetCount();
4676 08 Feb 18 nicklas 481                 metrics.hetPercentage = vcfData.getHetPercentage();
4622 16 Nov 17 nicklas 482                 Annotationtype.QC_GENOTYPE_COUNT.setAnnotationValue(dc, aligned, metrics.genotypeCount);
4658 26 Jan 18 nicklas 483                 if (metrics.genotypeCount > 0)
4658 26 Jan 18 nicklas 484                 {
4676 08 Feb 18 nicklas 485                   Annotationtype.QC_GENOTYPE_HET_PCT.setAnnotationValue(dc, aligned, metrics.hetPercentage);
4658 26 Jan 18 nicklas 486                 }
4688 27 Feb 18 nicklas 487                 if (useExternalProjectArchive)
4688 27 Feb 18 nicklas 488                 {
4688 27 Feb 18 nicklas 489                   Annotationtype.QC_GENOTYPE_STATUS.setAnnotationValue(dc, aligned, "Disabled");
4688 27 Feb 18 nicklas 490                 }
4622 16 Nov 17 nicklas 491               }
4619 13 Nov 17 nicklas 492             }
4590 25 Sep 17 nicklas 493           }
4590 25 Sep 17 nicklas 494         }
4590 25 Sep 17 nicklas 495         
4590 25 Sep 17 nicklas 496         dc.commit();
4590 25 Sep 17 nicklas 497       }
4590 25 Sep 17 nicklas 498       finally
4590 25 Sep 17 nicklas 499       {
4590 25 Sep 17 nicklas 500         if (dc != null) dc.close();
4590 25 Sep 17 nicklas 501       }
4590 25 Sep 17 nicklas 502   
4590 25 Sep 17 nicklas 503       return metrics;
4590 25 Sep 17 nicklas 504     }
4622 16 Nov 17 nicklas 505
4622 16 Nov 17 nicklas 506     /**
4622 16 Nov 17 nicklas 507       Helper method for copying the VCF file from the file server
4622 16 Nov 17 nicklas 508       while at the same time parsing it and extracting genotype
4622 16 Nov 17 nicklas 509       information and statistics.
4622 16 Nov 17 nicklas 510     */
5579 20 Aug 19 nicklas 511     private VcfData copyToBaseAndParse(File vcfFile, DerivedBioAssay alignment) 
4622 16 Nov 17 nicklas 512     {
4622 16 Nov 17 nicklas 513       // Stream for copying the vcfFile
4622 16 Nov 17 nicklas 514       InputStream fromFileServer = null;
4622 16 Nov 17 nicklas 515       OutputStream toBase = null;
4622 16 Nov 17 nicklas 516       
4622 16 Nov 17 nicklas 517       // The splitter will help us parse the file and at the same time copy it to BASE.
4622 16 Nov 17 nicklas 518       InputStreamSplitter splitter = null;
4622 16 Nov 17 nicklas 519       
4622 16 Nov 17 nicklas 520       VcfData vcfData = null;
4622 16 Nov 17 nicklas 521       try
4622 16 Nov 17 nicklas 522       {
4622 16 Nov 17 nicklas 523         fromFileServer = vcfFile.getDownloadStream(0);
4622 16 Nov 17 nicklas 524         toBase = vcfFile.getUploadStream(false);
4622 16 Nov 17 nicklas 525         splitter = new InputStreamSplitter(fromFileServer, true, true, toBase);
4622 16 Nov 17 nicklas 526         
4622 16 Nov 17 nicklas 527         VcfParser parser = new VcfParser();
4622 16 Nov 17 nicklas 528         vcfData = parser.parse(splitter, vcfFile.getName());
4622 16 Nov 17 nicklas 529         vcfFile.setFileServer(null);
4622 16 Nov 17 nicklas 530       }
4622 16 Nov 17 nicklas 531       catch (Exception ex)
4622 16 Nov 17 nicklas 532       {
5579 20 Aug 19 nicklas 533         throw new BaseException("Could not parse '" + vcfFile.getName() + "' for alignment: " + alignment.getName(), ex);
4622 16 Nov 17 nicklas 534       }
4622 16 Nov 17 nicklas 535       finally
4622 16 Nov 17 nicklas 536       {
4622 16 Nov 17 nicklas 537         FileUtil.close(splitter);
4622 16 Nov 17 nicklas 538         FileUtil.close(toBase);
4622 16 Nov 17 nicklas 539         FileUtil.close(fromFileServer);
4622 16 Nov 17 nicklas 540       }
4622 16 Nov 17 nicklas 541       return vcfData;
4622 16 Nov 17 nicklas 542     }
4622 16 Nov 17 nicklas 543
4590 25 Sep 17 nicklas 544   }
4590 25 Sep 17 nicklas 545
4622 16 Nov 17 nicklas 546   
4590 25 Sep 17 nicklas 547   static class Metrics
4590 25 Sep 17 nicklas 548   {
4656 24 Jan 18 nicklas 549     Long numReadsAfterMask = null;
4590 25 Sep 17 nicklas 550     Long numReadsAfterAlign = null;
4590 25 Sep 17 nicklas 551     Long readPairsExamined = null;
4590 25 Sep 17 nicklas 552     Long readPairDuplicates = null;
4590 25 Sep 17 nicklas 553     Float fractionDuplication = null;
4590 25 Sep 17 nicklas 554     int fragmentSizeAvg = -1;
4590 25 Sep 17 nicklas 555     int fragmentSizeStd = -1;
4590 25 Sep 17 nicklas 556     int fragmentSizeCount = -1;
4622 16 Nov 17 nicklas 557     int genotypeCount = -1;
4622 16 Nov 17 nicklas 558     int hetCount = -1;
4676 08 Feb 18 nicklas 559     float hetPercentage = Float.NaN;
4590 25 Sep 17 nicklas 560   }
4590 25 Sep 17 nicklas 561 }

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