notebook.community

Edit and run

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')

Content source: wtchg-kwiatkowski/pfx-paper-2015

Similar notebooks:

notebook.community | gallery | about