Metrics: variation


In [1]:
%run ../../shared_setup.ipynb


docker image cggh/biipy:v1.6.0

In [2]:
# load PASS variants for all three crosses
callsets = load_callsets(COMBINED_CALLSET_FN_TEMPLATE, 
                         variant_filter='FILTER_PASS')


2016-03-09 15:02:43.457160 :: loading /data/plasmodium/pfalciparum/pf-crosses/data/public/20141022/3d7_hb3.combined.final.npz
2016-03-09 15:02:43.816991 :: filter variants: excluding 157 (0.4%) retaining 42087 (99.6%) of 42244 variants
2016-03-09 15:02:43.834907 :: loading /data/plasmodium/pfalciparum/pf-crosses/data/public/20141022/hb3_dd2.combined.final.npz
2016-03-09 15:02:44.252911 :: filter variants: excluding 450 (1.2%) retaining 36461 (98.8%) of 36911 variants
2016-03-09 15:02:44.272729 :: loading /data/plasmodium/pfalciparum/pf-crosses/data/public/20141022/7g8_gb4.combined.final.npz
2016-03-09 15:02:44.683451 :: filter variants: excluding 304 (0.9%) retaining 34471 (99.1%) of 34775 variants

No. of independent progeny clones


In [3]:
for cross in CROSSES:
    samples = callsets[cross]['calldata'].dtype.names
    progeny = samples[2:]
    progeny_clones = set([p.split('/')[0] for p in progeny])
    print(cross, len(progeny_clones))


3d7_hb3 15
hb3_dd2 35
7g8_gb4 28

Coverage


In [4]:
tbl_samples = (etl
    .fromtsv(os.path.join(PUBLIC_DIR, 'samples.txt'))
    .convert('coverage', lambda v: int(v[:-1]))
)
tbl_samples


Out[4]:
0|cross 1|clone 2|sample 3|run 4|instrument 5|coverage
3d7_hb3 3D7 PG0051-C ERR019061 Illumina Genome Analyzer II 122
3d7_hb3 C01 PG0065-C ERR019064 Illumina Genome Analyzer II 163
3d7_hb3 C01 PG0062-C ERR019070 Illumina Genome Analyzer II 108
3d7_hb3 C02 PG0055-C ERR019066 Illumina Genome Analyzer II 102
3d7_hb3 C02 PG0053-C ERR019067 Illumina Genome Analyzer II 73

...


In [5]:
tbl_samples.valuecounts('cross')


Out[5]:
0|cross 1|count 2|frequency
7g8_gb4 40 0.40816326530612246
hb3_dd2 37 0.37755102040816324
3d7_hb3 21 0.21428571428571427

In [6]:
df_samples = tbl_samples.todataframe()
df_samples.groupby('cross').coverage.median()


Out[6]:
cross
3d7_hb3    102.0
7g8_gb4    106.5
hb3_dd2    110.0
Name: coverage, dtype: float64

In [7]:
df_samples.groupby('cross').coverage.min()


Out[7]:
cross
3d7_hb3    41
7g8_gb4    55
hb3_dd2    22
Name: coverage, dtype: int64

In [8]:
df_samples.groupby('cross').coverage.max()


Out[8]:
cross
3d7_hb3    173
7g8_gb4    250
hb3_dd2    637
Name: coverage, dtype: int64

In [9]:
tbl_samples.aggregate('cross', [('median', 'coverage', lambda g: np.median(list(g))),
                                ('min', 'coverage', min),
                                ('max', 'coverage', max)])


Out[9]:
0|cross 1|median 2|min 3|max
3d7_hb3 102.0 41 173
7g8_gb4 106.5 55 250
hb3_dd2 110.0 22 637

Count SNPs and INDELs


In [17]:
def count_variants(query):
    def f(row):
        callset = filter_variants(callsets[row.cross], query=query)
        return callset['variants'].size
    return f
        

tbl_variation = (etl
    .wrap([['cross']] + [[cross] for cross in CROSSES])
    .addfield('n_snps', count_variants('is_snp'))
    .addfield('n_indels', count_variants('~is_snp'))
    .addfield('n_snps_multiallelic', count_variants('is_snp & (num_alleles > 2)'))
    .addfield('n_indels_multiallelic', count_variants('~is_snp & (num_alleles > 2)'))
    .addfield('n_snps_multiallelic_gatk', count_variants('is_snp & (num_alleles > 2) & ((set == b"GATK") | (set == b"Intersection"))'))
    .addfield('n_indels_multiallelic_gatk', count_variants('~is_snp & (num_alleles > 2) & ((set == b"GATK") | (set == b"Intersection"))'))
    .addfield('n_snps_coding', count_variants('is_snp & (CDSAnnotationID != b".")'))
    .addfield('n_snps_noncoding', count_variants('is_snp & (CDSAnnotationID == b".")'))
    .addfield('n_indels_coding', count_variants('~is_snp & (CDSAnnotationID != b".")'))
    .addfield('n_indels_noncoding', count_variants('~is_snp & (CDSAnnotationID == b".")'))
    .addfield('ratio_snp_indel_coding', lambda row: row.n_snps_coding / row.n_indels_coding)
    .addfield('ratio_snp_indel_noncoding', lambda row: row.n_snps_noncoding / row.n_indels_noncoding)
    .melt(key='cross')
    .recast(variablefield='cross', valuefield='value')
)
tbl_variation.displayall()


