broadinstitute · lucidtronix · May 23, 2018 · Sep 17, 2018 · Sep 18, 2018 · Sep 18, 2018
diff --git a/api_tutorials/bqsr_cnn.py b/api_tutorials/bqsr_cnn.py
diff --git a/api_tutorials/howd_we_do.py b/api_tutorials/howd_we_do.py
diff --git a/api_tutorials/recal_utils.py b/api_tutorials/recal_utils.py
@@ -0,0 +1,69 @@
+import numpy as np
+import numpy.ma as ma
+from scipy.stats import entropy
+
+
+def get_next_batch(generator, quality_correction=0.0):
+    batch = next(generator)
+    tensor = batch[0][OQ_TENSOR_NAME]
+    bqsr = batch[0][BQSR_TENSOR_NAME]
+    label = batch[1]
+    pred = model.predict_on_batch(tensor)
+    pred_qscores = -10 * np.log10(
+        pred[:, :, args.labels['BAD_BASE']]) + quality_correction  # +10 only if the tensor is generated with a bias
+
+    orig_qscores = -10 * np.log10(1 - np.max(tensor[:, :, :4], axis=2))
+    annot = tensor[:, 0, (args.input_symbols['pair'], args.input_symbols['mq'])]
+
+    return pred_qscores, orig_qscores, bqsr, label, annot
+
+
+def tensor_to_quality_array(tensor):
+    '''
+    tensor : (batch_size, 151, 7)
+
+    returns : (batch_size, 151)
+    '''
+    return -10 * np.log10(1 - np.max(tensor[:, :, :4], axis=2))
+
+def KL_divergence_metric(y_true, y_pred):
+    '''	KL divergence metrics for Keras - still under construction '''
+
+    # maybe scikit learn
+    predicted_qs = -10*np.log10(y_pred[:,:, args.labels['BAD_BASE']])
+    match_qs = (predicted_qs[:,:,np.newaxis] * y_true)[:,:, args.labels['GOOD_BASE']]
+    mismatch_qs = (predicted_qs[:,:,np.newaxis] * y_true)[:,:, args.labels['BAD_BASE']]
+    match_qs = match_qs[match_qs > 0]
+    mismatch_qs = mismatch_qs[mismatch_qs > 0]
+
+    # bins are half open: 1 will go in the [1,2) bin
+    max_quality=50
+    match_hist, match_bins = np.histogram(np.round(match_qs), bins=max_quality, range = (0,max_quality))
+    mismatch_hist, mismatch_bins = np.histogram(np.round(mismatch_qs), bins=max_quality, range = (0,max_quality))
+
+    # compute the KL divergence KL(match||mismatch) - the order chosen arbitrariliy i.e. could've easily chosen KL(mismatch||match)
+    # mask bins with 0 probability mass because numpy doens't know 0*log(0)=0
+    ma_match_hist = ma.array(match_hist/np.sum(match_hist), mask=match_hist == 0)
+    ma_mismatch_hist = ma.array(mismatch_hist/np.sum(mismatch_hist), mask=match_hist == 0)
+    print(ma_match_hist)
+    print(ma_mismatch_hist)
+    print(entropy(ma_match_hist, ma_mismatch_hist))
+    KL = -ma.sum(ma_match_hist * ma.log(ma_mismatch_hist)) - ma.sum(- ma_match_hist * ma.log(ma_match_hist))
+    return KL
+
+def KL_divergence(match_qs, mismatch_qs):
+    ''' compute the KL divergence between the predicted qualities of bases that match the reference and those that don't
+    match_qs and mismatch_qs are both arrays of qualities, unsorted and unrounded, straight out of the CNN or SAM.
+    greater the KL divergence, the greater the separation between the two distributions
+    '''
+    # bins are half open: 1 will go in the [1,2) bin
+    max_quality=50
+    match_hist, match_bins = np.histogram(np.round(match_qs), bins=max_quality, range = (0,max_quality))
+    mismatch_hist, mismatch_bins = np.histogram(np.round(mismatch_qs), bins=max_quality, range = (0,max_quality))
+
+    # compute the KL divergence KL(match||mismatch) - the order chosen arbitrariliy i.e. could've easily chosen KL(mismatch||match)
+    # mask bins with 0 probability mass because numpy doens't know 0*log(0)=0
+    ma_match_hist = ma.array(match_hist, mask=match_hist == 0)
+    ma_mismatch_hist = ma.array(mismatch_hist, mask=match_hist == 0)
+    KL = -ma.sum(ma_match_hist * ma.log(ma_mismatch_hist)) - ma.sum(- ma_match_hist * ma.log(ma_match_hist))
+    return KL
diff --git a/arguments.py b/arguments.py
@@ -27,6 +27,8 @@ def is_broad_cluster():
 	import keras.backend as K
 
 def parse_args():
