# tileseqscore2.txt # notes for assessing tile array score transforms for gene/exon boundary detection # version 2; see ~/Desktop/dspp-work/daphwork/tilepatharrays/tilednascore.pl #.... redo with new exontab Data sets: Daphnia pulex sc17 (gnomon 400 genes, 2000 exons) and Daphnia sc1 (gnomon, 800 genes, 4000 exons) 47,000 Nimblegen tiles in sc17; 130,000 tiles in sc1 # tilednascore perl reads tile scores, then gene features; writes table # of tile scores at gene/exon boundaries (+/- 75bp now), and GC content cat dpx_sc17_nimsc901.gff dpx_sc17_qnorm.gff \ gdna/cDNA.l2r99sc17.gff gdna/gDNA.l2r99sc17.gff gdna/cDNA-gDNA.l2r99sc17.gff \ dpulex1_sc17_gnomodels.gff | perl ../tilednascore.pl > dpx_sc17_gnosc.gcin4.tab # Pos50 <- (-50:50) # wider look ... # this data has several transform scores and serial position score diffs # gcall<- read.table("dpx_sc17_pasasc.gcin2.tab",header=T) # gcall<- read.table("dpx_sc17_gnosc.gcin2.tab",header=T) gcall<- read.table("dpx_sc17_gnosc.gcin3.tab",header=T) # adds gene start,end 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)) ## fixme: order bound types not by lex sort; add gestart, geend 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 # > gcall[12:16,] # Pos A C G T N Ft sc.cDNA_gDNA sc.cDNA_l2r99 sc.cDNA_qnorm901 sc.cDNA_raw901 # 12 39 1 0 0 0 0 exend -0.195 1.69 0.829 8.72 # 13 38 0 1 0 0 0 exend -0.195 1.69 0.829 8.72 # 14 37 1 0 0 0 0 exend -0.353 1.43 0.816 8.46 # 15 36 0 1 0 0 0 exend -0.353 1.43 0.816 8.46 # 16 35 0 0 0 1 0 exend -0.353 1.43 0.816 8.46 # sc.gDNA_l2r99 df.cDNA_gDNA df.cDNA_l2r99 df.cDNA_qnorm901 df.cDNA_raw901 # 12 1.89 0.000 0.00 0.0000 0.00 # 13 1.89 0.000 0.00 0.0000 0.00 # 14 1.79 0.158 0.26 0.0129 0.26 # 15 1.79 0.000 0.00 0.0000 0.00 # 16 1.79 0.000 0.00 0.0000 0.00 # df.gDNA_l2r99 qb qe qpos tb tipos # 12 0.0000 23472 23367 23328 23305 23 # 13 0.0000 23472 23367 23329 23305 24 # 14 0.0983 23472 23367 23330 23330 0 # 15 0.0000 23472 23367 23331 23330 1 # 16 0.0000 23472 23367 23332 23330 2 #............... # Stat <- "sc.cDNA_raw901" Stats <- names(gcall)[8:12] # 13:17 for score-diffs (df.xxx) Stats <- names(gcall)[8:10] # 3 vars here; sc1 Stat <- Stats[2] # 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 #......... # combine scaffold data sets? sc1 (~800 genes), sc17 (~300 genes) BUT have diff col-score set; drop sc17 extras: names( gcall.gno4.sc1[, c(1:9,11:13,15:20)] ) names( gcall.gno4[, c(1:9,12:14,17:22)] ) gcall2<- rbind( gcall.gno4.sc1[, c(1:9,11:13,15:20)], gcall.gno4[, c(1:9,12:14,17:22)] ) # > dim(gcall2) # [1] 1738468 18 Stats <- names(gcall)[8:10] # 3 vars here; sc1 + sc17 Pos<- sort(unique(gcall$Pos)) #....... # Calculate average scores/position from gene/exon boundaries 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[,4] <- "none" # unfactorizes ?? gcall.gc$Ft <- factor(gcall.gc$Ft,levels=c("gestart","exstart","exend","geend","none"),ordered=T) # option to add sd, ci to add err bars to plot #....... 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_l2r99 != 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 <- 1 + ip - Pos[1] ; iprow <- iprow + length(Pos) * (iFt -1) # bad for empty Ft ** iprow <- iprow + 1 gcall.gc[ iprow, 1:3] <- c(ip, pctgc, mscore) # nums gcall.gc[ iprow, "Ft"] <- Ft # dang factor for (st in Stats) { gcall.gc[ iprow, st] <- mean( gcall[rpt, st]) # NOTE: rpt subset tile transition } # add se after all means for (st in Stats) { sdst <- paste("se.",st,sep="") sem<- sqrt( var( gcall[rpt, st]) / length(gcall[rpt, st]) ) gcall.gc[ iprow, sdst] <- sem } # test per-position cor( gc x score) ? for (st in Stats) { gcst <- paste("rgc.",st,sep="") rgc<- cor( gcall[rpt, st], (gcall[rpt,"G"]>0 | gcall[rpt,"C"]>0) ) # ,use="pairwise.complete.obs", method= pearson,kendall,spearman gcall.gc[ iprow, gcst] <- rgc } } rp <- gcall.gc$Ft == Ft if( sum(rp) > 0) { tscore <- gcall.gc[ rp, "tscore"] gcall.gc[ rp, "pctscore"] <- (tscore - min(tscore))/(max(tscore) - min(tscore)) 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 #.... gcall.gc[,"sc.cDNA_raw901"]<- gcall.gc[,"sc.cDNA_raw901"]-7 # minor scaling for plot #.... 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!) # scaf1 3 stats plStats<- Stats[c(2,1,3)] # cDNA, c-gDNA, gDNA, # scaf1b 4 stats plStats<- Stats[c(2,3,4)] # cDNA, qnorm, gDNA, # names(gcall.gc)[5:8] plFac<- levels(factor(gcall.gc$Ft)) layout(matrix(1:length(plFac),ncol=length(plFac),nrow=1)) par(cex=1) #? layout messes w/ these ylim<- c(0,1) dopct<- F doci<- F docor<- F needleg<- T for (Ft in plFac) { r<- gcall.gc$Ft == Ft if( sum(r) == 0) next ylab <- ifelse(dopct, "Percent", "Tile Score") ylab <- ifelse(docor, "Correlation GC x Score", ylab) if (dopct) ylim <- c(0,1) else { ylim <- range( gcall.gc[r,c("pctGC",plStats)] ); ylim[1] <- min(ylim[1], -1) ylim[2] <- max(ylim[2], 1) } gcoff<- 0 if (docor) { ylim <- c(-0.2,0.2); gcoff <- -0.5 } # fixme #plot( gcall.gc[r,XYgc], plot( gcall.gc[r,"Pos"], gcall.gc[r,"pctGC"]+gcoff, ylim=ylim, type="l", col="black", main=ftlabs[Ft,2], ylab=ylab,xlab="Base pos" ) abline( v=0, col="gray", lty=3) if(!dopct) abline( h=0, col="gray", lty=3) # crosshairs help ipst<- 1 for (st in plStats) { ipst <- ipst+1 ; pst<- ifelse(dopct, paste("p",st,sep=""), st) if(docor) pst<- paste("rgc.",st,sep="") pv <- gcall.gc[r,c("Pos",pst)] points( pv, col=ipst, pch=ipst+18) #? if(doci & !docor) { # not for docor pci<- 1.96 sest<- paste("se.",st,sep="") tmn<- gcall.gc[ r, st] tse<- pci * gcall.gc[ r, sest] tcihi<- (tmn + tse) tcilo<- (tmn - tse) if(dopct) { tmin<- min(tmn) trng<- max(tmn) - tmin tcihi<- (tcihi - tmin)/trng tcilo<- (tcilo - tmin)/trng } lines( gcall.gc[r,"Pos"], tcihi, col=ipst,lty=3) lines( gcall.gc[r,"Pos"], tcilo, col=ipst,lty=3) } # sm<- smooth(gcall.gc[r,pst]); lines(gcall.gc[r,"Pos"],sm, col=ipst, pch=ipst+18, type="l", lty=3) } legloc<- "topright"; legins<- c(-0.01,0); if(Ft == "exend") # for sc. # legloc<- "topleft"; legins<- c(0.0,0); if(Ft == "exend") # for df. { needleg<- F legend(legloc,c("GC",plStats), col=1:(length(plStats)+1), pch=c(NA,19:(18+length(plStats))), inset=legins, 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") } #........................ # # *** Puzzle: daphnia nimgen gDNA shows same gene-bound signals as cDNA, if less precise *** # why? -- should be more random w/ respect to exon/gene start/stop sites # Tests 1? create random tile DNA signal from gDNA ; is it my compute error? x 0. look at absolute signals, not just pctscore: > somewhat different story; gDNA signal is obviously weaker than cDNA however, precision is good (variance low); gDNA can discriminate exon/intron bounds x 2. *** test more scaffolds *** > added scaffold_1 (800 genes, 4000 exons); same effects, more precision x 3. tile score diff along w/ score; global effect DOESNT disappear per-gene serial changes > same puzzle with serial difference in score 4. ? separate groups of ON genes and OFF genes (from cDNA-exon ave scores) : same effect in both? x 5. look at distribution of Pos curves for each gene? not average; need curve from each gene; - is mean of curves any diff from curve of means? - score curve for -40 .. 40 is approx linear for exstart,exend; ave & sd of slope/exon? > this shows some effects for exon bounds (esp. gDNA, qnorm low; cDNA variants high) .. could be score normalizing effect? ~ 6. Is gDNA score tracking GC content? - pick out tiles by GC content variation in same gene part (exon, intron) to test? > some evidence for this, GC x gDNA in INTRON region is highest effect > EXON regions show lowest GC cor for both gDNA, cDNA #..... #6. looking at correlation per-base-position GC x tile score only: - there is small cor ( <= |0.1| ) - the cor for gene start,end is very noisy, not persuasive - the cor at exon/intron bounds is highest in *intron* (lower GC), it is cleaner, and drops to ~0 for exon (higher GC) - gDNA cor in introns is perhaps a bit higher than cDNA cor, e.g gDNA-intron r.mean= 0.08, cDNA-intron r.mean= 0.05 but little diff in gDNA, cDNA in exon parts ( r.gDNA ~ 0.024, r.cDNA ~ 0.015) * test signif diff here (should combine exstart, exend for better power) > Intron region shows strongest GC x gDNA and cDNA correlation. gDNA score has significantly greater correlation than cDNA to GC content. exon area has lower GC correlation for both. > gestart, geend show same mean correlation contrasts as this, with gestart < 0 and geend > 0 both significant, but less than exon start,end all of these cors (0.01 .. 0.08) are significantly > 0 * INTRON region > rp<- gcall.gc$Pos < 0 & gcall.gc$Ft == "exstart" > t.test( gcall.gc[rp,"rgc.sc.cDNA_l2r99"], gcall.gc[rp,"rgc.sc.gDNA_l2r99"]) Welch Two Sample t-test data: gcall.gc[rp, "rgc.sc.cDNA_l2r99"] and gcall.gc[rp, "rgc.sc.gDNA_l2r99"] t = -7.1, df = 147, p-value = 4.84e-11 # exstart t = -6.58, df = 148, p-value = 7.48e-10 # exend mean of x mean of y 0.0523 0.0804 # exstart ; p < 5.e-11 0.0563 0.0832 # exend ; p < 1.e-09 0.0385 0.0616 # gestart ; p < 5.e-04 0.0627 0.0874 # geend ; p < 5.e-05 * EXON region > rp<- gcall.gc$Pos > 0 & gcall.gc$Ft == "exstart" > t.test( gcall.gc[rp,"rgc.sc.cDNA_l2r99"], gcall.gc[rp,"rgc.sc.gDNA_l2r99"]) Welch Two Sample t-test data: gcall.gc[rp, "rgc.sc.cDNA_l2r99"] and gcall.gc[rp, "rgc.sc.gDNA_l2r99"] t = -2.67, df = 145, p-value = 0.008474 # exstart t = -2.20, df = 147, p-value = 0.02965 # exend mean of x mean of y 0.0156 0.0238 # extart ; p < 0.01 0.0220 0.0286 # exend ; p < 0.03 0.0535 0.0463 # gestart ; Not significant 0.0500 0.0641 # geend ; p < 0.02 #0. compare plots of absolute tile signal around gene bounds vs pct signal -- pct/relative signal shows same curve for cDNA, gDNA (esp exon bounds) -- but absolute shows large diff in scores. #2. scaffold_1 looks same as scaffold_17, if less noise (about 2x genes, exons) copied out cDNA,gDNA for dpx scaffold_1; oat:/c7/eugenes/daphnia/genomes/Daphnia_pulex/tilepatharrays/gdna/ # cat dpx_sc17_nimsc901.gff dpx_sc17_qnorm.gff < dont need these now? qnorm of below cdna? # BUT need dna now in _mimsc901.gff .... # scaffold_17 ns_cDNA_raw901 tile 482 539 857.00 . . dna=AAATCACAACTTCTATAAATTTGAAATGACGATTATGATTGATTTACAAA >> no IDs in cdna-tile.gff; need for fasta **; use nimgens original .gff? with probe_id=SCAFFOLD_17FS000002159 # dammit cdsgff2genbank is omitting due to MIN_CDS_LEN; need -a exon or other fix $bg/mb/augustus/scripts/cdsgff2genbank.pl -a exon -t tile -gff cDNA.l2r99sc1.gff -fa dpulex_sc1.fa \ > cDNA.l2r99sc1.pfa perl -pi -e's/\.1// if(/^>/);' cDNA.l2r99sc1a.pfa # dammit .1 from -a exon cat cDNA.l2r99sc1a.pfa cDNA.l2r99sc1a.gff | perl -ne\ 'if(/^#/){ print; $fa=0; } \ elsif(/^scaf/i) { $fa=0; ($id)=m/ID=([^;\s]+)/; $sq=$fa{$id}; s/$/;dna=$sq/ if($sq); print; }\ elsif(/^>(\S+)/){ $id=$1; $fa=1; } elsif($fa and /^[a-z]/i) { chomp; $fa{$id} .= $_; }'\ > cDNA.l2r99sc1b.gff # 139,000 tiles in sc1 vs 47,000 in sc17 cat gdna/cDNA.l2r99sc1.gff gdna/gDNA.l2r99sc1.gff gdna/cDNA-gDNA.l2r99sc1.gff \ dpulex1_sc1_gnomodels.gff | perl ../tilednascore.pl > dpx_sc1_gnosc.gcin4.tab gcall <- read.table("dpx_sc1_gnosc.gcin4.tab",header=T) # score weighted by dist from tile center >> add quantnorm of scaf1.cdna signal to cmpare >> ?? try on gdna signal ?? > dnasc1<- read.table("gdna/cDNA.l2r99sc1.ptab") > dim(dnasc1) [1] 138904 55 > dnasc1[1:4,1:9] V1 V2 V3 V4 V5 V6 V7 V8 V9 1 ns_cDNA_l2r99 5983 6033 -1.89 51 T G C T 2 ns_cDNA_l2r99 26953 27006 -1.97 54 A G A C 3 ns_cDNA_l2r99 27033 27095 -1.94 63 T T T C INDEX <- make.index.structure(dnasc1[,6:55]) dnasc1.qnorm <- quantilenorm( niter<-50, dnasc1[,4], nbases<-50) # signal # cbind(dnasc1[,1:3],qnorm=dnasc1.qnorm, tname=dnasc1[,5], dnasc1[,6:55] )[1:5,1:10] # or write.table.gff: cbind( ref="scaffold_1", src="ns_cDNA_l2r99", type="tile", beg=dnasc1[,2], end=dnasc1[,3], val=dnasc1.qnorm, or=".", ph=".")[1:5,] write.table( cbind( ref="scaffold_1", src="ns_cDNA_qnorm", type="tile", beg=dnasc1[,2], end=dnasc1[,3], val=dnasc1.qnorm, or=".", ph="."), file="dpx_sc1_qnorm.gff",quote=F,row.names=F,col.names=F,sep="\t") #... cat gdna/cDNA.l2r99sc1.gff gdna/gDNA.l2r99sc1.gff gdna/cDNA-gDNA.l2r99sc1.gff \ dpx_sc1_qnorm.gff \ dpulex1_sc1_gnomodels.gff | perl ../tilednascore.pl > dpx_sc1_gnosc.gcin4b.tab #3. done : same puzzle #5. gcall$qb marks each gene slope of xy line: y = a + b * x ; b = (y - a) / x ; scoreslope = (score - score-at-Pos0) / Pos for Pos in -40..40 score[n] = s0 + scoreslope * Pos[n] # do for all Stats of interest ... same loop Stat <- st <- "sc.cDNA_raw901" # all sr.n = 150 for Pos srlab<-"Exon start, tile score slope +-40bp"; r <- gcall$Ft == "exstart" & gcall$Pos > -40 & gcall$Pos < 40; qb <- unique(gcall[r,"qb"]) srlab<- "Gene start, tile score slope +-40bp"; r <- gcall$Ft == "gestart" & gcall$Pos > -40 & gcall$Pos < 40; qb <- unique(gcall[r,"qb"]) gcall.sr <- data.frame() for (iqb in qb) { gc1<- gcall[ r & gcall$qb == iqb, ] if(nrow(gc1) < 3) next sr1<- c(qb=iqb) for (st in Stats) { s0 <- mean( gc1[, st]) svec <- ( gc1[gc1$Pos != 0, st] - s0) / gc1[gc1$Pos != 0,"Pos"] # maybe want cor here: sum( (y-x)(y-x)? / (y2)*(x2) sr1<- cbind(sr1, mn=mean(svec, na.rm=T),sd=sd(svec, na.rm=T)) } gcall.sr <- rbind( gcall.sr, sr1) } names(gcall.sr) <- c("qb",paste( rep(Stats,each=2),names(gcall.sr)[-1],sep=".")) srset<- seq(2,10,2) # mean slope of scores ; sd = seq(3,11,2) i<-2; plot(density(gcall.sr[,i]),col=i,ylim=c(0,50),type="n",main=srlab) for (i in srset) lines(density(gcall.sr[,i]),col=i/2) legend("topright",names(gcall.sr[,srset]), col=srset/2,lwd=1,bty="n") # for gene starts: little difference? > dim(gcall.sr) [1] 346 11 # sc.cDNA_qnorm901.mn is outlier in very narrow, peaked density, r<> 0.00 - 0.02 others extend -0.02 .. 0.06 slope cDNA_l2r99.mn has highest mode/mean gDNA seems lowest slope mean ? #..... # s0 = 6.57 ; mean-s = 6.48 # sr <- mean((gc1[gc1$Pos != 0, st] - s0) / gc1[gc1$Pos != 0,"Pos"]) # sr = 0.00734 # s0 + sr * 40 = 6.86 # gc1[ gc1$Pos == 40, st] = 6.53 # > mean(svec) # [1] 0.00734 # > sd(svec) # [1] 0.0158 # # for (ip in Pos) { cat(ip, s0 + sr * ip, gc1[ gc1$Pos == ip, st], "\n") } # # > dim(gcall[r & gcall$qb == iqb,]) # [1] 151 22 # # > dim(gcall) # [1] 579507 22 # # > r<- gcall$Ft == "gestart"; length(unique(gcall[r,"qb"])) # [1] 353 # # > length(r) # [1] 579507 == all rows in gcall # > sum(r) # [1] 49411 # count is ngenes * Pos (-75..75 = length 151); 151 * 353 = 53303 # # dgbook2% cat dpulex1_sc17_gnomodels.gff | grep -c exon # 2189 # dgbook2% cat dpulex1_sc17_gnomodels.gff | grep -c mRNA # 376 # # > length(gcall$qb) # [1] 579507 # # > length(unique(gcall$qb)) # [1] 2089 << exons for dpx_sc17_gnomon, see above (missing 100?) # # > r<- gcall$Ft == "exstart"; length(unique(gcall[r,"qb"])) # [1] 1707 # inner exon starts # # > r<- gcall$Ft == "gestart"; length(unique(gcall[r,"qb"])) # [1] 353 # gene starts