Question: Custom extract reads from a sam/bam file into a valid sam/bam
0
gravatar for shwethacm
19 months ago by
shwethacm200
Seattle, WA
shwethacm200 wrote:

Hi there,

I am trying to "brute-force" split a bam file using grep. The reason why I am doing this is because I am trying to split by a custom tag that looks like these (bc would be the tag):
bc:B:S,1,1
bc:B:S,0,0
bc:B:S,2,2

(( But I cannot used bamtools because "bc" is not a recognized tag )). My file is small-sih (<3 GB) so I used grep to make my split files and appended the header to each sub-file.

But something is wrong because I cannot index my new sub-files. I ran ValidateSamFile from picard, and it is pointing to a missing tag in the header ("Error parsing SAM header. @RG line missing SM tag.") but not the one that I am interested in.

Does anyone know what I am doing wrong ?

alignment next-gen genome • 695 views
ADD COMMENTlink modified 19 months ago • written 19 months ago by shwethacm200
2
gravatar for Pierre Lindenbaum
19 months ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum122k wrote:

ran ValidateSamFile from picard, and it is pointing to a missing tag in the header ("Error parsing SAM header. @RG line missing SM tag."

that because you're ignoring the bam header

use

samtools view -h in.bam | grep  -E '(^@|bc\:B\:S,1,1)' | samtools view -S -b -o out.bam -

EDIT: when using ValidateSamFile you must also specify what should(not) be checked: http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

e.g:

IGNORE=RECORD_MISSING_READ_GROUP
ADD COMMENTlink modified 19 months ago • written 19 months ago by Pierre Lindenbaum122k
0
gravatar for shwethacm
19 months ago by
shwethacm200
Seattle, WA
shwethacm200 wrote:

Gotcha!

I realized that the mistake I was making was in appending the header separately in the end, which resulted in an ad-hoc mush of a samfile. Your method is so much more elegant!

Thanks, Pierre!

ADD COMMENTlink written 19 months ago by shwethacm200
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1034 users visited in the last hour