+	import sys
+	print(sys.argv)
 	parser = argparse.ArgumentParser()
 
 	# Required mode argument: what would you like to do?
@@ -207,7 +209,7 @@ def weight_path_from_args(args):
 	Arguments:
 		args: puts arguments into the file name skips args in the ignore array
 	'''		
-	save_weight_hd5 =  args.output_dir + args.id + '.hd5'
+	save_weight_hd5 = args.output_dir + args.id + '.hd5'
 	print('save weight path:' , save_weight_hd5)
 	return save_weight_hd5
 

diff --git a/quality/analyze.ipynb b/quality/analyze.ipynb
diff --git a/quality/quick_tests.ipynb b/quality/quick_tests.ipynb
diff --git a/quality/visualize.ipynb b/quality/visualize.ipynb
diff --git a/scripts/hc_vqsr_script.sh b/scripts/hc_vqsr_script.sh
@@ -64,7 +64,7 @@ SCATTER=200
 # REFERENCE=/seq/references/Homo_sapiens_assembly38/v0/Homo_sapiens_assembly38.fasta
 # INTERVAL=/seq/references/Homo_sapiens_assembly38/v0/variant_calling/wgs_calling_regions.v1.interval_list 
 
-# Clinical project g947x NA12878 HG38
+# Clinical project illumina NA12878 HG38
 BAM=/dsde/data/deep/vqsr/bams/NA12878_S1.bam
 BAMOUT=/dsde/data/deep/vqsr/bams/illumina_na12878_s1_bamout.bam
 VCF=/dsde/data/deep/vqsr/vcfs/illumina_na12878_s1.vcf.gz

diff --git a/scripts/write_tensors.sh b/scripts/write_tensors.sh
@@ -308,50 +308,50 @@ SPLIT_INTERVALS=/dsde/data/deep/vqsr/beds/wgs_10m_split_genome.interval_list
 # BED_FILE=/dsde/data/deep/vqsr/beds/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_nosomaticdel_noCENorHET7.bed
 
 # Clinical NA12878 1ug
-SAMPLE_NAME=SM-G9481
-DOWNSAMPLE_SNPS=0.003
-DOWNSAMPLE_INDELS=0.025
-DOWNSAMPLE_NOT_SNPS=0.5
-CHANNEL_ORDER=channels_last
-MODE=write_paired_read_tensors
-DATA_DIR=/dsde/data/deep/vqsr/tensors/g947_balanced_paired_read_channels_last/
-TRAIN_VCF=/dsde/data/deep/vqsr/vcfs/nist_na12878_giab_hg38.vcf.gz
-BAM_FILE=/dsde/data/deep/vqsr/bams/g94781_lod_1ug_na12878_bamout.bam
-NEGATIVE_VCF=/dsde/data/deep/vqsr/vcfs/g94781_lod_1ug_na12878_hc4_merged.vcf.gz
-SPLIT_INTERVALS=/dsde/data/deep/vqsr/beds/wgs_10m_split_genome_hg38.interval_list
-REFERENCE=/seq/references/Homo_sapiens_assembly38/v0/Homo_sapiens_assembly38.fasta
-BED_FILE=/dsde/data/deep/vqsr/beds/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_nosomaticdel_noCENorHET7.bed
-
-# Clinical NA12878 g947m
+# SAMPLE_NAME=SM-G9481
 # DOWNSAMPLE_SNPS=0.003
 # DOWNSAMPLE_INDELS=0.025
 # DOWNSAMPLE_NOT_SNPS=0.5
-# MODE=write_paired_read_tensors
-# SAMPLE_NAME=SM-G947M
 # CHANNEL_ORDER=channels_last
+# MODE=write_paired_read_tensors
 # DATA_DIR=/dsde/data/deep/vqsr/tensors/g947_balanced_paired_read_channels_last/
 # TRAIN_VCF=/dsde/data/deep/vqsr/vcfs/nist_na12878_giab_hg38.vcf.gz
-# BAM_FILE=/dsde/data/deep/vqsr/bams/g947m_o1d2v1_na12878_bamout.bam
-# NEGATIVE_VCF=/dsde/data/deep/vqsr/vcfs/g947m_o1d2v1_na12878_hc4_merged.vcf.gz
+# BAM_FILE=/dsde/data/deep/vqsr/bams/g94781_lod_1ug_na12878_bamout.bam
+# NEGATIVE_VCF=/dsde/data/deep/vqsr/vcfs/g94781_lod_1ug_na12878_hc4_merged.vcf.gz
 # SPLIT_INTERVALS=/dsde/data/deep/vqsr/beds/wgs_10m_split_genome_hg38.interval_list
 # REFERENCE=/seq/references/Homo_sapiens_assembly38/v0/Homo_sapiens_assembly38.fasta
 # BED_FILE=/dsde/data/deep/vqsr/beds/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_nosomaticdel_noCENorHET7.bed
 
-# Clinical NA12878 g947x
+# Clinical NA12878 g947m
 # DOWNSAMPLE_SNPS=0.003
 # DOWNSAMPLE_INDELS=0.025
 # DOWNSAMPLE_NOT_SNPS=0.5
 # MODE=write_paired_read_tensors
-# SAMPLE_NAME=SM-G947X
+# SAMPLE_NAME=SM-G947M
 # CHANNEL_ORDER=channels_last
 # DATA_DIR=/dsde/data/deep/vqsr/tensors/g947_balanced_paired_read_channels_last/
 # TRAIN_VCF=/dsde/data/deep/vqsr/vcfs/nist_na12878_giab_hg38.vcf.gz
-# BAM_FILE=/dsde/data/deep/vqsr/bams/g947x_o2d1v1_na12878_bamout.bam
-# NEGATIVE_VCF=/dsde/data/deep/vqsr/vcfs/g947x_o2d1v1_na12878_cnn_scored.vcf.gz
+# BAM_FILE=/dsde/data/deep/vqsr/bams/g947m_o1d2v1_na12878_bamout.bam
+# NEGATIVE_VCF=/dsde/data/deep/vqsr/vcfs/g947m_o1d2v1_na12878_hc4_merged.vcf.gz
 # SPLIT_INTERVALS=/dsde/data/deep/vqsr/beds/wgs_10m_split_genome_hg38.interval_list
 # REFERENCE=/seq/references/Homo_sapiens_assembly38/v0/Homo_sapiens_assembly38.fasta
 # BED_FILE=/dsde/data/deep/vqsr/beds/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_nosomaticdel_noCENorHET7.bed
 
+# Clinical NA12878 g947x
+DOWNSAMPLE_SNPS=0.005
+DOWNSAMPLE_INDELS=0.025
+DOWNSAMPLE_NOT_SNPS=0.5
+MODE=write_paired_read_tensors
+SAMPLE_NAME=SM-G947X
+DATA_DIR=/dsde/data/deep/vqsr/tensors/g947x_balanced_paired_read_ws128_na12878/
+TRAIN_VCF=/dsde/data/deep/vqsr/vcfs/nist_na12878_giab_hg38.vcf.gz
+BAM_FILE=/dsde/data/deep/vqsr/bams/g947x_o2d1v1_na12878_bamout.bam
+NEGATIVE_VCF=/dsde/data/deep/vqsr/vcfs/g947x_o2d1v1_na12878_cnn_scored.vcf.gz
+SPLIT_INTERVALS=/dsde/data/deep/vqsr/beds/wgs_10m_split_genome_hg38.interval_list
+REFERENCE=/seq/references/Homo_sapiens_assembly38/v0/Homo_sapiens_assembly38.fasta
+BED_FILE=/dsde/data/deep/vqsr/beds/HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_nosomaticdel_noCENorHET7.bed
+
+
 # Clinical NA12878 g947z
 # SAMPLE_NAME=SM-G947Z
 # DOWNSAMPLE_SNPS=0.003