0
$\begingroup$

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.

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

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.

$\endgroup$
5
  • $\begingroup$ "bamPEFragmentSize" will see this $\endgroup$
    – kcm
    Jan 6 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 at 15:51
  • $\begingroup$ is there any way by using this "bamPEFragmentSize" i can get smooth curve ? $\endgroup$
    – kcm
    Jan 9 at 15:32
  • 1
    $\begingroup$ Save the values to a file and replot however you like. $\endgroup$
    – Devon Ryan
    Jan 9 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 at 16:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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