### R code from vignette source 'vignettes/AllelicImbalance/inst/doc/AllelicImbalance.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: intro
###################################################
library(AllelicImbalance)


###################################################
### code chunk number 2: createExampleASEset
###################################################

searchArea <- GRanges(seqnames = c("17"),ranges = IRanges(79478301,79478361))

pathToFiles <- system.file("extdata/ERP000101_subset", package="AllelicImbalance")

reads <- impBamGAL(pathToFiles,searchArea,verbose=FALSE)

heterozygotePositions <- scanForHeterozygotes(reads,verbose=FALSE)

countList <- getAlleleCounts(reads, heterozygotePositions, verbose=FALSE)

a.simple <- ASEsetFromCountList(heterozygotePositions,countList)

a.simple


###################################################
### code chunk number 3: gettingASEsetFromBcf
###################################################
BcfGR <- impBcfGR(pathToFiles,searchArea,verbose=FALSE)

countListBcf <- getAlleleCounts(reads, BcfGR,verbose=FALSE)

a.bcf <- ASEsetFromCountList(BcfGR, countListBcf)


###################################################
### code chunk number 4: creating stranded ASEset
###################################################
plus <- getAlleleCounts(reads, heterozygotePositions, strand="+",verbose=F) 
minus <- getAlleleCounts(reads, heterozygotePositions, strand="-",verbose=F)

a.stranded <-
ASEsetFromCountList(
heterozygotePositions,
countListPlus=plus,
countListMinus=minus
)
a.stranded


###################################################
### code chunk number 5: highlightgetAreaFromGeneNames
###################################################
#Getting searchArea from genesymbol
library(org.Hs.eg.db )
searchArea<-getAreaFromGeneNames("ACTG1",org.Hs.eg.db)

#Getting rs-IDs
library(SNPlocs.Hsapiens.dbSNP.20120608)
#seqlevels(a.simple) <- "chr17"
gr <- rowRanges(a.simple)
updatedGRanges<-getSnpIdFromLocation(gr, SNPlocs.Hsapiens.dbSNP.20120608)
rowRanges(a.simple)<-updatedGRanges



###################################################
### code chunk number 6: creatingphenotypedASEset
###################################################

#simulate phenotype data
pdata <- DataFrame(
	Treatment=sample(c("ChIP", "Input"),length(reads),replace=TRUE),
	Gender=sample(c("male", "female"),length(reads),replace=TRUE), 
	row.names=paste("individual",1:length(reads),sep=""))

#make new ASEset with pdata
a.new <- ASEsetFromCountList(
		heterozygotePositions,
		countList,
		colData=pdata)

#add to existing object
colData(a.simple) <- pdata



###################################################
### code chunk number 7: usingStatisticsTests
###################################################

#use a subset for tests
a2 <- a.stranded[,5:10]
#two types of tests
binom.test(a2,"+")
chisq.test(a2,"-")


###################################################
### code chunk number 8: plottingDemonstration1
###################################################
barplot(a.stranded[1],strand="+")

#use other test
btp <- binom.test(a.stranded[1],"+")

barplot(a.stranded[1],strand="+", testValue=btp)



###################################################
### code chunk number 9: plottingDemonstration2
###################################################
barplot(a.simple,type="fraction")


###################################################
### code chunk number 10: plottingDemonstration3
###################################################
sampleColour<-rep("palevioletred",ncol(a.simple))
sampleColour[colData(a.simple)[,"Gender"]%in%"male"] <- "blue"

barplot(a.simple[1],type="fraction",sampleColour=sampleColour)



###################################################
### code chunk number 11: useAnnotationPlot
###################################################
library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

barplot(a.simple[1],OrgDb=org.Hs.eg.db,TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)



###################################################
### code chunk number 12: locationPlot
###################################################
#using count type
locationplot(a.simple,type="count")

#use annotation
locationplot(a.simple,OrgDb=org.Hs.eg.db,TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)



###################################################
### code chunk number 13: sessioninfo
###################################################

sessionInfo()