Biostar Beta. Not for public use.
Questions about bam to fastq
0
Entering edit mode
13 months ago
oghzzang • 40

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**

This is 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.

ADD COMMENTlink
0
Entering edit mode

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).

ADD REPLYlink
0
Entering edit mode

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.

ADD REPLYlink
0
Entering edit mode
13 months ago
5heikki 8.4k
Finland

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]
ADD COMMENTlink
0
Entering edit mode

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.

ADD REPLYlink

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.1