### R code from vignette source 'vignettes/spliceSites/inst/doc/spliceSites.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: spliceSites.Rnw:151-158 ################################################### library(spliceSites) bam<-character(2) bam[1]<-system.file("extdata","rna_fem.bam",package="spliceSites") bam[2]<-system.file("extdata","rna_mal.bam",package="spliceSites") reader<-bamReader(bam[1],idx=TRUE) gs<-getGapSites(reader,seqid=1) gs ################################################### ### code chunk number 2: spliceSites.Rnw:163-165 ################################################### ga<-alignGapList(reader) ga ################################################### ### code chunk number 3: spliceSites.Rnw:172-174 ################################################### mbg<-readMergedBamGaps(bam) mbg ################################################### ### code chunk number 4: spliceSites.Rnw:178-182 ################################################### prof<-data.frame(gender=c("f","m")) rtbg<-readTabledBamGaps(bam,prof=prof,rpmg=TRUE) rtbg getProfile(rtbg) ################################################### ### code chunk number 5: spliceSites.Rnw:216-229 ################################################### ljp<-lJunc(ga,featlen=6,gaplen=6,strand='+') ljp ljm<-lJunc(ga,featlen=6,gaplen=6,strand='-') ljm rjp<-rJunc(ga,featlen=6,gaplen=6,strand='+') rjp rjm<-rJunc(ga,featlen=6,gaplen=6,strand='-') rjm lrjp<-lrJunc(ga,lfeatlen=6,rfeatlen=6,strand='+') lrjp lrjm<-lrJunc(ga,lfeatlen=6,rfeatlen=6,strand='-') lrjm ################################################### ### code chunk number 6: spliceSites.Rnw:243-248 ################################################### ljp1<-lCodons(ljp,frame=1) ljp1 ljp2<-lCodons(ljp,frame=2) rjm1<-rCodons(ljm,frame=1) rjm2<-rCodons(ljm,frame=2) ################################################### ### code chunk number 7: spliceSites.Rnw:254-257 ################################################### lr1<-lrCodons(lrjp,frame=1,strand='+') lr2<-lrCodons(lrjp,frame=2,strand='+') lr3<-lrCodons(lrjp,frame=3,strand='+') ################################################### ### code chunk number 8: spliceSites.Rnw:263-266 ################################################### ljpc<-c(ljp1,ljp2) rjmc<-c(rjm1,rjm2) lrj<-c(lr1,lr2,lr3) ################################################### ### code chunk number 9: spliceSites.Rnw:270-273 ################################################### ljpc<-sortTable(ljpc) rjmc<-sortTable(rjmc) lrj<-sortTable(lrj) ################################################### ### code chunk number 10: spliceSites.Rnw:279-283 ################################################### trim_left(lrj,3) trim_right(lrj,3) resize_left(lrj,8) resize_right(lrj,8) ################################################### ### code chunk number 11: spliceSites.Rnw:294-297 ################################################### ucf<-system.file("extdata","uc_small.RData",package="spliceSites") uc<-loadGenome(ucf) annotation(ga)<-annotate(ga,uc) ################################################### ### code chunk number 12: spliceSites.Rnw:301-302 ################################################### annotation(rtbg)<-annotate(rtbg,uc) ################################################### ### code chunk number 13: spliceSites.Rnw:307-308 ################################################### strand(ga)<-getAnnStrand(ga) ################################################### ### code chunk number 14: spliceSites.Rnw:312-314 ################################################### gap<-addGeneAlignPart(ga) gap ################################################### ### code chunk number 15: spliceSites.Rnw:321-324 ################################################### dnafile<-system.file("extdata","dna_small.RData",package="spliceSites") load(dnafile) dna_small ################################################### ### code chunk number 16: spliceSites.Rnw:328-330 ################################################### ljpcd<-dnaRanges(ljpc,dna_small) seqlogo(ljpcd) ################################################### ### code chunk number 17: spliceSites.Rnw:334-335 ################################################### ljpca<-translate(ljpcd) ################################################### ### code chunk number 18: spliceSites.Rnw:342-347 ################################################### # A) For gapSites extractRange(ga,seqid="chr1",start=14000,end=30000) # B) For cRanges lj<-lJunc(ga,featlen=3,gaplen=6,strand='+') extractRange(lj,seqid="chr1",start=14000,end=30000) ################################################### ### code chunk number 19: spliceSites.Rnw:352-357 ################################################### lj<-lJunc(ga,featlen=6,gaplen=3,strand='+') ljw<-extractByGeneName(ljp,geneNames="POLR3K",src=uc) ljw ljw<-extractByGeneName(ljpcd,geneNames="POLR3K",src=uc) ljw ################################################### ### code chunk number 20: spliceSites.Rnw:364-372 ################################################### l<-12 lj<-lJunc(mbg,featlen=l,gaplen=l,strand='+') ljc<-c(lCodons(lj,1),lCodons(lj,2),lCodons(lj,3)) lrj<-lrJunc(mbg,lfeatlen=l,rfeatlen=l,strand='+') lrjc<-c(lrCodons(lrj,1),lrCodons(lrj,2),lrCodons(lrj,3)) jlrd<-dnaRanges(ljc,dna_small) ################################################### ### code chunk number 21: spliceSites.Rnw:376-378 ################################################### jlrt<-translate(jlrd) jlrt ################################################### ### code chunk number 22: spliceSites.Rnw:384-386 ################################################### jlrtt<-truncateSeq(jlrt) jlrtt ################################################### ### code chunk number 23: spliceSites.Rnw:391-393 ################################################### jtry<-trypsinCleave(jlrtt) jtry<-sortTable(jtry) ################################################### ### code chunk number 24: spliceSites.Rnw:402-404 (eval = FALSE) ################################################### ## annotation(rtbg)<-annotate(rtbg,uc) ## write.annDNA.tables(rtbg,dna_small,"gapSites.csv",featlen=3,gaplen=8) ################################################### ### code chunk number 25: spliceSites.Rnw:409-410 (eval = FALSE) ################################################### ## write.files(jtry,path=getwd(),filename="proteomic",quote=FALSE) ################################################### ### code chunk number 26: spliceSites.Rnw:434-435 ################################################### al<-alt_left_ranks(ga) ################################################### ### code chunk number 27: spliceSites.Rnw:440-441 ################################################### ar<-alt_ranks(ga) ################################################### ### code chunk number 28: spliceSites.Rnw:445-446 ################################################### plot_diff_ranks(ga) ################################################### ### code chunk number 29: spliceSites.Rnw:450-452 ################################################### aga<-annotation(ga) plot_diff(aga) ################################################### ### code chunk number 30: spliceSites.Rnw:459-471 ################################################### mes<-load.maxEnt() score5(mes,"CCGGGTAAGAA",4) score3(mes,"CTCTACTACTATCTATCTAGATC",pos=20) sq5<-scoreSeq5(mes,seq="ACGGTAAGTCAGGTAAGT") sq3<-scoreSeq3(mes,seq="TTTATTTTTCTCACTTTTAGAGACTTCATTCTTTCTCAAATAGGTT") gae<-addMaxEnt(ga,dna_small,mes) gae table(getMeStrand(gae)) sae<-setMeStrand(gae) sae ################################################### ### code chunk number 31: spliceSites.Rnw:475-479 ################################################### # hb<-load.hbond() seq<-c("CAGGTGAGTTC","ATGCTGGAGAA","AGGGTGCGGGC","AAGGTAACGTC","AAGGTGAGTTC") hbond(hb,seq,3) ################################################### ### code chunk number 32: spliceSites.Rnw:482-487 ################################################### gab<-addHbond(ga,dna_small) # D) cdRanges lj<-lJunc(ga,featlen=3,gaplen=8,strand='+') ljd<-dnaRanges(lj,dna_small) ljdh<-addHbond(ljd) ################################################### ### code chunk number 33: spliceSites.Rnw:497-504 ################################################### prof<-data.frame(gender=c("f","m")) rtbg<-readTabledBamGaps(bam,prof=prof,rpmg=TRUE) getProfile(rtbg) meta<-data.frame(labelDescription=names(prof),row.names=names(prof)) pd<-new("AnnotatedDataFrame",data=prof,varMetadata=meta) es<-readExpSet(bam,phenoData=pd) ################################################### ### code chunk number 34: spliceSites.Rnw:508-511 ################################################### ann<-annotate(es,uc) ucj<-getSpliceTable(uc) uja<-uniqueJuncAnn(es,ucj) ################################################### ### code chunk number 35: spliceSites.Rnw:517-524 ################################################### n<-10 cuff<-system.file("extdata","cuff_files",paste(1:n,"genes","fpkm_tracking",sep="."),package="spliceSites") gr<-system.file("extdata","cuff_files","groups.csv",package="spliceSites") groups<-read.table(gr,sep="\t",header=TRUE) meta<-data.frame(labelDescription=c("gender","age-group","location"),row.names=c("gen","agg","loc")) phenoData<-new("AnnotatedDataFrame",data=groups,varMetadata=meta) exset<-readCuffGeneFpkm(cuff,phenoData) ################################################### ### code chunk number 36: spliceSites.Rnw:537-541 ################################################### prof<-data.frame(gen=factor(c("w","m","w","w"),levels=c("m","w")), loc=factor(c("thx","thx","abd","abd"),levels=c("thx","abd")), ag =factor(c("y","y","m","o"),levels=c("y","m","o"))) prof ################################################### ### code chunk number 37: spliceSites.Rnw:545-569 ################################################### key1<-data.frame(id=1:5, seqid=c(1,1,2,2,3), lend=c(10,20,10,30,10), rstart=c(20,30,20,40,20), nAligns=c(11,21,31,41,51)) key2<-data.frame(id=1:5, seqid=c(1,1,2,2,4), lend=c(10,20,10,30,50), rstart=c(20,30,20,40,70), nAligns=c(21,22,23,24,25)) key3<-data.frame(id=1:5, seqid=c(1,2,4,5,5), lend=c(10,10,60,10,20), rstart=c(20,20,80,20,30), nAligns=c(31,32,33,34,35)) key<-rbind(key1,key2,key3) # Group positions ku<-aggregate(data.frame(nAligns=key$nAligns), by=list(seqid=key$seqid,lend=key$lend,rstart=key$rstart), FUN=sum) ################################################### ### code chunk number 38: spliceSites.Rnw:573-578 ################################################### # Count probes kpc<-new("keyProfiler",keyTable=key1[,c("seqid","lend","rstart")],prof=prof) addKeyTable(kpc,keyTable=key2[,c("seqid","lend","rstart")],index=2) addKeyTable(kpc,keyTable=key3[,c("seqid","lend","rstart")],index=4) ################################################### ### code chunk number 39: spliceSites.Rnw:582-584 ################################################### cp<-appendKeyTable(kpc,ku,prefix="c.") cp ################################################### ### code chunk number 40: spliceSites.Rnw:588-595 ################################################### # Count aligns kpa<-new("keyProfiler",keyTable=key1[,c("seqid","lend","rstart")],prof=prof,values=key1$nAligns) kpa@ev$dtb addKeyTable(kpa,keyTable=key2[,c("seqid","lend","rstart")],index=2,values=key2$nAligns) addKeyTable(kpa,keyTable=key3[,c("seqid","lend","rstart")],index=4,values=key3$nAligns) ca<-appendKeyTable(kpa,ku,prefix="aln.") ca