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

could you give us an example of cactus? #4

Closed
xuxingyubio opened this issue May 22, 2023 · 7 comments
Closed

could you give us an example of cactus? #4

xuxingyubio opened this issue May 22, 2023 · 7 comments

Comments

@xuxingyubio
Copy link

No description provided.

@xuxingyubio xuxingyubio changed the title Sorry to bother you, could you give me an example of cactus? I downloaded the pan genome constructed using cactus from HPRC, and it seems that there are some differences between it and pggb when using this software. As I have just come into contact with this field, I don't know how to solve it. could you give us an example of cactus? May 22, 2023
@xuxingyubio
Copy link
Author

Sorry to bother you, could you give me an example of cactus? I downloaded the pan genome constructed using cactus from HPRC, and it seems that there are some differences between it and pggb when using this software. As I have just come into contact with this field, I don't know how to solve it.

@danydoerr
Copy link
Member

Sure! I guess it breaks when grouping paths together, because cactus uses both "P" and "W" lines. I'll soon make some convenience options available in panacus that will allow grouping by sample/haplotype without providing an explicit group list. But until then, here's a recipe that works with the HPRC MC graph and other graphs that use W-lines:

  1. Download and unpack the graph:
https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.0-mc-grch38.gfa.gz
gunzip hprc-v1.0-mc-grch38.gfa.gz
  1. Prepare file to group paths by sample: <-- this is the crucial step!
grep '^P' hprc-v1.0-mc-grch38.gfa | cut -f2 > hprc-v1.0-mc-grch38.paths.txt
grep -e '^W' hprc-v1.0-mc-grch38.gfa | cut -f2-6 | awk '{ print $1 "#" $2 "#" $3 ":" $4 "-" $5 }' >> hprc-v1.0-mc-grch38.paths.txt
cut -f1 -d\# hprc-v1.0-mc-grch38.paths.txt > hprc-v1.0-mc-grch38.groupnames.txt
paste hprc-v1.0-mc-grch38.paths.txt hprc-v1.0-mc-grch38.groupnames.txt > hprc-v1.0-mc-grch38.groups.txt
  1. Prepare file to select subset of paths corresponding to haplotypes: <-- in the MC graph, reference paths are upper-cased
grep -ive 'grch38\|chm13' hprc-v1.0-mc-grch38.paths.txt > hprc-v1.0-mc-grch38.paths.haplotypes.txt
  1. Run panacus histgrowth to calculate pangenome growth for nodes (default) with quorum tresholds 0, 1, 0.5, and 0.1 using up to 16 threads:
RUST_LOG=info panacus histgrowth -t16 -q 0,1,0.5,0.1 -g hprc-v1.0-mc-grch38.groups.txt -s hprc-v1.0-mc-grch38.paths.haplotypes.txt hprc-v1.0-mc-grch38.gfa > hprc-v1.0-mc-grch38.histgrowth.node.txt

hprc-v1.0-mc-grch38.histgrowth.node.pdf
hprc-v1.0-mc-grch38.histgrowth.node.txt

@danydoerr danydoerr pinned this issue May 22, 2023
@xuxingyubio
Copy link
Author

Thank you for your answer.
How to quantify the amount of non-reference(GRCh38) sequence has always troubled me.
Does the above step indicate this information, or does it require some processing of the constructed pangenome?

@danydoerr
Copy link
Member

danydoerr commented May 22, 2023

No processing required, panacus can do that for you if you hand it an exclusion file. The part of the graph that is traversed by any of the specified paths in the exclusion file will be omitted from the computation. So, staying with this example, to quantify non-reference (GRCh38) sequence (bp), do the following:

  1. Produce exclusion file
grep -ie 'grch38' hprc-v1.0-mc-grch38.paths.txt > hprc-v1.0-mc-grch38.paths.grch38.txt
  1. Run panacus histrowth:
RUST_LOG=info panacus histgrowth -t16 -c bp -q 0,1,0.5,0.1 -e hprc-v1.0-mc-grch38.paths.grch38.txt -g hprc-v1.0-mc-grch38.groups.txt -s hprc-v1.0-mc-grch38.paths.haplotypes.txt hprc-v1.0-mc-grch38.gfa > hprc-v1.0-mc-grch38.histgrowth.nogrch38.bp.txt

@xuxingyubio
Copy link
Author

In ‘A draft human pangenome reference’,Fig.3g shows that when the first sample is added, the total length is about 30Mb. When I used the above method, its length seemed to be a haplotype length.How has it been processed?

@danydoerr
Copy link
Member

Fig. 3g is an ordered growth histogram (which you can also produce with this tool). But you can also do this by looking at all possible permutations, which is what panacus histgrowth does. If you follow the steps from above, you should be able to get the attached data / plot. As expected, the average 1st genome adds about 30Mb of non-reference sequence to the pangenome. The data is consistent with Fig. 3g from the paper.

hprc-v1.0-mc-grch38.histgrowth.nogrch38.bp.pdf
hprc-v1.0-mc-grch38.histgrowth.nogrch38.bp.txt

@danydoerr
Copy link
Member

A new release is out that fixes several bugs in the software. If your issue still persists, please open a new issue--I'm closing this one.

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

No branches or pull requests

2 participants