How to extract specific region from multiple bam files and merge the outputs into single bam?
I have more than 50 bam files and wanted to extract the reads from a region from all those bam files and merge the outputs into a single bam.
I'm using samtools for this:
samtools view -b sample1.bam "Chr10:19000-46500" > sample1_intregion.bam
I tried doing like below for all bam files:
#!/bin/bash
for f in *.bam
do
samtools view -b $f "Chr10:19000-46500" > $f_intregion.bam
done
it is not working and doesn't have the output named above. And how can I include the command into above loop for also merging all the output bam files into single bam?
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
paired-end reads, idxstats, and samtools and grep
written by Bogdan
-
Combine multiple samtools view
written by rah •
-
Splitting BAM file from paired end into forward and reverse strands?
written by Susmita Mandal
-
Consensus sequence from multiple BAM files on a single position
written by Federico Giorgi
-
Chimeric read insert length in SAM $9
written by Yijun Tian •
-
Writing a bash loop to convert sams to bams
written by raf.marcondes •
-
Error when split WES data into sub-files
written by Laven9 •
-
Larger number of reads after filtering 'chrM'??
written by cpc5 •
-
SAMTOOLS flagstat outputs
written by jt358 •
-
Samtools: append output
written by bongbang •
What is not working? Do you get an error?
Instead of
use
You could merge the extracted BAM files after all extractions are complete using
samtools merge.Doing
$f_intregion.bamalso will make the names odd likefile.bam_integration.bam.yes it says
[main_samview] random alignment retrieval only works for indexed BAM or CRAM files.NOTE: I had assumed that the pipe in block 1 below would work (without testing it at the time I wrote this up). But it does not work.
samtools viewextraction requires an index file and that can't be done on fly by the first step writing to STDOUT. So this operation will need to be done in two steps as far as I can tell (if someone has a better answer please post).So index the files first (in addition to sorting them) using
samtools sort --write-index.Following should work. Ideally the sorting should be done once outside this loop.
In the end
samtools mergeall the extracted files.I already have index files. I did with
samtools indexI tried this but have
sort: unrecognized option '--write-index'You must be using an old version of
samtools. What do you get when you dosamtools --version?1.6 is the version I'm using
Please upgrade.
v.1.6is ancient and was released in 2017. We are currently on1.18. (this numbering goes like 1.8,1.9,1.10,1.11 etc),I got the new version and extracted the tar.gz file and added the path to bashrc and did source it but still it says like below:
Did you compile the program? You will need to add the directory which contains the new
samtoolsto your$PATHby doingexport $PATH=/new_samtools_dir:$PATH.Be sure to confirm the
samtools --versionshows the new version (i.e. it is one that is seen first before your default1.6).Okay I installed the new version and it is working now. Thanq !!
But if I use the
|in the loop it is not wokring. So, I removed itThat's not how "thank you" is spelled in a professional setting. Please write professionally.
the
--write-indexis a quite new command, and due to the various library dependencies, people don't get the latest samtools version when installing with condaso the recommendation would be best to avoid that flag for a few years, maybe,
especially since one can index separately with ease
samtools index align.bamthat will always workalso for the OP, the best would be to sort their existing BAM files once and for all because sooner or later, they need the sorted BAM file, so resorting on the fly each time is not worth it.