bowtie2_qc.R 3.09 KB
Newer Older
Holger Brandl's avatar
Holger Brandl committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
#!/usr/bin/env Rscript
#+ echo=FALSE, message=FALSE

suppressMessages(library(docopt))

# retrieve and parse the command-line arguments
doc <- '
Create a small summary for bam-files
Usage: bams_qc.R <base_directory>

Options:
-c        Peform correlation analysis
'

opts <- docopt(doc, commandArgs(TRUE))
#opts <- docopt(doc, ".")

if(!exists("opts")){ stop(doc); return }
#print(opts)

baseDir <- normalizePath(opts$base_directory)

if(is.na(file.info(baseDir)$isdir)){
    stop(paste("base directory", baseDir, "does not exist"))
}

27 28
devtools::source_url("https://git.mpi-cbg.de/bioinfo/datautils/raw/v1.26/R/core_commons.R")
#devtools::source_url("https://git.mpi-cbg.de/bioinfo/datautils/raw/v1.26/R/ggplot_commons.R")
Holger Brandl's avatar
Holger Brandl committed
29 30 31 32 33 34 35 36 37


########################################################################################################################
#' # Mapping Summary for: `r baseDir`



logSuffix=".bowtie.log"
parseAlgnSummary <- function(alignSummary){
Holger Brandl's avatar
Holger Brandl committed
38 39 40 41
    # pe example
    #alignSummary="/lustre/projects/plantx/smed_genome/schMed2/alignments_SMED_P297/CP-1318_S1_L001.bowtie.logs"
    # se example
    # alignSummary= "/lustre/projects/plantx/smed_genome/schMed2/alignments_SMED_P297/SRR959562.bowtie.logs"
Holger Brandl's avatar
Holger Brandl committed
42 43
    algnData <- readLines(alignSummary)

Holger Brandl's avatar
Holger Brandl committed
44 45 46 47 48 49 50
    findLine <- function(lines, pattern)lines[str_detect(lines, pattern)] %>% str_trim()
#    findLine(algnData, "overall alignment rate")


    multiLine = findLine(algnData, "aligned >1 times")
    uniqueLine = findLine(algnData, "aligned exactly 1 time")

Holger Brandl's avatar
Holger Brandl committed
51
    data.frame(
Holger Brandl's avatar
Holger Brandl committed
52
        sample=trimEnd(basename(alignSummary), logSuffix),
Holger Brandl's avatar
Holger Brandl committed
53
        num_reads=as.numeric(str_split_fixed(algnData[1], " ", 2)[1]),
Holger Brandl's avatar
Holger Brandl committed
54
        mapping_efficiency=as.numeric(str_replace(str_split_fixed(findLine(algnData, "overall alignment rate"), " ", 2)[1], "%", "")),
Holger Brandl's avatar
Holger Brandl committed
55

Holger Brandl's avatar
Holger Brandl committed
56 57
        unique_mappers=as.numeric(str_split_fixed(uniqueLine, " ", 2)[1]),
        unique_mapper_prop=str_match(uniqueLine, "[(]([0-9.]*)[%)]*")  %>% subset(select=2) %>% as.numeric(),
Holger Brandl's avatar
Holger Brandl committed
58

Holger Brandl's avatar
Holger Brandl committed
59 60
        multi_mappers=as.numeric(str_split_fixed(multiLine, " ", 2)[1]),
        multi_mappers_prop=str_match(multiLine, "[(]([0-9.]*)[%)]*")  %>% subset(select=2) %>% as.numeric()
Holger Brandl's avatar
Holger Brandl committed
61 62 63 64 65 66 67 68 69 70
    )
}

algnSummary <- list.files(baseDir, logSuffix, full.names=TRUE, recursive=T) %>% ldply(parseAlgnSummary)
write.delim(algnSummary, file="bowtie2_qc.txt")

scale_fill_discrete <- function (...){ scale_color_brewer(..., type = "seq", palette="Set1", "fill", na.value = "grey50") }


#+ fig.height=nrow(algnSummary)/3
Holger Brandl's avatar
Holger Brandl committed
71
ggplot(algnSummary, aes(sample, mapping_efficiency)) +
Holger Brandl's avatar
Holger Brandl committed
72 73 74 75 76
    geom_bar(stat="identity") +
    coord_flip() +
    ylim(0,100) +
    ggtitle("mapping efficiency")

Holger Brandl's avatar
Holger Brandl committed
77
ggplot(algnSummary, aes(sample, num_reads)) +
Holger Brandl's avatar
Holger Brandl committed
78 79 80 81
    geom_bar(stat="identity") +
    coord_flip() +
    ggtitle("read counts") +scale_y_continuous(labels=comma)

Holger Brandl's avatar
Holger Brandl committed
82
ggplot(algnSummary, aes(sample, unique_mapper_prop)) +
Holger Brandl's avatar
Holger Brandl committed
83 84 85 86
    geom_bar(stat="identity") + coord_flip() +
    ggtitle("unique-mapper proportions") +
    scale_y_continuous(labels=comma)

Holger Brandl's avatar
Holger Brandl committed
87
ggplot(algnSummary, aes(sample, multi_mappers_prop)) +
Holger Brandl's avatar
Holger Brandl committed
88 89 90 91
    geom_bar(stat="identity") + coord_flip() +
    ggtitle("multi-mapper proportions") +
    scale_y_continuous(labels=comma)