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

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

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