plot get overridden rather to overlapped
1
0
Entering edit mode
mictadlo • 0
@mictadlo-10885
Last seen 20 months ago

I have the below code. Which executes twice a plot_denisities function. Unfortunately, when I execute this function again it plots a completely new picture rather overlap the two plots into one plot.

    #install if necessary
    source("http://bioconductor.org/biocLite.R")
    biocLite("Rsamtools")

    #load library
    library(Rsamtools)

    extracting_pos_neg_reads <- function(bam_fn) {

      #read in entire BAM file
      bam <- scanBam(bam_fn)

      #names of the BAM fields
      names(bam[[1]])
      # [1] "qname"  "flag"   "rname"  "strand" "pos"    "qwidth" "mapq"   "cigar"
      # [9] "mrnm"   "mpos"   "isize"  "seq"    "qual"

      #distribution of BAM flags
      table(bam[[1]]$flag)

      #      0       4      16
      #1472261  775200 1652949

      #function for collapsing the list of lists into a single list
      #as per the Rsamtools vignette
      .unlist <- function (x) {
        ## do.call(c, ...) coerces factor to integer, which is undesired
        x1 <- x[[1L]]
        if (is.factor(x1)) {
          structure(unlist(x), class = "factor", levels = levels(x1))
        } else {
          do.call(c, x)
        }
      }

      #store names of BAM fields
      bam_field <- names(bam[[1]])

      #go through each BAM field and unlist
      list <- lapply(bam_field, function(y)
        .unlist(lapply(bam, "[[", y)))

      #store as data frame
      bam_df <- do.call("DataFrame", list)
      names(bam_df) <- bam_field

      dim(bam_df)
      #[1] 3900410      13

      #---------

      #use chr22 as an example
      #how many entries on the negative strand of chr22?
      ###table(bam_df$rname == 'chr22' & bam_df$flag == 16)
      # FALSE    TRUE
      #3875997   24413

      #function for checking negative strand
      check_neg <- function(x) {
        if (intToBits(x)[5] == 1) {
          return(T)
        } else {
          return(F)
        }
      }

      #test neg function with subset of chr22
      test <- subset(bam_df)#, rname == 'chr22')
      dim(test)
      #[1] 56426    13
      table(apply(as.data.frame(test$flag), 1, check_neg))
      #number same as above
      #FALSE  TRUE
      #32013 24413

      #function for checking positive strand
      check_pos <- function(x) {
        if (intToBits(x)[3] == 1) {
          return(F)
        } else if (intToBits(x)[5] != 1) {
          return(T)
        } else {
          return(F)
        }
      }

      #check pos function
      table(apply(as.data.frame(test$flag), 1, check_pos))
      #looks OK
      #FALSE  TRUE
      #24413 32013

      #store the mapped positions on the plus and minus strands
      neg <- bam_df[apply(as.data.frame(bam_df$flag), 1, check_neg),
                    'pos']
      length(neg)
      #[1] 24413
      pos <- bam_df[apply(as.data.frame(bam_df$flag), 1, check_pos),
                    'pos']
      length(pos)
      #[1] 32013

      #calculate the densities
      neg_density <- density(neg)
      pos_density <- density(pos)

      #display the negative strand with negative values
      neg_density$y <- neg_density$y * -1

      return (list(neg_density, pos_density))
    }

    plot_densities <- function(density) {
      neg_density <- density[[1]]
      pos_density <- density[[2]]

      plot(
        pos_density,
        ylim = range(c(neg_density$y, pos_density$y)),
        main = "Coverage plot of Sample 5",
        xlab = "lenght 21",
        col = 'blue',
        type = 'h'
      )
      lines(neg_density, type = 'h', col = 'red')
    }


    filenames <- c("~/sample5-21.sam-uniq.sorted.bam", "~/sample5-24.sam-uniq.sorted.bam")

    for ( i in filenames){ 
      print(i)
      density <- extracting_pos_neg_reads(i)
      plot_densities(density)
    }

How could I overlap the two plots together?

Thank you in advance

rsamtools ggplot2 plot • 373 views
ADD COMMENT
0
Entering edit mode
@james-w-macdonald-5106
Last seen 21 hours ago
United States

This has nothing to do with any Bioconductor package, but instead is a question about base R (or more correctly your own code). You could ask on R-help (r-help@r-project.org), or you could do a Google search for something like 'multiple lines R plot', or you could look at the code for say plotDensities in limma.

ADD COMMENT

Login before adding your answer.

Traffic: 384 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6