count_alignments: Get alignment counts from BAM files

In e-myers/rnaseq: Process, Analyze and Visualize RNA-seq Data

Description

Given a list of BAM files, return an array of either primary (most likely) alignments or all alignments that passed QC.

Usage

count_alignments(inFiles, samtoolsPath = "~/anaconda2/bin/",
  what = "primary", verbose = TRUE, outFile = "alignmentCounts.txt")

Arguments

inFiles	Character - BAM file list
samtoolsPath	String - path to samtools directory
what	String - What type of reads to include
verbose	Logical - whether to print progress / results
outFile	String - Name of file to write counts to, or NULL (this function can take awhile, and there's no safeguard against overwriting this file here)

Details

Might make sense to give it an option to include, say, multi-mapped alignments. Requires samtools to be installed. For explanations of samtools flags (like "-f 0" and "-F 0x904" below), see this tool: http://broadinstitute.github.io/picard/explain-flags.html Note that "-f" means to include only things with these flags, and "-F" means to exclude anything with these flags. TIME: ~1m / BAM.

Value

Numeric vector? - Alignment count for each BAM

Author(s)

Emma Myers

Examples

readCounts = count_alignments(bamFileNames)
pdf("countsBar.pdf")
barplot(readCounts/1000000, las=2, main="Primary alignment counts", cex.names = 0.5, names.arg=gsub("_mapped_Aligned.sortedByCoord.out.bam", "", bamFileNames), ylab="Millions")
dev.off()

e-myers/rnaseq documentation built on May 20, 2019, 9:14 p.m.