Question: Pipe BBMap reformat output into SAMtools
0
gravatar for bruce.moran
19 months ago by
bruce.moran510
Ireland
bruce.moran510 wrote:

Hi all,

I am wondering if I can pipe output from reformat.sh in BBMap into SAMtools view. My goal is to take a set of fastq and return an unaligned BAM headed with sequence dictionary (this is for use in GATK4). Using some GATK test data this part works (of course!):

bbduk.sh in1=6484_snippet_1.fastq in2=6484_snippet_2.fastq trimq=30 out=stdout.fq | reformat.sh in=stdin.fq out=stdout.bam

I suppose the question then is: can I force reformat.sh to produce a BAM without specifying ".bam" ending? Or, what form does 'stdout.bam' take and how can SAMtools recognise it? Or, can reformat.sh do the reheader itself?

All help much appreciated,

Bruce.

pipe samtools bbmap • 723 views
ADD COMMENTlink modified 19 months ago • written 19 months ago by bruce.moran510
1

See Piping and Screen Output at https://github.com/BioInfoTools/BBMap/blob/master/docs/UsageGuide.txt#L18

ADD REPLYlink written 19 months ago by Santosh Anand4.6k

OK, I had read that, just tried a few things and it seems you have to stream through samtools view:

reformat.sh in=6484_snippet_1.fastq out=stdout.bam | samtools reheader 6484_snippet_header.sam - > 6484_snippet.rh.bam

returns [W::hts_close] EOF marker is absent. The input is probably truncated.

reformat.sh in=6484_snippet_1.fastq out=stdout.bam | samtools view -hb - | reheader 6484_snippet_header.sam - > 6484_snippet.rh.bam

returns a reheaded BAM! Thanks for the response and making me think a bit more.

ADD REPLYlink modified 19 months ago • written 19 months ago by bruce.moran510
3
gravatar for bruce.moran
19 months ago by
bruce.moran510
Ireland
bruce.moran510 wrote:

To answer my own Q:

bbduk.sh in1=6484_snippet_1.fastq \
          in2=6484_snippet_2.fastq \
          trimq=30 out=stdout.fq \
| reformat.sh in=stdin.fq \
               out=stdout.bam \
| samtools view -hb - | samtools reheader 6484_snippet_header.sam - > 6484_snippet.rh.bam
ADD COMMENTlink modified 19 months ago • written 19 months ago by bruce.moran510
1

Incidentally, you can skip the reformat step - bbduk.sh will allow "out=stdout.bam" also :)

ADD REPLYlink written 19 months ago by Brian Bushnell16k

Cool, I had to use addslash and underscore hence the reformat.sh, but good to know.

ADD REPLYlink written 19 months ago by bruce.moran510

Ah, I see. Well, alternately, reformat.sh supports quality-trimming, so I guess it's bbduk that was unnecessary (unless you were doing a kmer-based operation).

ADD REPLYlink written 19 months ago by Brian Bushnell16k

Thought about that also, but yes, using kmers.

ADD REPLYlink written 19 months ago by bruce.moran510

You have learned everything I have to teach,

ADD REPLYlink written 19 months ago by Brian Bushnell16k
1

My RTFM is strong (also helps that there is a good FM)

ADD REPLYlink written 19 months ago by bruce.moran510
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: 1346 users visited in the last hour