How could I extract reads from bam files by their read groups @RG ?, I've got one file containing all reads for 5 samples, the SM is appears NONE, So, I want to extract each sample reads by the @RG tags.
-
$\begingroup$ Hi and welcome, please have a think about targeting your question better. What are the genome sizes and machine? What did you try command wise? $\endgroup$– M__ ♦Oct 9, 2019 at 23:34
-
$\begingroup$ Okey thank u for answering. $\endgroup$– MalokiOct 9, 2019 at 23:39
-
$\begingroup$ So, I use ion torrent to sequence 3 bactrial whole genome and 2 samples of fungi. I used 5 codebars but when I launched the plan I Forgot to check if I put the names of each samples. When I've got the ubam file I found only one file containing all the reads. So I Want To know if I could get each sample by its @RG or ID. Knowing that the SM tag appears "None". $\endgroup$– MalokiOct 9, 2019 at 23:47
2 Answers
Typically I use samtools for operations like this. Specifically I use samtools view
with either -r
or -R
flag depending on the use case.
-r STR
Output alignments in read group STR [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.
-R FILE
Output alignments in read groups listed in FILE [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.
This sounds more like a job for samtools split
if you want to split out all the read groups into separate bams in one go.
http://www.htslib.org/doc/samtools-split.html
samtools split [options] merged.sam|merged.bam|merged.cram
Splits a file by read group.
Options:
-u FILE1
Put reads with no RG tag or an unrecognised RG tag into FILE1
-u FILE1:FILE2
As above, but assigns an RG tag as given in the header of FILE2
-f STRING
Output filename format string (see below) ["%*_%#.%."]
-v
Verbose output
Alternatively if you do want to proceed one read group at a time, in addition to samtools view
, you can also do this with GATK using the PrintReads tool with the ReadGroupReadFilter as documented here:
https://gatk.broadinstitute.org/hc/en-us/articles/360041415172-ReadGroupReadFilter
The command line would look like this:
gatk PrintReads \
-I input.bam \
--read-filter ReadGroupReadFilter \
--keep-read-group <RG:readgroup> \
-O output.bam
There's also a filter called ReadGroupBlackListReadFilter that allows you to exclude a list of readgroups.
https://gatk.broadinstitute.org/hc/en-us/articles/360041851151-ReadGroupBlackListReadFilter
The PrintReads tool doc is here:
https://gatk.broadinstitute.org/hc/en-us/articles/360041849911-PrintReads