密度错误.默认值(x = neg):找不到对象'neg' [英] Error in density.default(x = neg) : object 'neg' not found

问题描述

为了使用 ggplot2 ，我试图重写以下 plot_densities 功能.

  plot_densities<-函数(密度){neg_density<-密度[[1]]pos_density<-密度[[2]]阴谋(pos_density，ylim = range(c(neg_density $ y，pos_density $ y))，main =样本5的覆盖图"，xlab =长度21"，col ='blue'，类型='h')行数(neg_density，type ='h'，col ='red')} 

不幸的是，下面的新功能导致 density.default(x = neg)错误:未找到对象'neg'

  plot_densities2<-功能(密度){neg_density<-密度[[1]]pos_density<-密度[[2]]密度=附加(neg_density，pos_density)ggplot(as.data.frame(密度)，aes(x = x，y = y))+theme_bw()+geom_density(alpha = 0.5)} 

完整代码可在下面找到，数据可以从

I tried to rewrite the below plot_densities fuction in order to use ggplot2.

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')
}

Unfurtunately the new function below caused Error in density.default(x = neg) : object 'neg' not found

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

  densities = append(neg_density, pos_density)

  ggplot(as.data.frame(densities), aes(x=x, y=y)) + 
    theme_bw() +
    geom_density(alpha=0.5)
}

The full code can be found below and the data can be downloaded from here

#apt update && apt install zlib1g-dev

#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')
}


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

  densities = append(neg_density, pos_density)
  densities


  ggplot(as.data.frame(densities), aes(x=x, y=y)) + 
    theme_bw() +
    geom_density(alpha=0.5)
}

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_densities2(density)
}

解决方案

Density objects seem to be not the best ones to be used with append and as.data.frame. In particular, they contain some elements that caused problems but at the same time are unnecessary. What we may do is to pick only x and y elements as to construct the relevant data frame:

plot_densities2 <- function(density) {
  densities <- cbind(rbind(data.frame(density[[1]][1:2]), data.frame(density[[2]][1:2])), 
                     id = rep(c("neg", "pos"), each = length(density[[1]]$x)))
  print(ggplot(data = densities, aes(x = x, y = y, fill = id)) + 
          theme_bw() + geom_area(alpha = 0.5))
}

这篇关于密度错误.默认值(x = neg):找不到对象'neg'的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！