Split bam files from ATAC-seq PE data into nucleosome free and nucleosomal regions

Question

[anlo]

Hi guys,

I am looking for some options to split a BAM file based on the insert size with the purpose of separating nucleosome free regions (i.e. regions below approx. 100bp) in ATAC-seq data from larger nucleosomal regions. As far as I can see from discussions at this site, some use deeptools/bamCoverage with specification of --maxFragmentLength to generate a bedgraph/bigwig file that only contains nucleosome free region. I would preferably like to keep the bam file format, so do any of you have some ideas to a commando/tool that do the job of splitting a bam file based on insert sizes.

Thanks, Anne

Answer 1

[ATpoint]

This piece of code can extract fragments of a certain length and write to a new bam file:

#!/bin/bash
samtools view -h $1 | \
  awk -v LEN=$2 '{if ($9 <= LEN && $9 >= -(LEN) && $9 != 0 || $1 ~ /^@/) print $0}' | \
  samtools view -bh -o out.bam -

## Example to get fragments of 147bp or smaller:
./code.sh in.bam 147

For matters of completeness, you may check out the NucleoATAC software, which is designed to infer nucleosome-free and nucleosome-occupied regions within open chromatin domains. Documentation is (at least for me) quiet bulky and uneasy to digest, but maybe it can help you.

Answer 2

[Friederike]

haven't used it myself, but I remember that ATACseqQC did something along those lines -- check out their vignette