GENOMEWIDE.READ.COUNT does not equal the sum of REGION.READ.COUNT

[boxiangliu]

Hello Bryce and Graham,

I notice that the sum of REGION.READ.COUNT is less than GENOMEWIDE.READ.COUNT. Shouldn't they be equal? If not, why? Sorry if I am missing something obvious!

From reading extract_haplotype_read_counts.py, GENOMEWIDE.READ.COUNT is obtained directly from the h5 file, so naively I would assume that some of the reads does not overlap linked regions, therefore GENOMEWIDE.READ.COUNT is greater than the sum of REGION.READ.COUNT.

More importantly, if these two quantities are not equal, how should I calculate GENOMEWIDE.READ.COUNT after permuting REGION.READ.COUNTS (for test calibration)? Thanks!

[gmcvicker]

Hi Bosh,

You are correct that the GENOME.WIDE.READ.COUNT is the total number of genome-wide filtered mapped reads, not simply the number of reads in the linked regions being tested. This is intentional, and the value is used as a scaling coefficient for each individual's read counts.

The way we do the shuffling is to shuffle the genotypes across individuals while preserving the counts. This way the region read counts are kept appropriately paired with the total read counts. The combined_test.py script has a --shuffle option that will do this for you.

Graham