I subset a BAM file by some specific inter-chromosomal regions (e.g. samtools view chrXYZ) resulting in the loss of some mates when the read pairs overlapped the edges of these regions. I want to remove any remaining singleton reads after this but the obvious the standard method based on the sam flags doesn't work.
I can identify the culprits using:
java -Xmx10g -jar picard.jar ValidateSamFile I=in.bam O=out.log
But I can't seem to find a tool that can remove them automatically? Was hoping there was a parameter somewhere in a picard/samtool tool but wasn't able to find it (e.g. IGNORE_MISSING_MATES=false). Rather than having to extract the read ids and doing it manually (these are huge BAMs).
Anyone have any suggestions? Thanks.
•
written
by
benformatics
Since most tools use the SAM flags to identify unpaired reads (whereas here this is not the case as the reads are removed manually)
I decided to do the following:
Here is the command that sorts, filters, converts/re-sorts (I might be using an older version of samtools):
samtools sort -n -O BAM -o input.bam -@ 5 input_sorted.bam; samtools view -h input_sorted.bam | awk -f script.awk - > final.sam; samtools sort -@ 5 -O BAM -o final.bam final.sam; rm final.sam; samtools index -@ 5 final.bam
For the script.awk, most of it was written by ChatGPT... :
#!/usr/bin/awk -f
BEGIN {
prev = ""
prevline = ""
}
{
if ($0 ~ "^@") {
print $0
next
}
if ($1 == prev) {
print prevline
print $0
} else {
prev = $1
prevline = $0
}
}
NOTE: This requires that the read names for paired mates are actually identical (e.g. no \1 \2 etc...)
Log in to answer this question.
With BBTools
actually this doesn't work... it still seems to do this based on the SAM flags. i'm trying repair.sh but it is extremely slow
Thanks for the observation.
How are you using
repair.sh? If you know the read id's of the orphan reads perhaps it would be easier to filter the BAM file to remove those lines?