perl/etc/justonce/sperm/sperm_1_startplot.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4 use IO::Unread qw(unread);
   5 use Data::Dump qw(ddx);
   6
   7 die "Usage: $0 <input sam.gz/bam> <output>\n" if @ARGV < 2;
   8 my ($inf,$outf)=@ARGV;
   9
  10 =pod
  11 open CHR,'<','chr.lst' or die "[x]Cannot read chr.lst, use `samtools faidx` and append ChrID manually to get one!\n";
  12 my %ChrGI2ID;
  13 while (<CHR>) {
  14         my ($gi,$id) = split /\s+/,$_;
  15         $ChrGI2ID{$gi}=$id;
  16 }
  17 close CHR;
  18 $ChrGI2ID{'='}='=';
  19 ddx \%ChrGI2ID;
  20 =cut
  21
  22 my $FH;
  23 if ($inf =~ /\.sam\.gz$/i) {
  24         open $FH,'-|',"gzip -dc $inf" or die "Error opening $inf: $!\n";
  25 } elsif ($inf =~ /\.bam$/i) {
  26         open $FH,'-|',"samtools view -h $inf" or die "Error opening $inf: $!\n";
  27 } elsif ($inf =~ /\.sam$/i) {
  28         open $FH,'<',$inf or die "Error opening $inf: $!\n";
  29 } else { die; }
  30
  31 my $line;
  32 my %ChrLen;
  33 while ($line = <$FH>) {
  34         unless ($line =~ /^\@/) {
  35                 unread $FH, $line;
  36                 last;
  37         }
  38         if ($line =~ /^\@SQ/) {
  39                 my (undef,$gi,$len) = split /\s+/,$line;
  40                 $gi =~ s/^SN://;
  41                 $len =~ s/^LN://;
  42                 $ChrLen{$gi} = $len;
  43         }
  44 }
  45 ddx \%ChrLen;
  46
  47 my %DatbyChrM;
  48 #$DatbyChrM{$_}=[] for values %ChrGI2ID;
  49 while ($line = <$FH>) {
  50         my ($id, $flag, $ref, $pos, $mapq, $CIAGR, $mref, $mpos, $isize, $seq, $qual, @OPT) = split /\t/,$line;
  51         #print "$id, $flag, Chr$ChrGI2ID{$ref}, $pos, $mapq, $CIAGR, Chr$ChrGI2ID{$mref}, $mpos, $isize\n";
  52         #my $commonChrID = $ChrGI2ID{$ref};
  53         my $commonChrID = $ref;
  54         my $posto1m = int($pos/1000000);
  55         ++$DatbyChrM{$commonChrID}->[$posto1m];
  56 }
  57 #ddx \%DatbyChrM;
  58
  59 close $FH;
  60
  61 open OUT,'>',$outf or die "Error opening $outf: $!\n";
  62 print OUT "# $inf\n";
  63 for my $ccid (sort { "$a$b"=~/^\d+$/ ? $a<=>$b : $a cmp $b } keys %DatbyChrM) {
  64         print OUT "[$ccid]\n";
  65         my $ArrayRef = $DatbyChrM{$ccid};
  66         for my $i ( 0 .. $#$ArrayRef ) {
  67                 my $t = $$ArrayRef[$i];
  68                 print OUT (defined $t)?$t:0,"\t",1+$i,"\n";
  69         }
  70 }
  71
  72 close OUT;
  73
  74 __END__
  75 perl startplot.pl t.sam t.out
  76 grep \\[ t.out
  77
  78 perl sperm_1_startplot.n.pl Blood-MAL.sort.bam Blood-MAL.nstat
  79 perl sperm_1_startplot.n.pl Blood-MDA.sort.bam Blood-MDA.nstat
  80 perl sperm_1_startplot.n.pl Sperm23-MDA.sort.bam Sperm23-MDA.nstat
  81 perl sperm_1_startplot.n.pl Sperm24-MDA.sort.bam Sperm24-MDA.nstat
  82 perl sperm_1_startplot.n.pl Sperm28-MDA.sort.bam Sperm28-MDA.nstat
  83 perl sperm_1_startplot.n.pl SpermS01.sort.bam SpermS01.nstat
  84 perl sperm_1_startplot.n.pl SpermS02.sort.bam SpermS02.nstat
  85 perl sperm_1_startplot.n.pl SpermS03.sort.bam SpermS03.nstat
  86
  87 $ll -L *.bam *.nstat
  88 -rw-r--r-- 1 gaoshengjie bc_tumor  31K Sep  2 00:08 Blood-MAL.nstat
  89 -rw-r--r-- 1 gaoshengjie bc_tumor 131G Aug 12 04:37 Blood-MAL.sort.bam
  90 -rw-r--r-- 1 gaoshengjie bc_tumor  31K Sep  2 02:50 Blood-MDA.nstat
  91 -rw-r--r-- 1 gaoshengjie bc_tumor  56G Aug 12 18:30 Blood-MDA.sort.bam
  92 -rw-r--r-- 1 gaoshengjie bc_tumor  29K Sep  2 03:24 Sperm23-MDA.nstat
  93 -rw-r--r-- 1 gaoshengjie bc_tumor  20G Aug  8 11:52 Sperm23-MDA.sort.bam
  94 -rw-r--r-- 1 gaoshengjie bc_tumor  29K Sep  2 03:56 Sperm24-MDA.nstat
  95 -rw-r--r-- 1 gaoshengjie bc_tumor  20G Aug  8 12:10 Sperm24-MDA.sort.bam
  96 -rw-r--r-- 1 gaoshengjie bc_tumor  29K Sep  2 04:30 Sperm28-MDA.nstat
  97 -rw-r--r-- 1 gaoshengjie bc_tumor  20G Aug  8 12:02 Sperm28-MDA.sort.bam
  98 -rw-r--r-- 1 gaoshengjie bc_tumor  29K Sep  2 05:03 SpermS01.nstat
  99 -rw-r--r-- 1 gaoshengjie bc_tumor  22G Aug 11 13:06 SpermS01.sort.bam
 100 -rw-r--r-- 1 gaoshengjie bc_tumor  29K Sep  2 05:36 SpermS02.nstat
 101 -rw-r--r-- 1 gaoshengjie bc_tumor  22G Aug  8 12:37 SpermS02.sort.bam
 102 -rw-r--r-- 1 gaoshengjie bc_tumor  29K Sep  2 06:05 SpermS03.nstat
 103 -rw-r--r-- 1 gaoshengjie bc_tumor  19G Aug 11 15:12 SpermS03.sort.bam