Questions about bam to fastq

Question

[oghzzang]

Dear Biostar users.

I wanna convert bam to fastq. So I used "samtools bam2fq".

This is my script.

fastq1=`echo $(basename ${input_file%%.gz.*})`
fastq2=`echo ${fastq1/R1/R2}`

samtools bam2fq $input/$input_file > $output_dir/${input_file}_fastq

cat $output_dir/${input_file}_fastq | grep '^@.*/1$' -A 3 > $fastq1
cat $output_dir/${input_file}_fastq | grep '^@.*/2$' -A 3 > $fastq2

Is this right?

---

And I wanna identify whether my bam file includes unmapped reads. so I run samtools flagstat.

This is the result. In this result, can I think that my bam file includes unmapped reads?

40304229 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
40304225 + 0 mapped (100.00% : N/A)
39496046 + 0 paired in sequencing
19799602 + 0 read1
19696444 + 0 read2
39219752 + 0 properly paired (99.30% : N/A)
39486133 + 0 with itself and mate mapped
9909 + 0 singletons (0.03% : N/A)
219437 + 0 with mate mapped to a different chr
219437 + 0 with mate mapped to a different chr (mapQ>=5)

Thank you so much

[GenoMax]

I don't completely understand your script since it appears to be circular. I am going to assume that you are feeling a BAM file to samtools bam2fq.

Considering that the total number of reads is almost identical to mapped (only 4 less) it looks like your BAM must have only mapped reads (or your dataset was almost entirely mapped, rare but possible).

[oghzzang]

Dear genomax. Thanks for your reply.

My script is ..

samtools bam2fq [in.bam] > [out.fastq]

cat [out.fastq] | grep '^@.*/1$' -A 3 > [fastq.1]

cat [out.fastq]| grep '^@.*/2$' -A 3 > [fastq.2]

And as you said, my flagstat result looks like all reads are mapped.

You've been great help to me. Thank you.

Answer 1

[5heikki]

Using the program properly would save a lot of time..

Usage: samtools bam2fq [options...] <in.bam>
Options:
  -0 FILE   write paired reads flagged both or neither READ1 and READ2 to FILE
  -1 FILE   write paired reads flagged READ1 to FILE
  -2 FILE   write paired reads flagged READ2 to FILE
  -f INT    only include reads with all bits set in INT set in FLAG [0]
  -F INT    only include reads with none of the bits set in INT set in FLAG [0]
  -n        don't append /1 and /2 to the read name
  -O        output quality in the OQ tag if present
  -s FILE   write singleton reads to FILE [assume single-end]
  -t        copy RG, BC and QT tags to the FASTQ header line
  -v INT    default quality score if not given in file [1]
      --input-fmt-option OPT[=VAL]
               Specify a single input file format option in the form
               of OPTION or OPTION=VALUE
      --reference FILE
               Reference sequence FASTA FILE [null]

[oghzzang]

Dear 5heikki. Thank you so much.

But I can't find how to convert bam to fastq.1 and fastq.2 at a time using bam2fq. In other tutorials using bam2fq, bam is converted to fastq and separated to fastq1, fastq2. In fact, among samtools bam2fq options, i don't know whether "-1" option mean fastq1 file.

Do you mean "samtools bam2fq -1 [fastq1] -2 [fastq2] [in.bam] ?

I'm very much obliged to your comment.

[hiraksarkar.cs]

Absolutely, although the options are listed, they are not really usable. samtools bam2fq -1 [fastq1] -2 [fastq2] [in.bam] this also did not work.

Answer 2

[GenoMax]

Now you can do samtools fastq -1 R1.fastq.gz -2 R2.fastq.gz in.bam with latest samtools (v.1.11). Be sure to name sort (samtools sort -n) your BAM file before you do the conversion. See samtools manual page for option information.