Entering edit mode
5.4 years ago
Yijun Tian
▴
20
Hello everyone,
I am recently doing a script to dissect chimeric alignment from BAM file, here is my code:
#!/bin/bash
samtools=`which samtools`
bn=$(basename $1 .bam)
echo processing $bn
$samtools view -@ 20 -h -b -f97 -F22 $1 > $bn.flag97.bam
$samtools view -@ 20 -h -b -f145 -F42 $1 > $bn.flag145.bam
$samtools view -@ 20 -h -b -f161 -F22 $1 > $bn.flag161.bam
$samtools view -@ 20 -h -b -f81 -F42 $1 > $bn.flag81.bam
echo merging tmp BAM files
$samtools merge $bn.uns.pp.bam $bn.flag97.bam $bn.flag145.bam $bn.flag161.bam $bn.flag81.bam
echo sorting final BAM file
$samtools view -H $bn.uns.pp.bam -o SamHeader
$samtools view $bn.uns.pp.bam |awk '$7=="=" && $9^2>150^2 {print $0}' > SamChimeric
cat SamHeader SamChimeric |$samtools view -bS > $bn.Chimeric.bam
$samtools sort $bn.Chimeric.bam > $bn.Chimeric.sorted.bam
I have a question about the insert length. Since I was looking for some alignment located in circular DNA, so I assume the insert lengths for such read pair are incalculable. But seeing from the BAM file, the insert length are all given for each alignment:
$ samtools view s4_E.PE.sorted.Chimeric.sorted.bam | cut -f 9 | wc -l
454883
$ samtools view s4_E.PE.sorted.Chimeric.sorted.bam | cut -f 9 | grep -v '\*' | wc -l
454883
How does aligner give the insert size? (I'm using bowtie2)