Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

v3 subset script #566

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

v3 subset script #566

wants to merge 3 commits into from

Conversation

mike-w-wilson
Copy link
Contributor

This is the most updated subset script. I went ahead and removed subset specific callstats. We've deprioritized the old subset script for years and eventually moved it into a draft after it became outdated but we should get this in. This is the script, with the exception of the subset callstats removal, that was last run on v3 subset requests.

Copy link
Contributor

@KoalaQin KoalaQin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the subsetting only working on the subsets already exist in freq_meta or any list of samples? It seems missing to me if you are extracting call_stats for a subset in the code.

vd.row_key, variant_qc_annotations
)
)
logger.info("Dropping VRS as it is VCF incompatible as structured.")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
logger.info("Dropping VRS as it is VCF incompatible as structured.")
logger.info("Dropping VRS because its current structure is incompatible with VCF.")

Comment on lines +127 to +132
vd = vd.annotate_rows(
n_unsplit_alleles=hl.len(vd.alleles),
mixed_site=(hl.len(vd.alleles) > 2)
& hl.any(lambda a: hl.is_indel(vd.alleles[0], a), vd.alleles[1:])
& hl.any(lambda a: hl.is_snp(vd.alleles[0], a), vd.alleles[1:]),
)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why this extra step? Could you add a comment?

if vds_output:
vds.write(f"{output_path}/subset.vds", overwrite=args.overwrite)

if vcf: # TODO: Do we want to add in sharded vcf delivery using n_partitions?
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add an argument as an option?
But aren't we getting a VCF per Chr when subsetting?

"--variant-qc-annotations",
nargs="+",
type=str,
choices=VARIANT_QC_ANNOTATIONS,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is choices necessary here?


def compute_partitions(
mt, entry_size=3.5, partition_size=128000000, min_partitions=20
) -> int: # TODO: move to gnomad_methods?
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think this function should be in gnomad_methods/utils.py. It's very useful.

vds_output = args.vds
dense_mt = args.dense_mt
split_multi = args.split_multi
subset_call_stats = args.subset_call_stats
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How is this used? It seems that we are still copying all info fields?

if add_variant_qc:
logger.info("Adding variant QC annotations.")
vd = vd.annotate_rows(
**make_variant_qc_annotations_dict(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function from v4 is using v4 release_sites HT, but raw_qual_hists is now under histograms, maybe copy that function here, and use v3 release_sites HT?
It will run into error otherwise.

Copy link
Contributor

@KoalaQin KoalaQin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm checking a subsetting code for Cal and Bob, I found that if they are subsetting some non-releasable samples, I don't think we can still use release HT to get variant QC annotations, because some variants won't be in release HT>

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants