WorkingWithAlignedNucleotides.R
In GenomicAlignments: Representation and manipulation of short genomic alignments

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

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: bamfiles
###################################################
library(RNAseqData.HNRNPC.bam.chr14)
bamfiles <- RNAseqData.HNRNPC.bam.chr14_BAMFILES
names(bamfiles)  # the names of the runs


###################################################
### code chunk number 3: quickBamFlagSummary
###################################################
library(Rsamtools)
quickBamFlagSummary(bamfiles[1], main.groups.only=TRUE)


###################################################
### code chunk number 4: ScanBamParam
###################################################
flag1 <- scanBamFlag(isFirstMateRead=TRUE, isSecondMateRead=FALSE,
                     isDuplicate=FALSE, isNotPassingQualityControls=FALSE)
param1 <- ScanBamParam(flag=flag1, what="seq")


###################################################
### code chunk number 5: readGAlignments
###################################################
library(GenomicAlignments)
gal1 <- readGAlignments(bamfiles[1], use.names=TRUE, param=param1)


###################################################
### code chunk number 6: read_sequences
###################################################
mcols(gal1)$seq


###################################################
### code chunk number 7: original-query-sequences
###################################################
oqseq1 <- mcols(gal1)$seq
is_on_minus <- as.logical(strand(gal1) == "-")
oqseq1[is_on_minus] <- reverseComplement(oqseq1[is_on_minus])


###################################################
### code chunk number 8: is_dup
###################################################
is_dup <- duplicated(names(gal1))
table(is_dup)


###################################################
### code chunk number 9: same-name-implies-same-seq-in-U1-oqseq
###################################################
dup2unq <- match(names(gal1), names(gal1))
stopifnot(all(oqseq1 == oqseq1[dup2unq]))


###################################################
### code chunk number 10: oqseq1
###################################################
oqseq1 <- oqseq1[!is_dup]


###################################################
### code chunk number 11: most_frequent_cigars
###################################################
head(sort(table(cigar(gal1)), decreasing=TRUE))


###################################################
### code chunk number 12: cigarOpTable
###################################################
colSums(cigarOpTable(cigar(gal1)))


###################################################
### code chunk number 13: table_njunc
###################################################
table(njunc(gal1))


###################################################
### code chunk number 14: readGAlignmentPairs
###################################################
library(pasillaBamSubset)
flag0 <- scanBamFlag(isDuplicate=FALSE, isNotPassingQualityControls=FALSE)
param0 <- ScanBamParam(flag=flag0)
U3.galp <- readGAlignmentPairs(untreated3_chr4(), use.names=TRUE, param=param0)
head(U3.galp)


###################################################
### code chunk number 15: first-and-last-U3.galp
###################################################
head(first(U3.galp))
head(last(U3.galp))


###################################################
### code chunk number 16: isProperPair
###################################################
table(isProperPair(U3.galp))


###################################################
### code chunk number 17: keep-only-proper-pairs
###################################################
U3.GALP <- U3.galp[isProperPair(U3.galp)]


###################################################
### code chunk number 18: U3.GALP_names_is_dup
###################################################
U3.GALP_names_is_dup <- duplicated(names(U3.GALP))
table(U3.GALP_names_is_dup)


###################################################
### code chunk number 19: U3.GALP_qnames
###################################################
U3.uqnames <- unique(names(U3.GALP))
U3.GALP_qnames <- factor(names(U3.GALP), levels=U3.uqnames)


###################################################
### code chunk number 20: U3.GALP_dup2unq
###################################################
U3.GALP_dup2unq <- match(U3.GALP_qnames, U3.GALP_qnames)


###################################################
### code chunk number 21: gaps-in-U3.GALP
###################################################
head(unique(cigar(first(U3.GALP))))
head(unique(cigar(last(U3.GALP))))
table(njunc(first(U3.GALP)), njunc(last(U3.GALP)))


###################################################
### code chunk number 22: no-indels-in-U3.GALP
###################################################
colSums(cigarOpTable(cigar(first(U3.GALP))))
colSums(cigarOpTable(cigar(last(U3.GALP))))


###################################################
### code chunk number 23: sessionInfo
###################################################
sessionInfo()