My tutorial link


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_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)

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)) + 

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.

  • 1
    $\begingroup$ sorry don't quite get the question. is it about knowing the chromosome length from the bam file, or it's something to do with the plotting $\endgroup$
    – StupidWolf
    Jan 6, 2021 at 11:09
  • $\begingroup$ Now Im too bit confused ,after reading paper where they don't mention which chromosome they taking and plotting. It would be essentially this one " is it about knowing the chromosome length from the bam file" can i like "chr1-chr20" and see the insert size distribution . I was not able to figure out how to define more that one chromosome in that Grange argument // $\endgroup$
    – kcm
    Jan 6, 2021 at 11:24
  • $\begingroup$ @StupidWolf for example in this paper nature.com/articles/s41597-019-0071-0#Sec1 Fig2(b). $\endgroup$
    – kcm
    Jan 6, 2021 at 11:26

1 Answer 1


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.

  • $\begingroup$ "bamPEFragmentSize" will see this $\endgroup$
    – kcm
    Jan 6, 2021 at 14:59
  • $\begingroup$ " you're seriously reinventing the wheel." yes since I didn;t read about your tool "bamPEFragmentSize" I used it , Works absolutely fine $\endgroup$
    – kcm
    Jan 6, 2021 at 15:51
  • $\begingroup$ is there any way by using this "bamPEFragmentSize" i can get smooth curve ? $\endgroup$
    – kcm
    Jan 9, 2021 at 15:32
  • 1
    $\begingroup$ Save the values to a file and replot however you like. $\endgroup$
    – Devon Ryan
    Jan 9, 2021 at 15:48
  • $\begingroup$ okay i will see the option to save the file ....right now i ran this "bamPEFragmentSize -hist fragmentSize.png --maxFragmentLength 100 -b HSC1.bam Mono1.bam LSC1.bam Blast1.bam --numberOfProcessors 10 --samplesLabel HSC Mono LSC Blast " $\endgroup$
    – kcm
    Jan 9, 2021 at 16:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.