Question: Sorting my .bam with samtools produces differences in the fastq files generated afterwards
0
gravatar for guillepalou4
17 months ago by
guillepalou40 wrote:

Hello everybody,

Let me give you some background so you can understand my problem.

1- I produced a normal .BAM file.

2- I extracted ONLY the unmapped reads, which are 29M reads. It looks like this (first line):

HX6_24184:8:1205:14783:73264    141     *       0       0       55S22M74S       *       0       0       TCAAGAAGTTTTAGCAGAAGAAATTCCAATGCTTTTATTATATGGAGAAATTGAAAATACAGTTTATAGACCAGAAAAATATGATTATTGGACAACTAGATATGACCATACTAAACTAGATCATCCTAAATTATCATATGTAATAAGACCA AAAAFJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJFFFFJJFJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJAJFJJFJJJFFJJJJJJJJJFFJAJFJJJJJJJ<JJJJFJJJJJJFFJAAFAJFJJA RX:Z:AAACACCCAATGAAAC   QX:Z:-AFFFJJJJFJJJJJJ   BC:Z:CGTATCGG   QT:Z:AAF<AFJJ   XS:f:-138       XC:Z:   AC:Z:   AS:f:-137       XM:A:0  AM:A:0  XT:i:0  BX:Z:AAACACCCAATGAAAC-1 RG:Z:HAP6977146:LibraryNotSpecified:1:unknown_fc:0

3- I wanted to convert the unmapped.bam to paired-end fastq files (R1, R2 and singletons).

To convert the .BAM to .fastq I used samtools fastq:

samtools fastq -T BX -s ./singletons.fastq ./phased_possorted_unmapped_bf0x4_bam.bam -1 phased_possorted_unmapped_bf0x4_R1.fastq -2 phased_possorted_unmapped_bf0x4_R2.fastq

From the output I have the following sizes and number of reads:

  • phased_possorted_unmapped_bf0x4_R1.fastq --> 2.1GB and 7M reads
  • phased_possorted_unmapped_bf0x4_R1.fastq --> 2.4GB and 7M reads
  • singletons.fastq --> 4.9GB and 15M reads

Then I wanted to do exactly the same, but with my .BAM sorted by read name before doing this step. I used samtools sort -n of my .bam:

samtools sort -n ./phased_possorted_unmapped_bf0x4_bam.bam -o ./phased_possorted_unmapped_bf0x4_sorted_bam.bam

Then I used samtools fasq again and these are the results:

  • phased_possorted_unmapped_bf0x4_sorted_R1.fastq --> 3.7GB and 12.5M reads
  • phased_possorted_unmapped_bf0x4_sorted_R1.fastq --> 4.2GB and 12.5M reads
  • sorted_singletons.fastq --> 1.4GB and 4M reads

I don't understand why I have these differences in the number of reads for the fastq files. Sorting should not modify anything and I should have the same number of reads regardless the sorting.

I would appreciate if someone knows what is happening.

sort bam samtools next-gen fastq • 841 views
ADD COMMENTlink modified 17 months ago by Devon Ryan91k • written 17 months ago by guillepalou40

guillepalou4 : Please don't change the tag on this question back to tool. tool tag is generally reserved for posts that are about new tools (not for questions, like this one, about existing tools).

ADD REPLYlink modified 17 months ago • written 17 months ago by genomax70k

Hey, sorry I didn't put it back, I just updated my post because I had an error in it. Thanks for the advice though.

ADD REPLYlink written 17 months ago by guillepalou40

Gah, you rebroke the formatting :(

ADD REPLYlink written 17 months ago by Devon Ryan91k
4
gravatar for Devon Ryan
17 months ago by
Devon Ryan91k
Freiburg, Germany
Devon Ryan91k wrote:

It's not actually documented as far as I can tell, but samtools fastq expects the BAM file to be name sorted. You'll otherwise get odd results, since part of what the tool does is compare read names between successive entries.

ADD COMMENTlink written 17 months ago by Devon Ryan91k

Oh, I didn't know that, thanks! I also checked the sorted vs unsorted bam and they have different number of lines. I wouldn't expect that.

ADD REPLYlink written 17 months ago by guillepalou40
1

If an answer was helpful you should upvote it, if the answer resolved your question you should mark it as accepted. Upvote|Bookmark|Accept

ADD REPLYlink written 17 months ago by WouterDeCoster40k

Do you also know, why do I have different number of characters in the same R1 file for example, after using samtools fastq?

I have 1495694956 number of nucleotides and 1494905216 number of quality scores.

It doesn't make sense.

ADD REPLYlink written 17 months ago by guillepalou40
1

That is strange indeed. Can you try using reformat.sh in=your.bam out1=new_R1.fq out2=new_R2.fq (from BBMap Suite) to see if that fares better?

If you know that there are singletons (missing corresponding pair) in the file then you may need to run repair.sh in1=new_R1.fq in2=out1=new_R2.fq out1=clean_R1.fq out2=clean_R2.fq outs=singleton.fq to re-sync the data.

ADD REPLYlink modified 17 months ago • written 17 months ago by genomax70k

Thanks for the answer. I would prefer using samtools fastq to generate my fastq as It's the only tool I've found I can retrieve an specific TAG from the bam (I need to retrieve a BX: tag which is a barcode).

Also I found this on BBMap:

To discard reads that have mismatching lengths of bases and qualities: reformat.sh in=reads.fq out=fixed.fq tossbrokenreads

Although it may be the solution I would prefer to know why do I have this fastq files corrupted after being generated by samtools fastq...

ADD REPLYlink modified 17 months ago • written 17 months ago by guillepalou40

You're likely miscounting.

ADD REPLYlink written 17 months ago by Devon Ryan91k

It's probably... I used this command:

cat R1.fastq | paste - - - - | cut -f 3 | tr -d '\n' | wc -m

1495694956

cat R1.fastq | paste - - - - | cut -f 5 | tr -d '\n' | wc -m

1494905216

Do you know another command that compares the length of each read to the length of the quality scores? Taking into account that the header of my fastq files are like this:

@HX6_24184:8:2213:9242:58356 BX:Z:ACGCTCTAGGCTGCAA-1

ADD REPLYlink written 17 months ago by guillepalou40
1
zcat foo_R1.fastq.gz | awk '{if(NR%4==1) {name=$1} else if (NR%4==2) {s = length($1)} else if(NR%4==0) { if(s != length($1)) print name}}'

That will print out the names of any reads with mismatches in their sequence and quality lengths.

ADD REPLYlink modified 17 months ago • written 17 months ago by Devon Ryan91k

Thanks a lot Devon. It showed me no results so I guess my files are correct. But I wonder if my command was wrong. I suppose it is!

ADD REPLYlink written 17 months ago by guillepalou40

You command will work properly if and only if cut -f 3 and cut -f 5 can be guaranteed to give you the sequence and quality lines. What I expect is happening is that you sometimes have an extra tab in the header line or something like that, which will cause the field numbers to change and the results to get thrown off.

ADD REPLYlink written 17 months ago by Devon Ryan91k

I see... Anyway I checked and both gives me the sequence and the quality lines respectively. It's strange, I will try again similarly today. Thanks for the help :).

ADD REPLYlink written 17 months ago by guillepalou40
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: 678 users visited in the last hour