flagstat-collect-info.R

# flagstat.fpath <- "~/testing/br-mgd-srt-srt.txt"

# content of the file (flagstat output):

51712430 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
46380832 + 0 mapped (89.69% : N/A)
51712430 + 0 paired in sequencing
25856215 + 0 read1
25856215 + 0 read2
38000240 + 0 properly paired (73.48% : N/A)
44578824 + 0 with itself and mate mapped
1802008 + 0 singletons (3.48% : N/A)
449718 + 0 with mate mapped to a different chr
408844 + 0 with mate mapped to a different chr (mapQ>=5)

##############################################



extract.flagstat.info <- function(flagstat.fpath) {
  flagstat.lines <- readLines(flagstat.fpath)

  extract.numbers <- function(flagstat.line) {
    regmatches(flagstat.line,
               regexec("(\\d+) \\+ (\\d+) (.*? *\\((\\d+\\.\\d+)% : (N/A|\\d+)\\)*.*|.*)",
                       flagstat.line,
                       perl = TRUE)
    )[[1]][c(2, 3, 5, 6)]
  }

  flagstat.number.mat <- Reduce(rbind, lapply(flagstat.lines, extract.numbers))
  flagstat.number.mat <- tapply(as.numeric, flagstat.number.mat)
  print(flagstat.number.mat)
  tot.reads <- flagstat.number.mat[1, 1] + flagstat.number.mat[1, 2]
  total.qc.passed.reads <- flagstat.number.mat[1, 1]
  all.mapped.reads <- flagstat.number.mat[5, 1]
  all.mapped.as.pair.reads <- flagstat.number.mat[10, 1]
  all.mapped.reads.percent <- flagstat.number.mat[5, 3]
  all.mapped.as.pair.percent <- round(flagstat.number.mat[10, 1]/tot.reads, digits = 2)
  duplicates <- flagstat.number.mat[4, 1]
  res <- c(tot.reads,total.qc.passed.reads, all.mapped.reads,
           all.mapped.reads.percent, all.mapped.reads.percent, duplicates)
  names(res) <- c("total.reads", "qc-passed", "mapped", "mapped-%", "mapped.as.pair", "mapped.as.pair-%", "duplicates")
  res
}

# usage:
# extract.flagstat.info(flagstat.fpath)