# NOT RUN {
# Open (indexed) BAM file
bam<-system.file("extdata", "accepted_hits.bam", package="rbamtools")
reader<-bamReader(bam,idx=TRUE)
# Extract reads from BAM file
xlim = c(10000, 30000)
coords<-c(0,xlim[1], xlim[2])
range<-bamRange(reader,coords)
bamClose(reader)
# Calculate align depth
ad<-alignDepth(range)
ad
getParams(ad)
# Prepare plotting parameter
gene<-"WASH7P"
ensg_id <- "ENSG00000227232"
enst_id <- "ENST00000538476"
# Get exon positions
start<-c(14411,15000,15796,15904,16607,16748,16858,17233,17602,17915,18268,24737,29534)
end <- c(14502,15038,15901,15947,16745,16765,17055,17364,17742,18061,18366,24891,29806)
# Do plot
plotAlignDepth(ad, lwd = 2, xlim = xlim,
main = paste("Align depth for gene",gene),
ylab = "Align depth", start = start,
end = end, strand = "-",
transcript = paste("Chromosome 1",
"\tGene", ensg_id,
"\tTranscript ",enst_id)
)
# }
Run the code above in your browser using DataLab