perl/etc/WoodyMiaoLin/PbeBefore2015/addRG.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 die "\nAdd HEADER and TAG of READ GROUP to bam file.\nAuther: Woody\nUsage: $0 <in.bam> <out.bam>\n\n" if @ARGV<2;
   6
   7 my $in = shift;
   8 my $out = shift;
   9 $in =~ /^(.+).bam$/;
  10 my $ip = $1;
  11 $out =~ /^(.+).bam$/;
  12 my $op = $1;
  13
  14 open I, "-|", "samtools view -h $in";
  15 open HEAD, ">", "$op.header.sam";
  16 open READ, "|-", "samtools view -bS - >$op.reads.bam";
  17
  18 warn "Adding TAG...\n";
  19 my %RG;
  20 while (<I>) {
  21         if (/^@/) {
  22                 print HEAD;
  23                 print READ;
  24         } else {
  25                 chomp;
  26                 my @a = split /\t/;
  27                 my @b = split /:/, $a[0];
  28                 my $rg = "${ip}_$b[2]_$b[3]";
  29                 $RG{$rg} = 1;
  30                 print READ "$_\tRG:Z:$rg\n";
  31         }
  32 }
  33
  34 close I;
  35 close READ;
  36
  37 foreach (sort keys %RG) {
  38         print HEAD "\@RG\tID:$_\tPL:ILLUMINA\tSM:$ip\n";
  39 }
  40
  41 warn "Printing HEADER...\n";
  42 my $time = localtime();
  43 print HEAD "\@CO\tHEADER and TAG of READ GROUP added at $time. ID: FilePre_FlowId_LaneNum, SM: FilePre.\n";
  44 close HEAD;
  45
  46 warn "Reheadering...\n";
  47 system "samtools reheader $op.header.sam $op.reads.bam >$op.bam";
  48
  49 warn "Removing $op.header.sam $op.reads.bam...\n";
  50 system "rm $op.header.sam $op.reads.bam";
  51 warn "Done!\n";