2016-03-09 15:09:32.794681 :: filter variants: excluding 26699 (63.4%) retaining 15388 (36.6%) of 42087 variants
2016-03-09 15:09:32.804509 :: filter variants: excluding 15388 (36.6%) retaining 26699 (63.4%) of 42087 variants
2016-03-09 15:09:32.816272 :: filter variants: excluding 42075 (100.0%) retaining 12 (0.0%) of 42087 variants
2016-03-09 15:09:32.817924 :: filter variants: excluding 41708 (99.1%) retaining 379 (0.9%) of 42087 variants
2016-03-09 15:09:32.825786 :: filter variants: excluding 42075 (100.0%) retaining 12 (0.0%) of 42087 variants
2016-03-09 15:09:32.832307 :: filter variants: excluding 41708 (99.1%) retaining 379 (0.9%) of 42087 variants
2016-03-09 15:09:32.835129 :: filter variants: excluding 33219 (78.9%) retaining 8868 (21.1%) of 42087 variants
2016-03-09 15:09:32.843661 :: filter variants: excluding 35567 (84.5%) retaining 6520 (15.5%) of 42087 variants
2016-03-09 15:09:32.848935 :: filter variants: excluding 37981 (90.2%) retaining 4106 (9.8%) of 42087 variants
2016-03-09 15:09:32.853397 :: filter variants: excluding 19494 (46.3%) retaining 22593 (53.7%) of 42087 variants
2016-03-09 15:09:32.865127 :: filter variants: excluding 21576 (59.2%) retaining 14885 (40.8%) of 36461 variants
2016-03-09 15:09:32.876137 :: filter variants: excluding 14885 (40.8%) retaining 21576 (59.2%) of 36461 variants
2016-03-09 15:09:32.890091 :: filter variants: excluding 36454 (100.0%) retaining 7 (0.0%) of 36461 variants
2016-03-09 15:09:32.891714 :: filter variants: excluding 32627 (89.5%) retaining 3834 (10.5%) of 36461 variants
2016-03-09 15:09:32.898444 :: filter variants: excluding 36454 (100.0%) retaining 7 (0.0%) of 36461 variants
2016-03-09 15:09:32.902711 :: filter variants: excluding 32627 (89.5%) retaining 3834 (10.5%) of 36461 variants
2016-03-09 15:09:32.908679 :: filter variants: excluding 27853 (76.4%) retaining 8608 (23.6%) of 36461 variants
2016-03-09 15:09:32.915458 :: filter variants: excluding 30184 (82.8%) retaining 6277 (17.2%) of 36461 variants
2016-03-09 15:09:32.921221 :: filter variants: excluding 32782 (89.9%) retaining 3679 (10.1%) of 36461 variants
2016-03-09 15:09:32.926253 :: filter variants: excluding 18564 (50.9%) retaining 17897 (49.1%) of 36461 variants
2016-03-09 15:09:32.937502 :: filter variants: excluding 20079 (58.2%) retaining 14392 (41.8%) of 34471 variants
2016-03-09 15:09:32.948193 :: filter variants: excluding 14392 (41.8%) retaining 20079 (58.2%) of 34471 variants
2016-03-09 15:09:32.962656 :: filter variants: excluding 34468 (100.0%) retaining 3 (0.0%) of 34471 variants
2016-03-09 15:09:32.964265 :: filter variants: excluding 30739 (89.2%) retaining 3732 (10.8%) of 34471 variants
2016-03-09 15:09:32.970926 :: filter variants: excluding 34468 (100.0%) retaining 3 (0.0%) of 34471 variants
2016-03-09 15:09:32.974892 :: filter variants: excluding 30739 (89.2%) retaining 3732 (10.8%) of 34471 variants
2016-03-09 15:09:32.979990 :: filter variants: excluding 26266 (76.2%) retaining 8205 (23.8%) of 34471 variants
2016-03-09 15:09:32.987572 :: filter variants: excluding 28284 (82.1%) retaining 6187 (17.9%) of 34471 variants
2016-03-09 15:09:32.995863 :: filter variants: excluding 30740 (89.2%) retaining 3731 (10.8%) of 34471 variants
2016-03-09 15:09:33.000212 :: filter variants: excluding 18123 (52.6%) retaining 16348 (47.4%) of 34471 variants
2016-03-09 15:09:33.013067 :: filter variants: excluding 26699 (63.4%) retaining 15388 (36.6%) of 42087 variants
2016-03-09 15:09:33.022449 :: filter variants: excluding 15388 (36.6%) retaining 26699 (63.4%) of 42087 variants
2016-03-09 15:09:33.034443 :: filter variants: excluding 42075 (100.0%) retaining 12 (0.0%) of 42087 variants
2016-03-09 15:09:33.036097 :: filter variants: excluding 41708 (99.1%) retaining 379 (0.9%) of 42087 variants
2016-03-09 15:09:33.040794 :: filter variants: excluding 42075 (100.0%) retaining 12 (0.0%) of 42087 variants
2016-03-09 15:09:33.045353 :: filter variants: excluding 41708 (99.1%) retaining 379 (0.9%) of 42087 variants
2016-03-09 15:09:33.048946 :: filter variants: excluding 33219 (78.9%) retaining 8868 (21.1%) of 42087 variants
2016-03-09 15:09:33.057261 :: filter variants: excluding 35567 (84.5%) retaining 6520 (15.5%) of 42087 variants
2016-03-09 15:09:33.062726 :: filter variants: excluding 37981 (90.2%) retaining 4106 (9.8%) of 42087 variants
2016-03-09 15:09:33.068621 :: filter variants: excluding 19494 (46.3%) retaining 22593 (53.7%) of 42087 variants
2016-03-09 15:09:33.080517 :: filter variants: excluding 21576 (59.2%) retaining 14885 (40.8%) of 36461 variants
2016-03-09 15:09:33.090563 :: filter variants: excluding 14885 (40.8%) retaining 21576 (59.2%) of 36461 variants
2016-03-09 15:09:33.105211 :: filter variants: excluding 36454 (100.0%) retaining 7 (0.0%) of 36461 variants
2016-03-09 15:09:33.107869 :: filter variants: excluding 32627 (89.5%) retaining 3834 (10.5%) of 36461 variants
2016-03-09 15:09:33.115376 :: filter variants: excluding 36454 (100.0%) retaining 7 (0.0%) of 36461 variants
2016-03-09 15:09:33.121624 :: filter variants: excluding 32627 (89.5%) retaining 3834 (10.5%) of 36461 variants
2016-03-09 15:09:33.127136 :: filter variants: excluding 27853 (76.4%) retaining 8608 (23.6%) of 36461 variants
2016-03-09 15:09:33.133621 :: filter variants: excluding 30184 (82.8%) retaining 6277 (17.2%) of 36461 variants
2016-03-09 15:09:33.140690 :: filter variants: excluding 32782 (89.9%) retaining 3679 (10.1%) of 36461 variants
2016-03-09 15:09:33.145140 :: filter variants: excluding 18564 (50.9%) retaining 17897 (49.1%) of 36461 variants
2016-03-09 15:09:33.157373 :: filter variants: excluding 20079 (58.2%) retaining 14392 (41.8%) of 34471 variants
2016-03-09 15:09:33.167518 :: filter variants: excluding 14392 (41.8%) retaining 20079 (58.2%) of 34471 variants
2016-03-09 15:09:33.180130 :: filter variants: excluding 34468 (100.0%) retaining 3 (0.0%) of 34471 variants
2016-03-09 15:09:33.181623 :: filter variants: excluding 30739 (89.2%) retaining 3732 (10.8%) of 34471 variants
2016-03-09 15:09:33.189432 :: filter variants: excluding 34468 (100.0%) retaining 3 (0.0%) of 34471 variants
2016-03-09 15:09:33.192408 :: filter variants: excluding 30739 (89.2%) retaining 3732 (10.8%) of 34471 variants
2016-03-09 15:09:33.199867 :: filter variants: excluding 26266 (76.2%) retaining 8205 (23.8%) of 34471 variants
2016-03-09 15:09:33.208768 :: filter variants: excluding 28284 (82.1%) retaining 6187 (17.9%) of 34471 variants
2016-03-09 15:09:33.213852 :: filter variants: excluding 30740 (89.2%) retaining 3731 (10.8%) of 34471 variants
2016-03-09 15:09:33.220260 :: filter variants: excluding 18123 (52.6%) retaining 16348 (47.4%) of 34471 variants
0|variable 1|3d7_hb3 2|7g8_gb4 3|hb3_dd2
n_indels 26699 20079 21576
n_indels_coding 4106 3731 3679
n_indels_multiallelic 379 3732 3834
n_indels_multiallelic_gatk 379 3732 3834
n_indels_noncoding 22593 16348 17897
n_snps 15388 14392 14885
n_snps_coding 8868 8205 8608
n_snps_multiallelic 12 3 7
n_snps_multiallelic_gatk 12 3 7
n_snps_noncoding 6520 6187 6277
ratio_snp_indel_coding 2.159766195811008 2.199142321093541 2.3397662408263113
ratio_snp_indel_noncoding 0.28858495994334527 0.37845608025446537 0.350729172487009

In [15]:
callsets['3d7_hb3']['variants']['set']


Out[15]:
array([b'GATK', b'Intersection', b'Intersection', ..., b'GATK',
       b'Intersection', b'Intersection'], 
      dtype='|S40')