## ------------------------------------------------------------------------
library(SPLINTER)

## ---- message=FALSE------------------------------------------------------
library(BSgenome.Mmusculus.UCSC.mm9)
library(GenomicFeatures)
bsgenome <- BSgenome.Mmusculus.UCSC.mm9

## ----message=FALSE-------------------------------------------------------
data_path<-system.file("extdata",package="SPLINTER")
gtf_file<-paste(data_path,"/Mus_musculus.Ensembl.NCBIM37.65.partial.gtf",sep="")
txdb <- makeTxDbFromGFF(file=gtf_file,chrominfo = Seqinfo(genome="mm9"))

# txdb generation can take quite long, you can save the object and load it the next time
# saveDb(txdb,file="txdb_hg19.sqlite")
# txdb<-loadDb(file="txdb_hg19.sqlite")

# extract CDS and exon information from TxDb object
thecds<-cdsBy(txdb,by="tx",use.names=TRUE)
theexons<-exonsBy(txdb,by="tx",use.names=TRUE)

## ------------------------------------------------------------------------
typeofAS="SE"
mats_file<-paste(data_path,"/skipped_exons.txt",sep="")
splice_data <-extractSpliceEvents(data=mats_file, filetype='mats', splicetype=typeofAS)
splice_data$data[,c(1:10)]

## ------------------------------------------------------------------------
splice_data<-addEnsemblAnnotation(data=splice_data,species="mmusculus")

# (Optional) Sorting the dataframe, if you have supporting statistical information
splice_data$data<-splice_data$data[with(splice_data$data,order(FDR,-IncLevelDifference)),]
head(splice_data$data[,c(1:10)])

## ------------------------------------------------------------------------
single_id='Prmt5'
pp<-which(grepl(single_id,splice_data$data$geneSymbol)) # Prmt5 has 1 record

splice_data$data[pp,c(1:6)] # show all records

single_record<-splice_data$data[pp[1],]

## ------------------------------------------------------------------------
valid_tx <- findTX(id=single_record$GeneID,tx=theexons,db=txdb)

valid_cds<- findTX(id=single_record$GeneID,tx=thecds,db=txdb)

## ------------------------------------------------------------------------
roi <- makeROI(splice_data,pp[1])
roi

## ------------------------------------------------------------------------
compatible_tx<- findCompatibleEvents(valid_tx,roi=roi,verbose=TRUE)

compatible_cds<- findCompatibleEvents(valid_cds,valid_tx,roi=roi,verbose=TRUE)


## ------------------------------------------------------------------------
region_minus_exon <-removeRegion(compatible_cds$hits[[1]],roi)

## ------------------------------------------------------------------------
# Not relevant for this Prmt5 skipped exon example
region_plus_exon <-insertRegion(region_minus_exon,roi)

## ------------------------------------------------------------------------
event<-eventOutcomeCompare(seq1=compatible_cds$hits[[1]],seq2=region_minus_exon,
                    genome=bsgenome,direction=TRUE,fullseq=FALSE)

event

## ------------------------------------------------------------------------
aa<-getRegionDNA(roi,bsgenome) 
aa

## ---- eval=FALSE---------------------------------------------------------
#  primers<-callPrimer3(seq=aa$seq,sequence_target = aa$jstart,size_range='100-500')

## ------------------------------------------------------------------------
primers[,c(1:4)]

## ------------------------------------------------------------------------
primers <- data.frame(PRIMER_LEFT_SEQUENCE="ACTTTCGGACTCTGTGTGACT",
                      PRIMER_RIGHT_SEQUENCE="TCATAGGCATTGGGTGGAGG",
                      stringsAsFactors=FALSE)

## ------------------------------------------------------------------------
cp<-checkPrimer(primers[1,],bsgenome,roi)

cp

## ------------------------------------------------------------------------
pcr_result1<-getPCRsizes(cp,theexons)
pcr_result1

tx_minus_exon <-removeRegion(compatible_tx$hits[[1]],roi)
pcr_result2<-getPCRsizes(cp,tx_minus_exon)
pcr_result2

## ------------------------------------------------------------------------
relevant_pcr_bands<-splitPCRhit(pcr_result1,compatible_tx)

relevant_pcr_bands

## ------------------------------------------------------------------------
# reading in BAM files
mt<-paste(data_path,"/mt_chr14.bam",sep="")
wt<-paste(data_path,"/wt_chr14.bam",sep="")

# Plotting genomic range, read density and splice changes
eventPlot(transcripts=theexons,roi_plot=roi,bams=c(wt,mt),names=c('wt','mt'),
          annoLabel=single_id,rspan=2000)

# Barplot of PSI values if provided
psiPlot(single_record)

## ----sessionInfo, echo=FALSE---------------------------------------------
sessionInfo()