# ATAC seq : Plotting insert sizes

library(GenomicAlignments)

flag = scanBamFlag(isPaired = TRUE, isProperPair = TRUE), what = c("qname",
"mapq", "isize"), which = GRanges("chr20", IRanges(1, 63025520))))

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.

• 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 Jan 6 at 11:09
• 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 //
– kcm
Jan 6 at 11:24
• @StupidWolf for example in this paper nature.com/articles/s41597-019-0071-0#Sec1 Fig2(b).
– kcm
Jan 6 at 11:26

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.

• "bamPEFragmentSize" will see this
– kcm
Jan 6 at 14:59
• " you're seriously reinventing the wheel." yes since I didn;t read about your tool "bamPEFragmentSize" I used it , Works absolutely fine
– kcm
Jan 6 at 15:51
• is there any way by using this "bamPEFragmentSize" i can get smooth curve ?
– kcm
Jan 9 at 15:32
• Save the values to a file and replot however you like. Jan 9 at 15:48
• 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 "
– kcm
Jan 9 at 16:20