[Bioc-sig-seq] "identical" on PhredQuality objects - I'm confused

Janet Young jayoung at fhcrc.org
Wed Jun 15 03:58:01 CEST 2011


Hi there,

I'm investigating a mistake I made when running BWA (a little embarrassing, but it happens).  I had initially failed to use the -I flag on BWA to tell it to convert our Illumina ASCII-64 qualities appropriately.  Now I'm taking a look to see how much differnce it makes in our downstream analysis.

Anyway, that's not why I'm writing - while delving into this I've found something unexpected about the "identical" function, and I'd love to understand that better if possible.

I used scanBam to read in my two BWA output files - with and without the "-I" option when I ran BWA (same input seqs, everything else the same). As I expected, with the -I option BWA converted the qualities (my lane1read2_bwaWithI_bamscan object), and in the other it didn't (my lane1read2_bwaNotI_bamscan object).  That makes sense to me - here's how the quals look after reading BWA output in with scanBam.

### when I use the -I option with BWA:
> head(lane1read2_bwaWithI_bamscan[[1]][["qual"]])
  A PhredQuality instance of length 6
    width seq
[1]    50 HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHFHHHHFHFF
[2]    36 GFFEEHHFHHHHHHHEFEEEHHHHHFHHHEEFEE at E
[3]    50 ##################################################
[4]    31 HHHHHHHHHHHHHHHHHHHHHHEHHHHFHFH
[5]    24 HHHHHHHHHHHHHHDFHHHHHHHH
[6]    50 FFFFCD=8BGEBHGHHHHEHGFGHHDHEHFHEHHHHGHHHHHHHHHGHHH

### when I don't use the -I option with BWA:
> head(lane1read2_bwaNotI_bamscan[[1]][["qual"]]) 
  A PhredQuality instance of length 6
    width seq
[1]    50 gggggggggggggggggggggggggggggggggggggggggeggggegee
[2]    36 feeddggegggggggdedddgggggegggddedd_d
[3]    50 BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
[4]    31 ggggggggggggggggggggggdggggegeg
[5]    24 ggggggggggggggcegggggggg
[6]    50 eeeebc\Wafdagfggggdgfefggcgdgegdggggfgggggggggfggg


The part I don't understand is when I use identical to test whether those two qual objects are the same, it says they are identical:

> identical( lane1read2_bwaWithI_bamscan[[1]][["qual"]], lane1read2_bwaNotI_bamscan[[1]][["qual"]]) 
[1] TRUE

In one sense they're the same (if you convert from Illumina qualities to normal Phred qualities), but how did identical know to do that?  I might be misunderstanding identical - I thought that EVERYTHING about the object had to be identical.  

I want to understand this as it seems like a dangerous misunderstanding - I use identical a fair amount to make sure I'm not screwing things up.

Really I should probably have converted the "qual" item in my scamBam-generated list to a SolexaQuality object to keep track of the fact they're offset by 64 and not 33, but that's OK. At this point I just want to take a look at what's changed about the BWA output when I use that -I flag.

Any thoughts?

thanks,

Janet



More information about the Bioc-sig-sequencing mailing list