Question: Pipe BBMap reformat output into SAMtools
0
gravatar for bruce.moran
2.1 years ago by
bruce.moran620
Ireland
bruce.moran620 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 • 870 views
ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by bruce.moran620
1

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

ADD REPLYlink written 2.1 years ago by Santosh Anand4.9k

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 2.1 years ago • written 2.1 years ago by bruce.moran620
3
gravatar for bruce.moran
2.1 years ago by
bruce.moran620
Ireland
bruce.moran620 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 2.1 years ago • written 2.1 years ago by bruce.moran620
1

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

ADD REPLYlink written 2.1 years ago by Brian Bushnell16k

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

ADD REPLYlink written 2.1 years ago by bruce.moran620

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 2.1 years ago by Brian Bushnell16k

Thought about that also, but yes, using kmers.

ADD REPLYlink written 2.1 years ago by bruce.moran620

You have learned everything I have to teach,

ADD REPLYlink written 2.1 years ago by Brian Bushnell16k
1

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

ADD REPLYlink written 2.1 years ago by bruce.moran620
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: 1307 users visited in the last hour