inst/doc/customProDB.R

In customProDB: Generate customized protein database from NGS data, with a focus on RNA-Seq data, for proteomics search

### R code from vignette source 'customProDB.Rnw'

###################################################
### code chunk number 1: dbSNPversion
###################################################
library('rtracklayer')
session <- browserSession()
genome(session) <- 'hg19'
dbsnps <- trackNames(session)[grep('snp', trackNames(session), fixed=T)]
dbsnps


###################################################
### code chunk number 2: options
###################################################
options(width=70)


###################################################
### code chunk number 3: loadpkg
###################################################
library(customProDB)


###################################################
### code chunk number 4: PrepareAnnoRef
###################################################
transcript_ids <- c("NM_001126112", "NM_033360", "NR_073499", "NM_004448", 
			"NM_000179", "NR_029605", "NM_004333", "NM_001127511")
pepfasta <- system.file("extdata", "refseq_pro_seq.fasta", 
				package="customProDB")
CDSfasta <- system.file("extdata", "refseq_coding_seq.fasta",
				package="customProDB")
annotation_path <- tempdir()
PrepareAnnotationRefseq(genome='hg19', CDSfasta, pepfasta, annotation_path, 
		dbsnp = NULL, transcript_ids=transcript_ids, 
        splice_matrix=FALSE, ClinVar=FALSE)


###################################################
### code chunk number 5: PrepareAnnoENSEMBL (eval = FALSE)
###################################################
## ensembl <- useMart("ENSEMBL_MART_ENSEMBL", dataset="hsapiens_gene_ensembl",
##     host="sep2015.archive.ensembl.org", path="/biomart/martservice", 
##     archive=FALSE)
## annotation_path <- tempdir()
## transcript_ids <- c("ENST00000234420", "ENST00000269305", "ENST00000445888", 
## 			"ENST00000257430", "ENST00000508376", "ENST00000288602",
## 			"ENST00000269571", "ENST00000256078", "ENST00000384871")
## PrepareAnnotationEnsembl(mart=ensembl, annotation_path=annotation_path, 
## 			splice_matrix=FALSE, dbsnp=NULL, 
##             transcript_ids=transcript_ids, COSMIC=FALSE)


###################################################
### code chunk number 6: calculate
###################################################
load(system.file("extdata/refseq", "exon_anno.RData", package="customProDB"))
bamFile <- system.file("extdata/bams", "test1_sort.bam", package="customProDB")
load(system.file("extdata/refseq", "ids.RData", package="customProDB"))
RPKM <- calculateRPKM(bamFile, exon, proteincodingonly=TRUE, ids)


###################################################
### code chunk number 7: ouputpro
###################################################
load(system.file("extdata/refseq", "proseq.RData", package="customProDB"))
outf1 <- paste(tempdir(), '/test_rpkm.fasta', sep='')
Outputproseq(RPKM, 1, proteinseq, outf1, ids)


###################################################
### code chunk number 8: inpvcf
###################################################
# single sample
vcffile <- system.file("extdata/vcfs", "test1.vcf", package="customProDB")
vcf <- InputVcf(vcffile)
length(vcf)
vcf[[1]][1:3]


###################################################
### code chunk number 9: inpvcfm
###################################################
# multiple samples in one VCF file
vcffile <- system.file("extdata", "test_mul.vcf", package="customProDB")
vcfs <- InputVcf(vcffile)


###################################################
### code chunk number 10: location
###################################################
table(values(vcf[[1]])[['INDEL']])
index <- which(values(vcf[[1]])[['INDEL']]==TRUE)
indelvcf <- vcf[[1]][index]

index <- which(values(vcf[[1]])[['INDEL']]==FALSE)
SNVvcf <- vcf[[1]][index]

load(system.file("extdata/refseq", "ids.RData", package="customProDB"))
txdb <- loadDb(system.file("extdata/refseq", "txdb.sqlite", package="customProDB"))
SNVloc <- Varlocation(SNVvcf,txdb,ids)
indelloc <- Varlocation(indelvcf,txdb,ids)
table(SNVloc[,'location'])


###################################################
### code chunk number 11: poscoding
###################################################
load(system.file("extdata/refseq", "exon_anno.RData", package="customProDB"))
load(system.file("extdata/refseq", "dbsnpinCoding.RData", package="customProDB"))
load(system.file("extdata/refseq", "cosmic.RData", package="customProDB"))
postable_snv <- Positionincoding(SNVvcf, exon, dbsnpinCoding, COSMIC=cosmic)
postable_snv
postable_indel <- Positionincoding(indelvcf, exon)
postable_indel


###################################################
### code chunk number 12: aavar
###################################################
load(system.file("extdata/refseq", "procodingseq.RData", package="customProDB"))
txlist <- unique(postable_snv[, 'txid'])
codingseq <- procodingseq[procodingseq[, 'tx_id'] %in% txlist,]
mtab <- aaVariation (postable_snv, codingseq)
mtab


