Skip to content

Commit

Permalink
cnvpytor track: bug fixes and use baf_i1 signal
Browse files Browse the repository at this point in the history
  • Loading branch information
@arpanda
arpanda committed Nov 12, 2023
1 parent ad7c742 commit 6ac2ae0
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 29 deletions.
22 changes: 13 additions & 9 deletions js/cnvpytor/HDF5IndexedReader.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`
}
Expand Down Expand Up @@ -49,19 +50,21 @@ 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
}

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);
Expand Down Expand Up @@ -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){
Expand Down Expand Up @@ -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)

}
Expand Down
47 changes: 27 additions & 20 deletions js/cnvpytor/cnvpytorVCF.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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--){
Expand All @@ -39,7 +48,6 @@ class CNVpytorVCF {
dp_count: 0,
hets_count:0,
hets: [],
//likelihood_score: [],
};
}

Expand All @@ -60,7 +68,6 @@ class CNVpytorVCF {
wigFeatures[chr][featureBin].hets.push({ref:ad_a, alt:ad_b})
}


}

// get the chromosome names
Expand Down Expand Up @@ -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 = []
Expand All @@ -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] = [] }
Expand All @@ -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

Expand Down Expand Up @@ -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
}
Expand All @@ -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{
Expand Down

0 comments on commit 6ac2ae0

Please sign in to comment.