[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