### R code from vignette source 'TransView.Rnw' ################################################### ### code chunk number 1: options ################################################### options(width=65) ################################################### ### code chunk number 2: install (eval = FALSE) ################################################### ## if (!requireNamespace("BiocManager", quietly=TRUE)) ## install.packages("BiocManager") ## BiocManager::install("TransView") ################################################### ### code chunk number 3: load_files ################################################### library("TransView") library("GenomicRanges") library("pasillaBamSubset") fn.chipseq.bam<-dir(system.file("extdata", package="TransView"),full=T,patt="bam$") fn.macs<-dir(system.file("extdata", package="TransView"),full=T,patt="xls$") fn.dm3.gtf<-dir(system.file("extdata", package="TransView"),full=T,patt="gtf.gz$")[1] fn.mm9.gtf<-dir(system.file("extdata", package="TransView"),full=T,patt="gtf.gz$")[2] fn.pas_paired<-untreated1_chr4() fn.pas_upaired<-untreated3_chr4() ################################################### ### code chunk number 4: generate_map_1 ################################################### dens.ind<-parseReads(fn.chipseq.bam[2],verbose=0,description="Induced") dens.wt<-parseReads(fn.chipseq.bam[1],verbose=0,description="Basal") ################################################### ### code chunk number 5: bamshow_1 ################################################### dens.ind ################################################### ### code chunk number 6: generate_map_1 ################################################### tvs<-tvStats(dens.ind) tvs$nreads ################################################### ### code chunk number 7: generate_hist ################################################### dens.ind.hist<-histogram(dens.ind) dens.wt.hist<-histogram(dens.wt) ################################################### ### code chunk number 8: barplot1 ################################################### barplot(rbind(dens.ind.hist[1:50]+1,dens.wt.hist[1:50]+1),xlab="Read count",ylab="Positions",col=c("blue","red"),beside=T,legend.text=c("Induced","Basal"),log="y",args.legend=list(x = "topright")) ################################################### ### code chunk number 9: generate_map_2 ################################################### peaks<-macs2gr(fn.macs,psize=500) dens.ind.filt<-parseReads(fn.chipseq.bam[2],verbose=0,description="ChIP",set_filter=peaks[1,]) us<-slice1(dens.ind,chrom=as.character(seqnames(peaks[1])),start=start(peaks[1]),end=end(peaks[1])) rs<-slice1(dens.ind.filt,chrom=as.character(seqnames(peaks[1])),start=start(peaks[1]),end=end(peaks[1])) all(us==rs) size(dens.ind.filt) size(dens.ind) ################################################### ### code chunk number 10: sliceN ################################################### slices.ind<-sliceN(dens.ind,ranges=peaks) slices.wt<-sliceN(dens.wt,ranges=peaks) ################################################### ### code chunk number 11: slice_2 (eval = FALSE) ################################################### ## slices.nobckgd<-sliceN(dens.ind,control=dens.wt,ranges=peaks,treads_norm=F) ## plot(slices.ind[[1]],ylab="Total reads") ## lines(slices.wt[[1]],type="p",col=4) ## lines(slices.nobckgd[[1]],type="p",col=2) ## legend(400,150,c("Induced","Basal","Induced corrected"),col=c(1,4,2),pch="o",bty="n") ## slices.nobckgd.fc<-sliceN(dens.ind,ranges=peaks,control=dens.wt,input_method="/",treads_norm=F) ## summary(slices.nobckgd.fc[[1]]) ################################################### ### code chunk number 12: generate_map_4 ################################################### dens.pas_upaired<-parseReads(fn.pas_upaired,spliced=T,verbose=0,description="Unpaired") dens.pas_paired<-parseReads(fn.pas_paired,spliced=T,verbose=0,description="Paired") ################################################### ### code chunk number 13: GTF ################################################### gtf.mm9<-gtf2gr(fn.mm9.gtf) gtf.dm3<-gtf2gr(fn.dm3.gtf) head(gtf.mm9) ################################################### ### code chunk number 14: spliced_1 ################################################### dens.pas_paired.filt<-parseReads(fn.pas_paired,spliced=T,verbose=0,description="RNA-Seq",set_filter=gtf.dm3) size(dens.pas_paired)/size(dens.pas_paired.filt) ################################################### ### code chunk number 15: spliced_2 ################################################### slices.exprs.pangolin<-slice1T(dens.pas_paired.filt,"NM_001014685",gtf.dm3,stranded=T,concatenate=F) pangolin.exon.12<-slices.exprs.pangolin[["NM_001014685.12"]] pangolin.exon.12 ################################################### ### code chunk number 16: spliced_3 ################################################### #Just for demonstration all refseq ids are taken, not recommended for a full sized GTF! all_ids<-unique(mcols(gtf.dm3)$transcript_id) slices.exprs<-sliceNT(dens.pas_paired,all_ids,gtf.dm3,stranded=T,concatenate=F) pangolin.all.exon.12<-slices.exprs[["NM_001014685.12"]] all(pangolin.exon.12,pangolin.all.exon.12) ################################################### ### code chunk number 17: spliced_4 ################################################### peaks.anno<-annotatePeaks(peaks=peaks,gtf=gtf.mm9,limit=2e3) ################################################### ### code chunk number 18: spliced_5 ################################################### peaks.tss<-peak2tss(peaks.anno, gtf.mm9,peak_len=1000) ################################################### ### code chunk number 19: pTV1 ################################################### cluster_results<-plotTV(dens.ind,dens.wt,regions=peaks.anno,show_names=T,norm_readc=F,scale="individual",verbose=0,label_size=.9) ################################################### ### code chunk number 20: pTV2 ################################################### cluster_results<-plotTV(dens.ind,regions=peaks.anno,control=c(dens.wt),show_names=F,norm_readc=F,verbose=0,label_size=.9) ################################################### ### code chunk number 21: pTV1_c (eval = FALSE) ################################################### ## cluster_results<-plotTV(dens.ind,dens.wt,regions=peaks.anno,show_names=T,norm_readc=F,scale="individual",verbose=0) ################################################### ### code chunk number 22: pTV2_c (eval = FALSE) ################################################### ## cluster_results<-plotTV(dens.ind,regions=peaks.anno,control=c(dens.wt),show_names=F,norm_readc=F,verbose=0) ################################################### ### code chunk number 23: pTV3 ################################################### genes2plot<-unique(mcols(gtf.dm3)$transcript_id) cluster_results<-plotTV(dens.pas_paired,dens.pas_upaired,regions=genes2plot,cluster=5,gtf=gtf.dm3,show_names=T,verbose=0,ex_windows=300) ################################################### ### code chunk number 24: tplots_1 ################################################### ngenes<-length(peaks.anno) fake.array<-matrix(rnorm(n=ngenes*8,mean=10,sd=2),nrow=ngenes,ncol=8,dimnames=list(paste(rep("Gene",ngenes),1:ngenes),paste("E",1:8,sep=""))) ################################################### ### code chunk number 25: pTV4 ################################################### cluster_results<-plotTV(fake.array,regions=peaks.anno,show_names=T,gclust="expression",cluster="hc_sp",label_size=0.7,verbose=0) ################################################### ### code chunk number 26: pTV3_c (eval = FALSE) ################################################### ## genes2plot<-unique(mcols(gtf.dm3)$transcript_id) ## cluster_results<-plotTV(dens.pas_paired,dens.pas_upaired,regions=genes2plot,gtf=gtf.dm3,show_names=T,cluster=5,verbose=0,ex_windows=300) ################################################### ### code chunk number 27: fplots_1_c (eval = FALSE) ################################################### ## ngenes<-length(peaks.anno) ## fake.array<-matrix(rnorm(n=ngenes*8,mean=10,sd=2),nrow=ngenes,ncol=8,dimnames=list(paste(rep("Gene",ngenes),1:ngenes),paste("E",1:8,sep=""))) ################################################### ### code chunk number 28: pTV4_c (eval = FALSE) ################################################### ## cluster_results2<-plotTV(fake.array,regions=peaks.anno,show_names=T,gclust="expression",cluster="hc_sp",label_size=0.7,verbose=0) ################################################### ### code chunk number 29: tplots_1_c ################################################### cluster_results ################################################### ### code chunk number 30: tplots_1_d (eval = FALSE) ################################################### ## cluster_order(cluster_results)[1:10] ## clusters(cluster_results)[1:10] ## summaryTV(cluster_results2) ################################################### ### code chunk number 31: tplots_1_e (eval = FALSE) ################################################### ## plotTV(dens.pas_paired,dens.pas_upaired,regions=genes2plot,gtf=gtf.dm3,show_names=T,rowv=cluster_results,verbose=0,ex_windows=300) ################################################### ### code chunk number 32: tplots_1_f ################################################### cluster_df<-plotTVData(cluster_results) ################################################### ### code chunk number 33: tplots_1_g (eval = FALSE) ################################################### ## ggplot(cluster_df,aes(x=Position,y=Average_scores,color=Sample))+geom_point()+facet_wrap(Plot~Cluster,scales="free",ncol=5) ################################################### ### code chunk number 34: tplots_1_e (eval = FALSE) ################################################### ## peak1.df<-meltPeak(dens.ind,dens.wt,region=peaks.tss["Peak.1"],peak_windows=100,rpm=F) ## ggplot(peak1.df,aes(x=Position,y=Reads,color=Label))+geom_line(size=2) ################################################### ### code chunk number 35: tplots_1_e (eval = FALSE) ################################################### ## peak1.df<-meltPeak(dens.ind,dens.wt,region=peaks.tss["Peak.1"],peak_windows=100,rpm=T) ## ggplot(peak1.df,aes(x=Position,y=NormalizedReads,color=Label))+geom_line(size=2)#Wrong! ################################################### ### code chunk number 36: sessionInfo ################################################### sessionInfo()