From 6ac2ae0cf1ba99429a7ef496ca49d62b0db1440d Mon Sep 17 00:00:00 2001 From: arpanda Date: Sun, 12 Nov 2023 10:17:48 -0600 Subject: [PATCH] cnvpytor track: bug fixes and use baf_i1 signal --- js/cnvpytor/HDF5IndexedReader.js | 22 +++++++++------ js/cnvpytor/cnvpytorVCF.js | 47 ++++++++++++++++++-------------- 2 files changed, 40 insertions(+), 29 deletions(-) diff --git a/js/cnvpytor/HDF5IndexedReader.js b/js/cnvpytor/HDF5IndexedReader.js index f028105f2..bab77d30e 100644 --- a/js/cnvpytor/HDF5IndexedReader.js +++ b/js/cnvpytor/HDF5IndexedReader.js @@ -17,6 +17,7 @@ class SignalNames{ 'gc_RD': `his_rd_p_${this.chrom}_${this.signal_bin_size}_GC`, 'gc_partition' : `his_rd_p_${this.chrom}_${this.signal_bin_size}_partition_GC_merge`, 'baf': `snp_likelihood_${this.chrom}_${this.signal_bin_size}_mask`, + 'baf_i1': `snp_i1_${this.chrom}_${this.signal_bin_size}_mask`, 'Mosaic_segments' : `his_rd_p_${this.chrom}_${this.signal_bin_size}_partition_GC_mosaic_segments_2d`, 'Mosaic_calls': `his_rd_p_${this.chrom}_${this.signal_bin_size}_partition_GC_mosaic_call_2d` } @@ -49,8 +50,11 @@ class HDF5Reader { return this.h5_obj } + /** + * + * @returns - a list of keys of the pytor file + */ async get_keys(){ - /* returns a list of keys of the pytor file*/ let h5_obj = await this.fetch(); return h5_obj.keys } @@ -58,10 +62,9 @@ class HDF5Reader { async get_rd_signal(bin_size=this.bin_size){ let h5_obj = await this.fetch(); - // let h5_obj_keys = h5_obj.keys; + this.h5_obj = h5_obj this.pytor_keys = h5_obj.keys - // console.log("keys", this.pytor_keys ) // get available bin sizes let signal_bin = new ParseSignals(this.pytor_keys); @@ -119,15 +122,17 @@ class HDF5Reader { // baf likelihood // let signal_baf_1 = `snp_likelihood_${chrom}_${bin_size}_mask` - let signal_baf_1 = signal_name_obj.signals['baf'] - let chr_wig_bafs = await this.get_baf_signals(chrom, bin_size, signal_baf_1) - + // let signal_baf_1 = signal_name_obj.signals['baf'] + // let chr_wig_bafs = await this.get_baf_signals(chrom, bin_size, signal_baf_1) + // let signal_baf_1 = `snp_i1_${chrom}_${bin_size}_mask` - // let chr_wig_bafs = await this.get_baf_signals_v2(h5_obj, h5_obj_keys, chrom, bin_size, signal_baf_1) + let signal_baf_1 = signal_name_obj.signals['baf_i1'] + let chr_wig_bafs = await this.get_baf_signals_v2(chrom, bin_size, signal_baf_1) + wigFeatures_baf1 = wigFeatures_baf1.concat(chr_wig_bafs[0]) wigFeatures_baf2 = wigFeatures_baf2.concat(chr_wig_bafs[1]) - // this.rd_call_combined(h5_obj, h5_obj_keys, chrom, bin_size, rd_stat) + } this.callers = [] if (wigFeatures_rd_call_combined.length != 0){ @@ -175,7 +180,6 @@ class HDF5Reader { const chrom_dataset = await this.h5_obj.get(mosaic_call_segments) const t0 = Date.now() let chrom_data = await chrom_dataset.value - //console.log(`rd_call_combined ${mosaic_call_segments} ${Date.now() - t0}`) segments = this.decode_segments(chrom_data) } diff --git a/js/cnvpytor/cnvpytorVCF.js b/js/cnvpytor/cnvpytorVCF.js index 2b9d95b10..2ee61e11a 100644 --- a/js/cnvpytor/cnvpytorVCF.js +++ b/js/cnvpytor/cnvpytorVCF.js @@ -7,18 +7,27 @@ function getMean(data) { return (data.reduce(function (a, b) { return a + b; }) / data.length); } + class CNVpytorVCF { + /** + * Create a class instance. + * + * @param {Array} allVariants - An array containing all variants + * @param {number} binSize - The bin size for processing variants. + */ constructor(allVariants, binSize) { this.allVariants = allVariants this.rowBinSize = 10000 this.binSize = binSize - this.binFactor = binSize / this.rowBinSize + this.binFactor = parseInt(binSize / this.rowBinSize) } + /** + * Read rd and BAF information from the vcf file and call accoring to the caller + */ async read_rd_baf(caller='ReadDepth'){ - /* Read rd and BAF information from the vcf file and call accoring to the caller*/ - + // Step1: Parse data from the vcf file; for a fixed rowBinSize var wigFeatures = {} for (let i = this.allVariants.length-1; i >= 0; i--){ @@ -39,7 +48,6 @@ class CNVpytorVCF { dp_count: 0, hets_count:0, hets: [], - //likelihood_score: [], }; } @@ -60,7 +68,6 @@ class CNVpytorVCF { wigFeatures[chr][featureBin].hets.push({ref:ad_a, alt:ad_b}) } - } // get the chromosome names @@ -137,19 +144,13 @@ class CNVpytorVCF { async getAllbins() { const bins = await this.computeDepthFeatures() - - //console.log('getAllbins', bins["value"]) - const fitter = new g_utils.GetFit(bins) - const distParams = fitter.fit_data() - // dconsole.log('rd list', distParams) + return bins } - - formatDataStructure_BAF(wigFeatures, feature_column, scaling_factor=-1) { /* Rescale the BAF level from 0:1 to scaling_factpr:0*/ const baf1 = [] @@ -176,8 +177,14 @@ class CNVpytorVCF { return [baf1, baf2] } + + /** + * Adjust the bin values to actual bin size + * @param {*} wigFeatures - wig features after processing the varaints + * @returns + */ adjust_bin_size(wigFeatures){ - /* adjust the bin values to actual bin size */ + var avgbin = {} for (let chr of this.chromosomes) { if (!avgbin[chr]) { avgbin[chr] = [] } @@ -192,14 +199,15 @@ class CNVpytorVCF { hets_count: 0, binScore: 0, likelihood_score: [], + dp_sum_score: 0 } } - for (var j = k * 10; j < 10 * k + 10; j++) { + for (var j = k * this.binFactor; j < this.binFactor * k + this.binFactor; j++) { if (wigFeatures[chr][j]) { - var tmp_score = parseInt(wigFeatures[chr][j].dp_sum_score / wigFeatures[chr][j].dp_count) * 100; - avgbin[chr][k].binScore += tmp_score; + + avgbin[chr][k].dp_sum_score += wigFeatures[chr][j].dp_sum_score; avgbin[chr][k].dp_count += wigFeatures[chr][j].dp_count avgbin[chr][k].hets_count += wigFeatures[chr][j].hets_count @@ -230,7 +238,7 @@ class CNVpytorVCF { } } } - + avgbin[chr][k].binScore = parseInt(avgbin[chr][k].dp_sum_score / avgbin[chr][k].dp_count) * 100; const updated_bin = this.get_max_min_score(avgbin[chr][k]) avgbin[chr][k].max_likelihood = updated_bin.value } @@ -240,9 +248,8 @@ class CNVpytorVCF { } -function beta(a, b, p, phased = true) { - //p ** a * (1 - p) ** b + p ** b * (1 - p) ** a; - return Math.pow(p, a) * Math.pow(1-p, b) + Math.pow(p, b) * Math.pow(1-p, a) +function beta(a, b, p, phased = true) { + return Math.pow(p, a) * Math.pow(1-p, b) + Math.pow(p, b) * Math.pow(1-p, a) } class MeanShiftCaller{