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
2 answers
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.
Log in to answer this question.