Changeset 5831


Ignore:
Timestamp:
Feb 21, 2020, 2:36:41 PM (22 months 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:
2 edited

Legend:

Unmodified
Added
Removed
  • extensions/net.sf.basedb.reggie/trunk/config/reggie-config.xml

    r5830 r5831  
    361361        </novoalign>
    362362       
     363        <genome-dict>${ReferenceDir}/gatk_bundle/2.8/b37/human_g1k_v37_decoy.dict</genome-dict>
     364        <genome-fasta>${ReferenceDir}/gatk_bundle/2.8/b37/human_g1k_v37_decoy.fasta</genome-fasta>
     365       
     366        <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>
     367       
    363368      </align-mips>
    364369     
  • extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/MipsAlignJobCreator.java

    r5830 r5831  
    148148    String novoalign_path = cfg.getRequiredConfig("programs/novoalign/path", alignParameterSet);
    149149    String novosort_path = cfg.getRequiredConfig("programs/novosort/path", alignParameterSet);
     150    String samtools_path = cfg.getRequiredConfig("programs/samtools/path", alignParameterSet);
    150151
    151152    // Options for the programs
     
    156157    String novoalign_index = cfg.getRequiredConfig("align-mips/novoalign/index", alignParameterSet);
    157158    String novoalign_options = cfg.getRequiredConfig("align-mips/novoalign/options", alignParameterSet);
     159    String genome_dict = cfg.getRequiredConfig("align-mips/genome-dict", alignParameterSet);
     160    String genome_fasta = cfg.getRequiredConfig("align-mips/genome-fasta", alignParameterSet);
     161    String mark_duplicate_options = cfg.getRequiredConfig("align-mips/mark-duplicate-options", alignParameterSet);
    158162
    159163    // Load common items
     
    278282      script.cmd("NovoAlign="+novoalign_path);
    279283      script.cmd("NovoSort="+novosort_path);
     284      script.cmd("samtools="+samtools_path);
    280285      script.cmd("ReferenceDir=" + referenceRoot);
    281286      script.cmd("AmpliconsDir=" + ampliconsRoot);
    282287      script.cmd("AmpliconsBed=${AmpliconsDir}/"+amplicons_bed);
    283288      script.cmd("NovoIndex="+novoalign_index);
     289      script.cmd("GenomeDict="+genome_dict);
     290      script.cmd("GenomeFasta="+genome_fasta);
    284291      script.cmd("ArchiveRoot="+archiveRoot);
    285292      script.cmd("FastqFolder=${ArchiveRoot}"+fastQFolder);
     
    309316      script.cmd("mkdir bam.sorted");
    310317      script.cmd("mkdir bam.merged");
     318      script.cmd("mkdir bam.final");
    311319      script.cmd("cp ${ScriptDir}/stdwrap.sh .");
    312320      script.cmd("cp ${ScriptDir}/stderrwrap.sh .");
     
    428436      mergeCmd += " > mergesam.out";
    429437      script.cmd(mergeCmd);
    430 
     438      script.newLine();
     439
     440      script.comment("Re-header BAM files");
     441      script.cmd("${samtools} view -H bam.merged/novo.bam | grep -v \"^@SQ\" > header.sam");
     442      script.cmd("grep \"^@SQ\" ${GenomeDict} >> header.sam");
     443      script.cmd("${samtools} reheader -P header.sam bam.merged/novo.bam > bam.merged/novo_reheaded.bam");
     444      script.newLine();
     445
     446      script.comment("Separate concordant, discordant and unmapped read pairs");
     447      script.cmd("${samtools} view -b -h -@ 2 -f 3 bam.merged/novo_reheaded.bam > bam.merged/concordant.bam");
     448      script.cmd("${samtools} view -b -h -@ 2 -G 12 -F 2 bam.merged/novo_reheaded.bam > bam.merged/discordant.bam");
     449      script.cmd("${samtools} view -b -h -@ 2 -f 12 bam.merged/novo_reheaded.bam > bam.merged/unmapped.bam");
     450     
     451      script.comment("Picard UmiAwareMarkDuplicatesWithMateCigar");
     452      String markDupCmd = "./stdwrap.sh ./picard2 UmiAwareMarkDuplicatesWithMateCigar";
     453      markDupCmd += " " + mark_duplicate_options;
     454      markDupCmd += " -UMI_METRICS bam.final/concordant.umi_metrics.txt";
     455      markDupCmd += " -METRICS_FILE bam.final/concordant.dedup_metrics.txt";
     456      markDupCmd += " -INPUT bam.merged/concordant.bam";
     457      markDupCmd += " -OUTPUT bam.final/concordant.bam";
     458      markDupCmd += " >> markduplicates.out";
     459      script.cmd(markDupCmd);
     460     
     461      markDupCmd = "./stdwrap.sh ./picard2 UmiAwareMarkDuplicatesWithMateCigar";
     462      markDupCmd += " " + mark_duplicate_options;
     463      markDupCmd += " -UMI_METRICS bam.final/discordant.umi_metrics.txt";
     464      markDupCmd += " -METRICS_FILE bam.final/discordant.dedup_metrics.txt";
     465      markDupCmd += " -INPUT bam.merged/discordant.bam";
     466      markDupCmd += " -OUTPUT bam.final/discordant.bam";
     467      markDupCmd += " >> markduplicates.out";
     468      script.cmd(markDupCmd);
     469     
    431470      /*
    432471      script.progress(95, "Copying result files to project archive");
Note: See TracChangeset for help on using the changeset viewer.