###################################################
### code chunk number 13: outvarpro
###################################################
outfile <-  paste(tempdir(), '/test_snv.fasta', sep='')
load(system.file("extdata/refseq", "proseq.RData", package="customProDB"))
OutputVarproseq(mtab, proteinseq, outfile, ids)


###################################################
### code chunk number 14: outabber
###################################################
txlist_indel <- unique(postable_indel[, 'txid'])
codingseq_indel <- procodingseq[procodingseq[, 'tx_id'] %in% txlist_indel, ]
outfile <-  paste(tempdir(), '/test_indel.fasta', sep='')
Outputaberrant(postable_indel, coding=codingseq_indel, proteinseq=proteinseq, 
			outfile=outfile, ids=ids)


###################################################
### code chunk number 15: junctype
###################################################
bedfile <- system.file("extdata/beds", "junctions1.bed", package="customProDB")
jun <-  Bed2Range(bedfile,skip=1,covfilter=5)
jun

load(system.file("extdata/refseq", "splicemax.RData", package="customProDB"))
load(system.file("extdata/refseq", "ids.RData", package="customProDB"))
junction_type <- JunctionType(jun, splicemax, txdb, ids) 
junction_type[10:19,]
table(junction_type[, 'jun_type'])


###################################################
### code chunk number 16: novjunc
###################################################
outf_junc <- paste(tempdir(), '/test_junc.fasta',sep='')
library('BSgenome.Hsapiens.UCSC.hg19')
OutputNovelJun <- OutputNovelJun(junction_type, Hsapiens, outf_junc, 
        proteinseq)


###################################################
### code chunk number 17: sharedPrO
###################################################
path <- system.file("extdata/bams", package="customProDB")
bamFile<- paste(path, '/', list.files(path,pattern="*bam$"), sep='')
rpkms <- sapply(bamFile, function(x) 
			calculateRPKM(x, exon, proteincodingonly=TRUE, ids))
#colnames(rpkms) <- c('1', '2', '3')
#rpkms
outfile <- paste(tempdir(), '/test_rpkm_share.fasta', sep='')
pro <- OutputsharedPro(rpkms, cutoff=1, share_sample=2, proteinseq, 
        outfile, ids)


###################################################
### code chunk number 18: mulvcf
###################################################
path <- system.file("extdata/vcfs", package="customProDB")
vcfFiles<- paste(path, '/', list.files(path, pattern="*vcf$"), sep='')
vcfs <- lapply(vcfFiles, function(x) InputVcf(x))
shared <- Multiple_VCF(vcfs, share_num=2)
shared


###################################################
### code chunk number 19: muljunc
###################################################
path <- system.file("extdata/beds", package="customProDB")
bedFiles<- paste(path, '/', list.files(path, pattern="*bed$"), sep='')
juncs <- lapply(bedFiles, function(x) Bed2Range(x, skip=1, covfilter=5))
sharedjun <- SharedJunc(juncs, share_num=2, ext_up=100, ext_down=100)
sharedjun


###################################################
### code chunk number 20: easyr
###################################################
bamFile <- system.file("extdata/bams", "test1_sort.bam", 
			package="customProDB")
vcffile <- system.file("extdata/vcfs", "test1.vcf", package="customProDB")
bedfile <- system.file("extdata", "junctions.bed", package="customProDB")
annotation_path <- system.file("extdata/refseq", package="customProDB")
outfile_path <- tempdir()
outfile_name='test'
easyRun(bamFile, RPKM=NULL, vcffile, annotation_path, outfile_path, 
        outfile_name, rpkm_cutoff=1, INDEL=TRUE, lablersid=TRUE, COSMIC=TRUE, 
        nov_junction=FALSE) 


###################################################
### code chunk number 21: easyrmul
###################################################
bampath <- system.file("extdata/bams", package="customProDB")
vcfFile_path <- system.file("extdata/vcfs",  package="customProDB")
annotation_path <- system.file("extdata/refseq", package="customProDB")
outfile_path <- tempdir()	
outfile_name <- 'mult'
easyRun_mul(bampath, RPKM_mtx=NULL, vcfFile_path, annotation_path, rpkm_cutoff=1,
            share_num=2, var_shar_num=2, outfile_path, outfile_name, INDEL=TRUE,
            lablersid=TRUE, COSMIC=TRUE, nov_junction=FALSE)


###################################################
### code chunk number 22: FASTAnor
###################################################
outfile_path <- system.file("extdata/tmp", package="customProDB")
readLines(file(paste(outfile_path, '/test_rpkm.fasta', sep=''), 'rt'), 1)


###################################################
### code chunk number 23: FASTAsnv
###################################################
readLines(file(paste(outfile_path, '/test_snv.fasta', sep=''), 'rt'), 1)


###################################################
### code chunk number 24: FASTAindel
###################################################
readLines(file(paste(outfile_path, '/test_indel.fasta', sep=''), 'rt'), 1)


###################################################
### code chunk number 25: FASTAjun
###################################################
readLines(file(paste(outfile_path, '/test_junc.fasta', sep=''), 'rt'), 1)


###################################################
### code chunk number 26: SessionInfo
###################################################
sessionInfo()