#!/usr/bin/perl # tilednascore.pl =head1 usage =head2 Version2: use tile gff inputs with several source-scores/tile ($doinbounds) - all tile.gff have gff.type = tile ; not in SO :( ; allow type=region? cat dpx_sc17_nimsc901.gff dpx_sc17_qnorm.gff \ gdna/cDNA.l2r99sc17.gff gdna/gDNA.l2r99sc17.gff gdna/cDNA-gDNA.l2r99sc17.gff \ dpulex1_sc17_pasamodels.gff | perl ../tilednascore.pl \ - NO NEED to sort tile*.gff: sort -k1,1 -k4,4n -k2,2 ? =head2 Version1: use gene model gff (mRNA, exon) cat dpx_sc17_nimsc.ptab | grep 10533901_532 | cat - dpulex1_sc17_pasamodels.gff |\ perl ../tilednascore.pl > dpx_sc17_nimsc901.gcin.tab cat dpx_sc17_qnorm.ptab | grep 10533901_532 | cat - dpulex1_sc17_pasamodels.gff | \ perl ../tilednascore.pl > dpx_sc17_nimsc901qnorm.gcin.tab =head2 Version0: drop use of gff exon,intron sets cat dpx_sc17_nimsc.ptab | grep 10533901_532 | \ cat dpulex1_sc17_pasamodels.{exons,introns} - | perl ../tilednascore.pl =head1 author Don Gilbert, gilbertd@indiana.edu, july 2008 =cut use strict; our $VERSION = 2; our $DEBUG = 1; # ....... options ......... my $FTBINSIZE= 500; # feature bin size (ignore) my $TIBINSIZE= 500; # tile bin size; smaller?, less redundant calls ? my $OFF= 75; # 50; # +/- OFF around boundary for put_exontab; want >tile size (35+bp) # ^note; put_exontab doesn't now check if offset is overlapping other signal (near exon/gene) my $showft= 1; my $doinbounds= 1; # FIXME: TRUE needs inputs: tiletab > gff; FALSE needs gff > tiletab my $doexontab= 1; # fixed default for doinbounds my $skipgenes= 0; # skip gene start/end; only inner exon start/end my $domeansig= 1; # mean value per qpos/boundtype versus all tilescores at that base my $dologscore= 0; # set from tile score max > 10000 my $doweight= 1; # test score weighted by dist from tile center # my $gffpatt="(\t(pasa|NCBI|JGI|AUG))|((ID|Parent)=)"; # fixme? distinguish gff input from ptab input my $gffpatt="(ID|Parent)="; # fixme? distinguish gff input from ptab input my $tilegffpatt="\ttile"; # fixme? distinguish gff input from ptab input # ..... local data ......... if ($VERSION>1) { $doinbounds=1; $doexontab= 1; $domeansig=1 }; # always my (%ftlist, %talist); my %tileset=(); # hash of tile signals/tile-location my %scorekeys=(); # global for many scores/tile my @scorekeys=(); my ($mintile,$maxtile)=(undef,0); my %dologscore = ( # vers2 hack; should input tile.gff w/ transformed score ns_cDNA_raw901 => 1, ); my %sqindex= ( A => "1 0 0 0 0", # ACGTN C => "0 1 0 0 0", # ACGTN G => "0 0 1 0 0", # ACGTN T => "0 0 0 1 0", # ACGTN N => "0 0 0 0 1", # ACGTN ); my $didheader= 0; # ..... input/output loop ........ if($doinbounds) { # put_exontabhead(); # put_featabhead(); while(<>) { next if(/^#/); chomp; my @v= split"\t"; if(/$tilegffpatt/) { get_tilegff(@v); } else { put_exontab(@v); } # vers2; assumes mRNA>exon sorting #Vers1: if(/$gffpatt/) { put_exontab(@v); } else { get_tiletab(@v); } # assumes mRNA>exon sorting #Vers0: if(/$gffpatt/) { put_featab(@v); } else { get_tiletab(@v); } } } else { # obsolete option? warn "OBSOLETE OPTION doinbound: $doinbounds\n"; while(<>) { next if(/^#/); chomp; my @v=split"\t"; if(/$gffpatt/) { get_featab(@v); } else { put_tiletab(@v); } } } info_exontab(); # debug info #......... end ............... =head2 put_exontab reads GFF exon,mRNA/gene lines only (no intron) writes table of positions from boundarys (intron-start, intron-end, gene-start, gene-end) with Pos ; A,C,G,T,N ; Boundtype: ib,ie,gb,ge; tilescore, ... =cut sub put_exontabhead { if($VERSION > 1) { my @scorecols= @_; print join("\t","Pos","A C G T N","Ft", @scorecols, # add several tscore, tdiff here?? need to know how many before print # drop qt, te as useless "qb qe","qpos","tb","tipos", # accessory fields; not needed? ),"\n"; } else { # print join("\t",$Pos,$Seqi,$Typ,$tv,$tdiff,$qt,"$qb $qe",$qpos,"$tb $te",$tipos),"\n"; print join("\t","Pos","A C G T N","Ft","tscore","tdiff","qt","qb qe","qpos","tb te","tipos"),"\n"; # header } } my($ntile, $ntiledup, $ngene, $nexon); my($gebegin,$geend,$gid); my($lxbegin,$lxend,$lxid); sub info_exontab { my $grtile= ($ntile<1) ? 0 : sprintf "%.2f", $ntiledup / $ntile; warn "# Kept ngene=$ngene, nexon=$nexon, ntile=$ntile, grp/tile=$grtile\n"; } sub put_exontab { #my @v= @_; my($qr,$qt,$qb,$qe,$qo,$qd)= @v[0,2,3,4,6,8]; my($qr,$qs,$qt,$qb,$qe,$qp,$qo,$qf,$qd)= @_; my $isgene= ($qt =~ /^mRNA/); #?? |gene my $isexon= ($qt =~ /^exon/); return unless($isexon || $isgene); $dologscore=1 if($maxtile>10000 && $mintile>0); #? here each call? unless($didheader++) { if($VERSION>1) { @scorekeys= sort keys %scorekeys; my @scorecols= (); # sort keys %scorekeys; push @scorecols, map { "score.$_"; } @scorekeys; push @scorecols, map { "diff.$_"; } @scorekeys; put_exontabhead(@scorecols); } else { put_exontabhead(); } } my $qlen1= ($qe - $qb); # length - 1 my ($id)= $qd =~ m/(?:ID|Parent)=([^;\s]+)/; warn "ERROR: no ID: @_ \n" unless($id); my $keep= 0; my %didta=(); my ($Typ,$Pos,$Seqi); # save gene bounds for exon1st .. exonInt .. exonEnd my $rev= ($qo eq "-")? 1: 0; if($rev) { ($qb,$qe)= ($qe,$qb); } my $whichexon=""; # single,first,inner,last if($isexon) { warn "ERROR: gene sort order: exon $id ne gene $gid\n" if($id ne $gid); $whichexon= ($qb == $gebegin && $qe == $geend) ? "single" : ($qb == $gebegin) ? "first" : ($qe == $geend) ? "last" : "inner"; $nexon++; #at end# ($lxbegin,$lxend,$lxid)= ($qb,$qe,$id); } elsif($isgene) { ($gebegin,$geend,$gid)= ($qb,$qe,$id); $ngene++; return if ($skipgenes); } ## return if ($DEBUG>1); # for counts only... # dang; these depend on $qo orientation strand my @pos1= (($qb - $OFF) .. ($qb + $OFF)); my @posE= (($qe - $OFF) .. ($qe + $OFF)); #?? if($rev) { @pos1= reverse @pos1; } #?? foreach my $part (qw(start end)) { next if($isexon && $part eq "start" && $whichexon =~ /first|single/); next if($isexon && $part eq "end" && $whichexon =~ /last|single/); my @pos = ($part eq "start") ? @pos1 : @posE; my $q0 = ($part eq "start") ? $qb : $qe; my $Typ = ($isgene) ? "ge" : "ex"; $Typ .= $part; %didta=(); my $lmeanv= undef; my %lmeanv=(); #v2 my $ipos= 0; foreach my $qpos (@pos) { $Pos= $qpos - $q0; $ipos++; $Pos= -$Pos if ($rev); # *!* my $ib= int($qpos / $TIBINSIZE); $talist{$ib} or next; # add option: mean of all @ft signals per qpos/Pos/Type (different tiles) my ($sumv,$nv,$tipos1,$tb1,$te1)= (0) x 10; my (%sumv); # ver2 my @ft= @{ $talist{$ib} }; foreach my $ft (@ft) { my($tb,$te,$tv,$tr,$typ,$tid,$sq,@scores)= @{$ft}; # VERSION 2 adds @scores my $over= ($qpos <= $te && $qpos >= $tb) ? 1 : 0; next unless($over); my $ld; if($VERSION > 1) { $ld= "$qpos.$tid"; # uniq id : must include qpos ! } else { $ld= "$qpos.$tr.$tb.$te.$tid"; # uniq id : must include qpos ! @scores=("score.$typ=$tv"); } next if($didta{$ld}); if($over) { # already checked $keep++; $didta{$ld}++; my @sq= split " ",uc($sq); my $nq= $#sq; # slide tile seq over qpos my $tipos = $qpos - $tb; # offset of tile from signal +/- steps (-20..20) next if($tipos > $nq || $tipos < 0); $Seqi= $sqindex{ $sq[$tipos] } || "0 0 0 0 0"; # version2 adds many scores/tile if($VERSION < 2) { $tv= log($tv) if($dologscore); my $tdiff= 0; # not for no-domeansig print join("\t",$Pos,$Seqi,$Typ,$tv,$tdiff,$qt,"$qb $qe",$qpos,"$tb $te",$tipos),"\n" unless($domeansig); } ($tb1,$te1,$tipos1)= ($tb,$te,$tipos); # should we add a score weight based on distance of qpos from tile center? just use tipos? my $tiwt= 1; if($doweight) { my $tiw = abs($te - $tb)/2; $tiwt = ($tiw - abs( $qpos - ($tb + $te)/2 ))/ $tiw; # how far from center of tile; overlapped tiles? $tv= $tv * $tiwt; } $sumv += $tv; $nv += $tiwt; # $nv++; if($VERSION>1) { foreach my $sc (@scores) { my ($sk,$sv)= split"=",$sc; $sk =~ s/^score\.//; if($doweight) { $sv= $sv * $tiwt; } $sumv{$sk} += $sv; } } } } # tipos, tb, te change over mean; qpos, qb, qe, rest are same my($meanv, $tdiff)= (0,0); if($nv) { # version2 adds many scores/tile if($VERSION>1) { my $nolast= (not defined $lmeanv); my (@meanv,@tdiff); foreach my $sk (@scorekeys) { # sort keys %sumv my $sv= $sumv{$sk}; $lmeanv= $lmeanv{$sk}; $meanv= sprintf "%.5g", $sv / $nv; $tdiff= sprintf "%.5g", (!defined $lmeanv) ? 0 : ($rev) ? ($lmeanv - $meanv) : ($meanv - $lmeanv); push @meanv, $meanv; push @tdiff, $tdiff; $lmeanv{$sk}= $meanv; } print join("\t",$Pos,$Seqi,$Typ,@meanv,@tdiff,"$qb $qe",$qpos,$tb1,$tipos1),"\n" if ($domeansig); } else { $meanv= sprintf "%.5g", $sumv / $nv; $tdiff= sprintf "%.5g", (! defined $lmeanv) ? 0 : ($rev) ? ($lmeanv - $meanv) # this may be wrong; : ($meanv - $lmeanv); $lmeanv= $meanv; print join("\t",$Pos,$Seqi,$Typ,$meanv,$tdiff,$qt,"$qb $qe",$qpos,"$tb1 $te1",$tipos1),"\n" if ($domeansig); } } } } if($isexon) { #? both exon, gene? ($lxbegin,$lxend,$lxid)= ($qb,$qe,$id); } } sub get_tiletab { # VERSION 2 ; use get_tilegff instead my($tr,$tb,$te,$tv,$tid,$sq,@moresq)= @_; $maxtile= $tv if ($tv>$maxtile); $mintile= $tv unless (defined $mintile && $tv > $mintile); if(@moresq) { $sq= join(" ",$sq,@moresq); } my $ft= [$tb,$te,$tv,$tr,"tile",$tid,$sq]; # $qb,$qe,$qt,$qd my @bins= (int($tb/$TIBINSIZE) .. int($te/$TIBINSIZE)); foreach my $ib (@bins) { push( @{$talist{$ib}}, $ft); } } # new version2; see below ; have separate rows at same tile for diff source/score transforms sub get_tilegff { my($tr,$ts,$tt,$tb,$te,$tv,$or,$ph,$at)= @_; # skip 6,7 strand,phase my $tid= "$tr.$tb.$te"; # location only for tile id? my $sq = ""; if ($at =~ m/dna=(\w+)/) { # only in 1 of several tile gff files $sq= $1; $sq= join " ", split "",$sq; # insert spaces for old version format } # my %at= map{ my($k,$v)=split"=",$_; $k=>$v; } split";",$at; # assume clean # my $sq= $at{'dna'}; # fixme? if( $dologscore{$ts} && $tv > 0) { $tv= log($tv); } # transform this source's raw score $scorekeys{$ts}++; my $ft = $tileset{$tid}; if( $ft ) { push( @$ft, "score.$ts=$tv"); $ft->[6]= $sq if($sq); $ntiledup++; } else { # $ft= [$tb,$te,$tv,$tr,"tile",$tid,$sq]; # orig; $ft= [$tb,$te,$tv,$tr,"tile",$tid,$sq,"score.$ts=$tv"]; $tileset{$tid}= $ft; $ntile++; $ntiledup++; my @bins= (int($tb/$TIBINSIZE) .. int($te/$TIBINSIZE)); foreach my $ib (@bins) { push( @{$talist{$ib}}, $ft); } } } =head1 OBSOLETE METHODS - get_featab - put_featab - put_tiletab =cut sub put_tiletab { # OBSOLETE my($tr,$tb,$te,$tv,$tid,$sq)= @_; my $keep=0; my %didft=(); my @sq= map {$_=" "} split " ",$sq; my $nq= $#sq; my @bins= (int($tb/$FTBINSIZE) .. int($te/$FTBINSIZE)); foreach my $ib (@bins) { $ftlist{$ib} or next; my @ft= @{ $ftlist{$ib} }; foreach my $ft (@ft) { my($lb,$le,$lt,$ld)= @{$ft}[0,1,2,3]; next if($didft{$ld}); my $over= ($tb <= $le && $te >= $lb) ? 1 : 0; if($over) { $keep++; $didft{$ld}++; my $tp= substr($lt,0,1); # for my $ti ($tb .. $te) for my $t0 (0..$nq) { my $ti = $t0 + $tb; # my $t0= $ti - $tb; if ($ti >= $lb && $ti <= $le) { my $sv= $tp; if ($sq[$t0] =~ /e|i|b/ && $sq[$t0] ne $sv) { $sq[$t0] = "b"; } else { $sq[$t0]= $sv; } } } } } } map{ $_="o" if(/ /) }@sq; $sq = join(" ",@sq) if($showft); print join("\t",$tr,$tb,$te,int($tv),$sq),"\n" if $keep; # rewrite dpx_sc17_*.ptab with e,i,o key for seq } sub get_featab { # OBSOLETE my @v= @_; my($qr,$qt,$qb,$qe,$qo,$qd)= @v[0,2,3,4,6,8]; return unless($qt =~ /exon|intron|CDS/); # filter by $qt here? $qd =~ s/Parent=//; my $ft= [$qb,$qe,$qt,$qd]; # no $qo my @bins= (int($qb/$FTBINSIZE) .. int($qe/$FTBINSIZE)); foreach my $ib (@bins) { push( @{$ftlist{$ib}}, $ft); } } =head2 put_featab write exon/intron boundary values; offsets: inbegin: exone/-20 .. intronb/+20; inend: introne/-20 .. exonb/+20 Ft: ib = intron begin; ie = intron end Pos A C G T N Ft tscore < scores or bool 1/0 ? -19 1 0 0 0 0 ib 4000 .. 0 .. 20 * problem, for put1 need only models w/ both exon,intron cat dpulex1_sc17_pasamodels.{exons,introns} | sort -k1,1 -k9,9 -k4,4n |\ perl -ne'@v=split"\t"; ($d)=m/Parent=Model.(\w+)/; $p=0 if($d ne $ld); \ if(/intron/){ $p=1; print $lv; } print if($p); $lv=$_; $ld=$d;' \ > dpulex1_sc17_pasamodels.exintrons ** change to use only exons, mRNA for gene start/stop =cut sub put_featabhead { # print join("\t","Pos","A C G T N","Ft","tscore"),"\n"; # header # print join("\t",$Pos,$Seqi,$Ft,$tv,$qt,"$qb $qe","$tb $te",$tin),"\n"; print join("\t","Pos","A C G T N","Ft","tscore","qt","qb qe","tb te","tin"),"\n"; # header } sub put_featab { # OBSOLETE my @v= @_; my($qr,$qt,$qb,$qe,$qo,$qd)= @v[0,2,3,4,6,8]; my $qlen1= ($qe - $qb); # length - 1 my $keep= 0; my %didta=(); my ($Ft,$Pos,$Seqi); my $isintron= ($qt =~ /^i/); my $isexon = ($qt =~ /^e/); my @bins= (int($qb/$TIBINSIZE) .. int($qe/$TIBINSIZE)); foreach my $ib (@bins) { $talist{$ib} or next; my @ft= @{ $talist{$ib} }; foreach my $ft (@ft) { my($tb,$te,$tv,$tr,$typ,$tid,$sq)= @{$ft}; my $ld= "$tr.$tb.$te.$tid"; #? uniq id next if($didta{$ld}); my $over= ($qb <= $te && $qe >= $tb) ? 1 : 0; if($over) { $keep++; $didta{$ld}++; my @sq= split " ",uc($sq); my $nq= $#sq; # always 50 now ? my $tin= $tb - $qb; # inset of tile into feature, may be - before ft if($isintron) { $Ft="ib"; foreach $Pos ( 1 .. 20 ) { # 1 == first intron base my $qi= $tb + $Pos - 1; # index to qfeat next if($qi > $qe || $qi < $qb); my $tipos= $tin + $Pos - 1; # 1..20 if tile starts at ft next if ($tipos < 0 || $tipos > 19); my $ti= $tipos;# - 1; # $Pos - 1; # index to @sq next if($ti > $nq); $Seqi= $sqindex{ $sq[$ti] } || "0 0 0 0 0"; print join("\t",$Pos,$Seqi,$Ft,$tv,$qt,"$qb $qe","$tb $te",$tin),"\n"; } $Ft="ie"; foreach $Pos ( -19 .. 0 ) { # 0 == last intron base my $qi= $te + $Pos; # index to qfeat next if($qi < $qb || $qi > $qe); # my $tipos= $qlen1 + ($tin + $Pos); # -19 .. 0 if tile starts at ft # next if ($tipos < 0 || $tipos > $qlen1); #??? my $tipos= $qe - $tb + $Pos; # tiend = tbegin + tlen-1 = $qe - $tb; my $ti= $tipos; # - 1; # index to @sq # my $ti= $qlen1 + $Pos; # index to @sq; may be >> @seq next if($ti > $nq); $Seqi= $sqindex{ $sq[$ti] } || "0 0 0 0 0"; print join("\t",$Pos,$Seqi,$Ft,$tv,$qt,"$qb $qe","$tb $te",$tin),"\n"; } } elsif($isexon) { # revise to use exons only; mark Ft : ib, ie, gb, ge (gene) # need @bins and tile scores from $qb - 20 .. $qe + 20 $Ft="ib"; #?? exon-ib not seen? no see it; see double exon-ie ? foreach $Pos ( -19 .. 0 ) { my $qi= $te + $Pos; # index to qfeat next if($qi < $qb || $qi > $qe); # my $tipos= $qlen1 + ($tin + $Pos); # -19 .. 0 if tile starts at ft # next if ($tipos < 0 || $tipos > $qlen1); #??? my $tipos= $qe - $tb + $Pos; # tiend = tbegin + tlen-1 = $qe - $tb; my $ti= $tipos; # - 1; # index to @sq # my $ti= $qlen1 + $Pos; # index to @sq; may be >> @seq next if($ti > $nq); $Seqi= $sqindex{ $sq[$ti] } || "0 0 0 0 0"; print join("\t",$Pos,$Seqi,$Ft,$tv,$qt,"$qb $qe","$tb $te",$tin),"\n"; } $Ft="ie"; foreach $Pos ( 1 .. 20 ) { # 1 == first exon base my $qi= $tb + $Pos - 1; # index to qfeat next if($qi > $qe || $qi < $qb); my $tipos= $tin + $Pos - 1; # 1..20 if tile starts at ft next if ($tipos < 0 || $tipos > 19); my $ti= $tipos;# - 1; # $Pos - 1; # index to @sq # my $ti= $Pos; # index to @sq next if($ti > $nq); $Seqi= $sqindex{ $sq[$ti] } || "0 0 0 0 0"; print join("\t",$Pos,$Seqi,$Ft,$tv,$qt,"$qb $qe","$tb $te",$tin),"\n"; } } } } } } =head1 tiletab data OLD: file variants of daphnia nimblegen tile location, signal, bases a. raw signal / 1 track; b. GN-quantnorm NEW: combine below a, b, c, d, e in one file (same loc, seq, diff signal transforms) do as GFF format? where to put all signals? gff.score a,b,c,d? or gff.9 a=1;b=1;c=1;... where to put seq: gff.9 seq=AGCTTTTT... dgbook2% wc -l *.ptab 286956 dpx_sc17_nimsc.ptab = all tracks 47826 dpx_sc17_nimsc901.ptab = a. raw signal sample track1 47826 dpx_sc17_qnorm.ptab = b. Gerstein quantilenorm() signal dgbook2% wc -l *DNA*sc17.gff : as gff format 47826 cDNA.l2r99sc17.gff = c. modest log, normal transform of all tracks raw signal 47826 gDNA.l2r99sc17.gff = d. gDNA signal log transform 47826 cDNA-gDNA.l2r99sc17.gff = e. cDNA - gDNA (after log transforms) ==> dpx_sc17_nimsc901.ptab <== source: ns_cDNA_raw901 scaffold_17 482 539 857.00 10533901_532 A A A T C A C A A C T T C T A T A A A T T T G A A A T G A C G A T T A T G A T T G A T T T A C A A A scaffold_17 1834 1894 501.44 10533901_532 C T T C A A A A A A T T T G A A A T G A C G A T T C T G A T A G A T A G A A A T A G A T T C A T A G scaffold_17 1854 1917 745.67 10533901_532 C G A T T C T G A T A G A T A G A A A T A G A T T C A T A G T A G C T G G G T C A T T C C A T T T T cat dpx_sc17_nimsc901.ptab | perl -ne\ 'chomp; ($r,$b,$e,$v,$src,$sq)=split"\t"; $sq=~s/\s+//g; \ print join("\t",$r,"ns_cDNA_raw901","tile",$b,$e,$v,".",".","dna=$sq"),"\n";' \ > dpx_sc17_nimsc901.gff ==> dpx_sc17_qnorm.ptab <== source: ns_cDNA_qnorm901 scaffold_17 482 539 0.683013269576169 10533901_532 A A A T C A scaffold_17 1834 1894 0.0484065494515341 10533901_532 C T T C A G scaffold_17 1854 1917 0.454126100256601 10533901_532 C G A T T T cat dpx_sc17_qnorm.ptab | perl -ne\ 'chomp; ($r,$b,$e,$v,$src,$sq)=split"\t"; $v=sprintf "%.5g",$v; \ print join("\t",$r,"ns_cDNA_qnorm901","tile",$b,$e,$v,".","."),"\n";' \ > dpx_sc17_qnorm.gff ==> cDNA-gDNA.l2r99sc17.gff <== scaffold_17 ns_cDNA_gDNA region 482 539 -3.8 . . scaffold_17 ns_cDNA_gDNA region 1834 1894 -3.05 . . scaffold_17 ns_cDNA_gDNA region 1854 1917 -2.65 . . ==> cDNA.l2r99sc17.gff <== scaffold_17 ns_cDNA_l2r99 region 482 539 -0.66 . . scaffold_17 ns_cDNA_l2r99 region 1834 1894 -1.96 . . scaffold_17 ns_cDNA_l2r99 region 1854 1917 -1.06 . . ==> gDNA.l2r99sc17.gff <== scaffold_17 ns_gDNA_l2r99 region 482 539 3.14 . . scaffold_17 ns_gDNA_l2r99 region 1834 1894 1.09 . . scaffold_17 ns_gDNA_l2r99 region 1854 1917 1.59 . . ...... cat dpx_sc17_nimsc901.gff dpx_sc17_qnorm.gff \ gdna/cDNA.l2r99sc17.gff gdna/gDNA.l2r99sc17.gff gdna/cDNA-gDNA.l2r99sc17.gff |\ sort -k1,1 -k4,4n -k2,2 \ scaffold_17 ns_cDNA_gDNA region 482 539 -3.8 . . scaffold_17 ns_cDNA_l2r99 region 482 539 -0.66 . . scaffold_17 ns_cDNA_qnorm901 tile 482 539 0.68301 . . scaffold_17 ns_cDNA_raw901 tile 482 539 857.00 . . dna=AAATCACAACTTCTATAAATTTGAAATGA CGATTATGATTGATTTACAAA scaffold_17 ns_gDNA_l2r99 region 482 539 3.14 . . scaffold_17 ns_cDNA_gDNA region 1834 1894 -3.05 . . scaffold_17 ns_cDNA_l2r99 region 1834 1894 -1.96 . . scaffold_17 ns_cDNA_qnorm901 tile 1834 1894 0.048407 . . scaffold_17 ns_cDNA_raw901 tile 1834 1894 501.44 . . dna=CTTCAAAAAATTTGAAATGACGATTCTGA TAGATAGAAATAGATTCATAG =head2 vers2 Rplots gcall<- read.table("dpx_sc17_pasasc.gcin3.tab",header=T) # adds gene start,end gcall<- read.table("dpx_sc17_gnosc.gcin4.tab",header=T) # score weighted by dist from tile center Pos<- sort(unique(gcall$Pos)) gcall$Ft<- factor(gcall$Ft,levels=c("gestart","exstart","exend","geend","none"),ordered=T) hd<- names(gcall) hd1<-gsub("score.","sc.",gsub("diff.","df.",gsub("ns_","",hd))) names(gcall)<- hd1 Stat <- "sc.cDNA_raw901" Stats <- names(gcall)[8:12] # stats: grep("(sc|df).",hd1,value=T) # [1] "sc.cDNA_gDNA" "sc.cDNA_l2r99" "sc.cDNA_qnorm901" "sc.cDNA_raw901" # [5] "sc.gDNA_l2r99" "df.cDNA_gDNA" "df.cDNA_l2r99" "df.cDNA_qnorm901" # [9] "df.cDNA_raw901" "df.gDNA_l2r99" FtLevels <- levels(factor(gcall$Ft)) # skip unused gcall.gc <- data.frame( matrix( data=c(0,0,0,0), ncol=4, nrow=( length(FtLevels) * length(Pos)), byrow=T)) colnames(gcall.gc) <- c("Pos","pctGC","tscore","Ft") # change tscore, pctscore to Stat,"%".Stat ? gcall.gc$Ft <- "none" # unfactorizes ?? gcall.gc$Ft <- factor(gcall.gc$Ft,levels=c("gestart","exstart","exend","geend","none"),ordered=T) iprow <- 0 for (Ft in FtLevels) { iFt <- grep(Ft, FtLevels) for (ip in Pos) { rp <- gcall[,"Pos"] == ip & gcall[,"Ft"] == Ft rpt <- rp & gcall$df.cDNA_raw901 != 0 # drop missing transitions? if(sum(rp) == 0) next # no empty Ft pctgc <- sum( gcall[rp,"G"]>0 | gcall[rp,"C"]>0) / sum(rp) mscore <- mean( gcall[rp, Stat] ) iprow <- iprow + 1 gcall.gc[ iprow, 1:3] <- c(ip, pctgc, mscore) # nums gcall.gc[ iprow, "Ft"] <- Ft for (st in Stats) { gcall.gc[ iprow, st] <- mean( gcall[rpt, st]) # NOTE: rpt subset tile transition } } rp <- gcall.gc$Ft == Ft if( sum(rp) > 0) { tscore <- gcall.gc[ rp, "tscore"] pctscore <- (tscore - min(tscore))/(max(tscore) - min(tscore)) gcall.gc[ rp, "pctscore"] <- pctscore for (st in Stats) { pst<- paste("p",st,sep="") tscore<- gcall.gc[ rp, st] gcall.gc[ rp, pst] <- (tscore - min(tscore))/(max(tscore) - min(tscore)) } } } # plot tilescore/gc at gene/exon boundaries ftlabs <- matrix( c( "gestart", "Gene Start Boundary", "exstart","Intron -> Exon Boundary", "exend", "Exon -> Intron Boundary", "geend", "Gene End Boundary", "none", "no Boundary" ), ncol=2, byrow=T); rownames(ftlabs)<- ftlabs[,1] XYgc<- c("Pos","pctGC") XYscore<- c("Pos","pctscore") # more choices here plStats<- Stats[c(2,5)] # cDNA, gDNA plStats<- Stats[c(2,1,5)] # cDNA, c-gDNA, gDNA # plStats<- Stats[c(2,1)] # cDNA, cDNA-gDNA plStats<- Stats[c(4,3)] # raw,qnorm plStats<- Stats[c(3,5,2)] # cDNA, gDNA, qnorm(of raw!) plFac<- levels(factor(gcall.gc$Ft)) layout(matrix(1:length(plFac),ncol=length(plFac),nrow=1)) par(cex=1) #? layout messes w/ these needleg<- T for (Ft in plFac) { r<- gcall.gc$Ft == Ft if( sum(r) == 0) next plot( gcall.gc[r,XYgc], ylim=0:1, type="l", col="black", main=ftlabs[Ft,2], ylab="Percent",xlab="Base pos" ) ipst<- 1 for (st in plStats) { ipst <- ipst+1 ; pst<- paste("p",st,sep="") points(gcall.gc[r,c("Pos",pst)], ylim=0:1, col=ipst, pch=ipst+18) #? sm<- smooth(gcall.gc[r,pst]) # , kind="3R" lines(gcall.gc[r,"Pos"],sm, ylim=0:1, col=ipst, pch=ipst+18, type="l", lty=3) } if(Ft == "geend") { needleg<- F legend("topright",c("GC",plStats), col=1:(length(plStats)+1), pch=c(NA,19:(18+length(plStats))), inset=c(-0.01,0), lwd=c(1,rep(NA,length(plStats))), bty="n") } } if(needleg) { didleg<- F legend("topright",c("GC",plStats), col=1:(length(plStats)+1), pch=c(NA,19:(18+length(plStats))), inset=c(-0.01,0), lwd=c(1,rep(NA,length(plStats))), bty="n") } =head2 vers1 Rplots > gcind<- read.table("dpx_sc17_nimsc901.gcin.tab",header=T) > dim(gcind) [1] 33764 14 > gcind[1:10,] Pos A C G T N Ft tscore qt qb qe tb te tin 1 11 0 0 0 1 0 ie 2500 exon 50845 50958 50835 50893 -10 2 12 1 0 0 0 0 ie 2500 exon 50845 50958 50835 50893 -10 3 13 0 0 0 1 0 ie 2500 exon 50845 50958 50835 50893 -10 4 14 1 0 0 0 0 ie 2500 exon 50845 50958 50835 50893 -10 > rp<- gcind[,"Pos"] == -19 & gcind[,"Ft"] == "ib" > sum(gcind[rp,"G"]>0 | gcind[rp,"C"]>0)/sum(gcind[rp,"tscore"]>0) > mean(log(gcind[rp,"tscore"])) [1] 6.8 Pos <- (-19:20) gcind.gc <- data.frame( matrix( data=c(0,0.1,0,0), nrow=length(Pos), ncol=4, byrow=T)) # Pos, pctGC, mtscore colnames(gcind.gc) <- c("Pos","pctGC","tscore","type") rownames(gcind.gc) <- Pos for (ip in Pos) { tp <- "ie" # ib rp<- gcind[,"Pos"] == ip & gcind[,"Ft"] == tp pctgc<- sum(gcind[rp,"G"]>0 | gcind[rp,"C"]>0)/sum(gcind[rp,"tscore"]>0) mtscore<- mean(log(gcind[rp,"tscore"])) type<- ifelse(tp == "ib", 1, 2) gcind.gc[ 1 + ip - Pos[1],] <- c(ip, pctgc, mtscore, type) } tscore <- gcind.ib.gc[,"tscore"] gcind.ib.gc[,"pctscore"] <- pctscore <- (tscore - min(tscore))/(max(tscore) - min(tscore)) tscore <- gcind.ie.gc[,"tscore"] gcind.ie.gc[,"pctscore"] <- pctscore <- (tscore - min(tscore))/(max(tscore) - min(tscore)) layout(matrix(c(1,2),ncol=2,nrow=1)) plot(gcind.ie.gc[,1:2],ylim=0:1,main="Intron->Exon Boundary",type="l") lines(ine.gc,col="gray",lty=2) # true... lines(gcind.ie.gc[,c(1,5)],ylim=0:1,col="red") plot(gcind.ib.gc[,1:2],ylim=0:1,main="Exon->Intron Boundary",type="l") lines(ins.gc,col="gray",lty=2) # true... lines(gcind.ib.gc[,c(1,5)],ylim=0:1,col="red") legend("topright",c("SampleGC%","FullGC%","TileScore%"),col=c("black","gray","red"),lwd=1) =cut __END__