Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • insert PG and CL tags into BAM header with BWA

    What do you recommend I use to insert full command line into the BAM header created with BWA. BWA already creates this:

    @PG ID:bwa PN:bwa VN:0.5.9-r16

    But I'd like to capture the whole process (aln and sampe) and insert them as CL tags.

    Tophat does it nicely like this for example:

    @PG ID:TopHat VN:1.2.0 CL:/packages/tophat-1.2.0/tophat -r 43 -p 8 -o /some/path/s_1_topHat -G /some/path/Homo_sapiens.NCBI36.54.noMT.gtf /some/path/human_b36_both_noM /some/path/s_1_1.fastq /some/path/s_1_3.fastq

  • #2
    One solution is a script in your language of choice (Perl, Python, etc) to edit the current header, combined with samtools reheader to update the BAM file.

    Comment


    • #3
      First way:

      samtools view -H file.bam > header.sam

      edit header.sam file

      samtools reheader header.sam file.bam > file_new_head.bam

      OR as a one liner, something like this:

      samtools view -h file.bam | perl -pe 's/^\@PG\told line$/@PG\tnew line/' | samtools view -Sb -o file_new_head.bam -

      Chris

      Comment


      • #4
        adding @PG @RG and more with BamUtil polishBam

        I disc-overed this great 'Umich' tooolbox after searching for quite some time (http://genome.sph.umich.edu/wiki/BamUtil:_polishBam)

        I first stored my full bwa command in a variable 'cmd' then looked for my bwa version using 'expr' and finally built the full @PG string taking care of keeping tabs (that was hard! ).

        The result was accepted by polishBam and creates a true @PG line.
        I guess similar approach will work for other header tags.

        Code:
        rgstring='@RG\tID:NA18507\tLB:lib-NA18507\tPU:unkn-0.0\tPL:ILLUMINA\tSM:GAIIx-chr21-BWA.mem'
        
        ... more parameters not replicated here (obvious)
        
        cmd="bwa mem -R \"${rgstring}\" \
        	-M \
        	-t ${nthr} \
        	${refgen} \
        	${fq1} ${fq2} > ${outfolder}/${bwape}.sam \
        	2>${outfolder}/${pfx}_bwa_mem.err"
         
        # execute the command
        eval ${cmd}
         
        #### add @PG group to results with **bamUtil**
        # important formatting issues here
        # $'\t' is adding a 'tab'
        # \' surrounding ${cmd} \' make sure the long text is not interpreted
        # finally, $pgstring should be quoted
        bwaver=$(expr "$(bwa 2>&1)" : '.*Version:\(.*\)Contact.*')
        pgstring=@PG$'\t'ID:BWA$'\t'VN:${bwaver}$'\t'CL:\'${cmd}\'
         
        bam polishBam --PG "${pgstring}" \
        	-i ${outfolder}/${bwape}.sam \
        	-o ${outfolder}/${bwape}pg.sam \
        	-l ${outfolder}/polishBam_${infile%%.bam}.log
        Hope this will help others!
        thanks to the BamUtil team for these great tools



        PS: the quick and dirty method consisting in adding the @PG line with the @RG (-R) option seems to work too.

        # substitute the following string in the upper command
        rgstring='@RG\tID:NA18507\tLB:lib-NA18507\tPU:unkn-0.0\tPL:ILLUMINA\tSM:GAIIx-chr21-BWA.mem\n@PG\tID:BWA\tVN:...'

        still, I would prefer an additional BWA parameter to do this and more about headers.
        http://www.bits.vib.be/index.php

        Comment

        Latest Articles

        Collapse

        • seqadmin
          Exploring the Dynamics of the Tumor Microenvironment
          by seqadmin




          The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
          07-08-2024, 03:19 PM
        • seqadmin
          Exploring Human Diversity Through Large-Scale Omics
          by seqadmin


          In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
          06-25-2024, 06:43 AM

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by seqadmin, 07-10-2024, 07:30 AM
        0 responses
        25 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 07-03-2024, 09:45 AM
        0 responses
        201 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 07-03-2024, 08:54 AM
        0 responses
        211 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 07-02-2024, 03:00 PM
        0 responses
        193 views
        0 likes
        Last Post seqadmin  
        Working...
        X