### 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) seqlevels(reads) <- "17" 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 <- rowData(a.simple) updatedGRanges<-getSnpIdFromLocation(gr, SNPlocs.Hsapiens.dbSNP.20120608) rowData(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()