Changeset 5831

Timestamp:

Feb 21, 2020, 2:36:41 PM (3 years ago)

Author:

Nicklas Nordborg

Message:

References #1218: Implement MIPs alignment

Final steps in the alignment pipeline. Re-headers, splitting and marking duplicates. All critical work has been done. There are still some collections of metrics to implement, as well as error handling and storing result files back to the project archive.

Location:

extensions/net.sf.basedb.reggie/trunk

Files:

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

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

@@ -361 +361 @@
         </novoalign>
         
+        <genome-dict>${ReferenceDir}/gatk_bundle/2.8/b37/human_g1k_v37_decoy.dict</genome-dict>
+        <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>
+        
       </align-mips>
       

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

@@ -148 +148 @@
     String novoalign_path = cfg.getRequiredConfig("programs/novoalign/path", alignParameterSet);
     String novosort_path = cfg.getRequiredConfig("programs/novosort/path", alignParameterSet);
+    String samtools_path = cfg.getRequiredConfig("programs/samtools/path", alignParameterSet);
 
     // Options for the programs
@@ -156 +157 @@
     String novoalign_index = cfg.getRequiredConfig("align-mips/novoalign/index", alignParameterSet);
     String novoalign_options = cfg.getRequiredConfig("align-mips/novoalign/options", alignParameterSet);
+    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);
 
     // Load common items
@@ -278 +282 @@
       script.cmd("NovoAlign="+novoalign_path);
       script.cmd("NovoSort="+novosort_path);
+      script.cmd("samtools="+samtools_path);
       script.cmd("ReferenceDir=" + referenceRoot);
       script.cmd("AmpliconsDir=" + ampliconsRoot);
       script.cmd("AmpliconsBed=${AmpliconsDir}/"+amplicons_bed);
       script.cmd("NovoIndex="+novoalign_index);
+      script.cmd("GenomeDict="+genome_dict);
+      script.cmd("GenomeFasta="+genome_fasta);
       script.cmd("ArchiveRoot="+archiveRoot);
       script.cmd("FastqFolder=${ArchiveRoot}"+fastQFolder);
@@ -309 +316 @@
       script.cmd("mkdir bam.sorted");
       script.cmd("mkdir bam.merged");
+      script.cmd("mkdir bam.final");
       script.cmd("cp ${ScriptDir}/stdwrap.sh .");
       script.cmd("cp ${ScriptDir}/stderrwrap.sh .");
@@ -428 +436 @@
       mergeCmd += " > mergesam.out";
       script.cmd(mergeCmd);
-
+      script.newLine();
+
+      script.comment("Re-header BAM files");
+      script.cmd("${samtools} view -H bam.merged/novo.bam | grep -v \"^@SQ\" > header.sam");
+      script.cmd("grep \"^@SQ\" ${GenomeDict} >> header.sam");
+      script.cmd("${samtools} reheader -P header.sam bam.merged/novo.bam > bam.merged/novo_reheaded.bam");
+      script.newLine();
+
+      script.comment("Separate concordant, discordant and unmapped read pairs");
+      script.cmd("${samtools} view -b -h -@ 2 -f 3 bam.merged/novo_reheaded.bam > bam.merged/concordant.bam");
+      script.cmd("${samtools} view -b -h -@ 2 -G 12 -F 2 bam.merged/novo_reheaded.bam > bam.merged/discordant.bam");
+      script.cmd("${samtools} view -b -h -@ 2 -f 12 bam.merged/novo_reheaded.bam > bam.merged/unmapped.bam");
+      
+      script.comment("Picard UmiAwareMarkDuplicatesWithMateCigar");
+      String markDupCmd = "./stdwrap.sh ./picard2 UmiAwareMarkDuplicatesWithMateCigar";
+      markDupCmd += " " + mark_duplicate_options;
+      markDupCmd += " -UMI_METRICS bam.final/concordant.umi_metrics.txt";
+      markDupCmd += " -METRICS_FILE bam.final/concordant.dedup_metrics.txt";
+      markDupCmd += " -INPUT bam.merged/concordant.bam";
+      markDupCmd += " -OUTPUT bam.final/concordant.bam";
+      markDupCmd += " >> markduplicates.out";
+      script.cmd(markDupCmd);
+      
+      markDupCmd = "./stdwrap.sh ./picard2 UmiAwareMarkDuplicatesWithMateCigar";
+      markDupCmd += " " + mark_duplicate_options;
+      markDupCmd += " -UMI_METRICS bam.final/discordant.umi_metrics.txt";
+      markDupCmd += " -METRICS_FILE bam.final/discordant.dedup_metrics.txt";
+      markDupCmd += " -INPUT bam.merged/discordant.bam";
+      markDupCmd += " -OUTPUT bam.final/discordant.bam";
+      markDupCmd += " >> markduplicates.out";
+      script.cmd(markDupCmd);
+      
       /*
       script.progress(95, "Copying result files to project archive");