TransWikia.com

ATAC seq : Plotting insert sizes

Bioinformatics Asked on January 26, 2021

My tutorial link

library(GenomicAlignments)

atacReads <- readGAlignmentPairs("MACS2_peak_call/ALL_ATAC_data/HSC1.bam", param = ScanBamParam(mapqFilter = 1, 
            flag = scanBamFlag(isPaired = TRUE, isProperPair = TRUE), what = c("qname", 
        "mapq", "isize"), which = GRanges("chr20", IRanges(1, 63025520))))
# length(atacReads)
atacReads


atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)


library(magrittr)
library(dplyr)
library(ggplot2)
fragLenPlot <- table(insertSizes) %>% data.frame %>% dplyr::rename(InsertSize = insertSizes, 
    Count = Freq) %>% mutate(InsertSize = as.numeric(as.vector(InsertSize)), 
 Count = as.numeric(as.vector(Count))) %>% ggplot(aes(x = InsertSize, y = Count)) + 
  geom_line()

fragLenPlot + theme_bw()

fragLenPlot + scale_y_continuous(trans = "log2") + theme_bw()

How to i specify the Grange to span the total chromosome. Now I use only "chr20" as shown in the example.

One Answer

You can exclude the which statement and get the entirety of all of the chromosomes. Note however that R is a very bad platform to use for this, since you end up reading everything into memory (I hope the BAM file isn't huge).

Note that you're seriously reinventing the wheel. There are a number of packages that will directly produce the plot you want, such as bamPEFragmentSize from deepTools.

Correct answer by Devon Ryan on January 26, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP