joonan-lab / cwas Goto Github PK

zcat jarvis.all_chr.099_.bed.gz | sort -k1,1 -k2,2n -k3,3n | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' > jarvis.all_chr.099_.bed
bedops --partition jarvis.all_chr.099_.bed | bedmap --echo --max --delim '\t' - jarvis.all_chr.099_.bed > nonOverlap.jarvis.all_chr.099_.bed
mv nonOverlap.jarvis.all_chr.099_.bed nonOverlap.jarvis.all_chr.099_.bedGraph
bedGraphToBigWig nonOverlap.jarvis.all_chr.099_.bedGraph hg38.chrom.sizes nonOverlap.jarvis.all_chr.099_.bw

vertebrate.phyloP46way.hg19ToHg38.bw: 2.0
vertebrate.phastCons46way.hg19ToHg38.bw: 0.2
nonOverlap.jarvis.all_chr.09_.bw : 0.9
nonOverlap.jarvis.all_chr.099_.bw : 0.99
nonOverlap.jarvis.all_chr.0999_.bw : 0.999
nonOverlap.Supplementary_Data_2_Constrained.Z4.bw : 4.0

@property
    def categorizer(self) -> Categorizer:
        categorizer = Categorizer(self.category_domain, self.gene_matrix)
        categorizer.phastcons_cutoff = self.annotation_bw_cutoff[
            "phastCons46way"
        ]
        categorizer.phylop_cutoff = self.annotation_bw_cutoff["phyloP46way"]
        categorizer.jarvis09_cutoff = self.annotation_bw_cutoff["JARVIS09"]
        categorizer.jarvis099_cutoff = self.annotation_bw_cutoff["JARVIS099"]
        categorizer.jarvis0999_cutoff = self.annotation_bw_cutoff["JARVIS0999"]
        categorizer.constraintZ4_cutoff = self.annotation_bw_cutoff["ConstraintZ4"]
        return categorizer

class Categorizer:
    def __init__(self, category_domain: dict, gene_matrix: dict) -> None:
        self._category_domain = category_domain
        self._gene_matrix = gene_matrix
        self._default_phylop = -2.0
        self._default_phastcons = 0.0
        self._phylop_cutoff = 2.0
        self._phastcons_cutoff = 0.2
        self._default_jarvis09 = 0.0
        self._default_jarvis099 = 0.0
        self._default_jarvis0999 = 0.0
        self._default_constraintZ4 = 0.0
        self._jarvis09_cutoff = 0.9
        self._jarvis099_cutoff = 0.99
        self._jarvis0999_cutoff = 0.999
        self._constraintZ4_cutoff = 4.0

...

    @property
    def jarvis09_cutoff(self):
        return self._jarvis09_cutoff

    @jarvis09_cutoff.setter
    def jarvis09_cutoff(self, value):
        self._jarvis09_cutoff = value

    @property
    def jarvis099_cutoff(self):
        return self._jarvis099_cutoff

    @jarvis099_cutoff.setter
    def jarvis099_cutoff(self, value):
        self._jarvis099_cutoff = value

    @property
    def jarvis0999_cutoff(self):
        return self._jarvis0999_cutoff

    @jarvis0999_cutoff.setter
    def jarvis0999_cutoff(self, value):
        self._jarvis0999_cutoff = value

    @property
    def constraintZ4_cutoff(self):
        return self._constraintZ4_cutoff

    @constraintZ4_cutoff.setter
    def constraintZ4_cutoff(self, value):
        self._constraintZ4_cutoff = value

...

        jarvis09_conv_func = (
            lambda x: self._default_jarvis09
            if x == ""
            else max(map(float, x.split("&")))
        )
        jarvis099_conv_func = (
            lambda x: self._default_jarvis099
            if x == ""
            else max(map(float, x.split("&")))
        )
        jarvis0999_conv_func = (
            lambda x: self._default_jarvis0999
            if x == ""
            else max(map(float, x.split("&")))
        )
        ConstraintZ4_conv_func = (
            lambda x: self._default_constraintZ4
            if x == ""
            else max(map(float, x.split("&")))
        )

...

        jarvis09_scores = np.vectorize(jarvis09_conv_func)(
            annotated_vcf["JARVIS09"].values
        )
        jarvis099_scores = np.vectorize(jarvis099_conv_func)(
            annotated_vcf["JARVIS099"].values
        )
        jarvis0999_scores = np.vectorize(jarvis0999_conv_func)(
            annotated_vcf["JARVIS0999"].values
        )
        ConstraintZ4_scores = np.vectorize(ConstraintZ4_conv_func)(
            annotated_vcf["ConstraintZ4"].values
        )

...

        is_jarvis09_conservation_arr = (
            jarvis09_scores >= self._jarvis09_cutoff
        ).astype(np.int32)
        is_jarvis099_conservation_arr = (
            jarvis099_scores >= self._jarvis099_cutoff
        ).astype(np.int32)
        is_jarvis0999_conservation_arr = (
            jarvis0999_scores >= self._jarvis0999_cutoff
        ).astype(np.int32)
        is_ConstraintZ4_conservation_arr = (
            ConstraintZ4_scores >= self._constraintZ4_cutoff
        ).astype(np.int32)

...

        annotation_ints += np.vectorize(
            lambda is_jarvis09_conservation: 2
            ** conservation_annotation_idx["JARVIS09"]
            if is_jarvis09_conservation
            else 0
        )(is_jarvis09_conservation_arr)
        annotation_ints += np.vectorize(
            lambda is_jarvis099_conservation: 2
            ** conservation_annotation_idx["JARVIS099"]
            if is_jarvis099_conservation
            else 0
        )(is_jarvis099_conservation_arr)
        annotation_ints += np.vectorize(
            lambda is_jarvis0999_conservation: 2
            ** conservation_annotation_idx["JARVIS0999"]
            if is_jarvis0999_conservation
            else 0
        )(is_jarvis0999_conservation_arr)
        annotation_ints += np.vectorize(
            lambda is_ConstraintZ4_conservation: 2
            ** conservation_annotation_idx["ConstraintZ4"]
            if is_ConstraintZ4_conservation
            else 0
        )(is_ConstraintZ4_conservation_arr)