Tutorial: How to make an arcplot

Are you dazzled by the beautiful colors and arcs of the arcplots in our preprint? Example commands are shown here to learn how to make your own for the gene HSPB1 with a 5kb window on either side with data from the 1000 Genomes Project.

1.) Reformat variants in your region of interest.

python ~/AlleleAnalyzer/preprocessing/generate_gens_dfs/get_gens_df.py ~/ALL.chr7_GRCh38.genotypes.20170504.bcf 7:76297558-76309301 hspb1_gens -v

2.) Annotate variants based on which type(s) of Cas you want to use.

python ~/AlleleAnalyzer/preprocessing/annotate_variants/annot_variants.py hspb1_gens.h5 SpCas9 ~/excisionFinderData_public/hg38_pams/ ~/hg38.fa hspb1_annots_spcas9

3.) Run ExcisionFinder in exhaustive mode to get all possible targetable pairs.

python ~/AlleleAnalyzer/scripts/ExcisionFinder.py -v ~/gene_list_hg38.tsv HSPB1 hspb1_annots_spcas9.h5 10000 SpCas9 ~/ALL.chr7_GRCh38.genotypes.20170504.bcf hspb1_combos --exhaustive --window=5000

4.) Generate arcplot input file.

Note: you can download the 1kg_allsamples.tsv here. It's the samples file from the 1000 Genomes Project.

python ~/AlleleAnalyzer/plotting_scripts/gen_arcplot_input.py hspb1_combos_exh.tsv hspb1_arcplot_in --sample_legend=~/1kg_allsamples.tsv

.) Plot arcplot.

Rscript ~/AlleleAnalyzer/plotting_scripts/arcplot_black.R hspb1_arcplot_in.tsv hspb1_arcplot_10 10 76302558 76304301 5000 hspb1_filt_10