inst/doc/TransView.R

In TransView: Read density map construction and accession. Visualization of ChIPSeq and RNASeq data sets

### 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()