#+TITLE: Cool Bioinformatics Scripts

* Table of Content :TOC:QUOTE:
#+BEGIN_QUOTE
- [[#qqplot][qqplot]]
- [[#downsampling-bam-files][downsampling bam files]]
- [[#genotype-refinement-using-beagle-3][genotype refinement using beagle 3]]
- [[#calculate-genotype-discordance][calculate genotype discordance]]
- [[#r2-vs-maf-plot][R2 vs MAF plot]]
- [[#phylogenetic-tree][Phylogenetic Tree]]
- [[#fixref][fixref]]
#+END_QUOTE

* [[file:qqplot.py][qqplot]]
You can use make a QQ plot in the following ways.

- one-liner for reading tons of millions of P values from the pipe

#+begin_src shell
# python 
zcat pval.txt.gz | qqplot.py -out test -title "QQ plot on the fly"
# julia (recommand to run it in the REPL)
zcat pval.txt.gz | qqplot.jl --out test --title "QQ plot on the fly"
#+end_src

*warning* : If you have 100 billion P values to process you should definitely use [[qqplot.jl]] instead of [[file:qqplot.py][qqplot.py]]. The hourly processed number of lines of julia version is 5 billion while python is only 700 million on my server.

- running in a julia REPL (recommanded)

#+begin_src julia
include("qqplot.jl")
cmd = pipeline(`zcat pval.gz`, `awk 'NR>1{print $10}'`)
sigp, expp = qqfly("test", cmd=cmd)
#+end_src


- use qqplot.py in your script

#+begin_src python
import numpy as np
from qqplot import qq
p = np.random.random(1000000)
qq(x=p, figname="test.png")
#+end_src

[[file:image/qqplot.png]]

* [[file:downsample.sh][downsampling bam files]]

#+begin_src shell
Usage: downsample.sh [-b <bamlist>] [-d <depth>] [-n <cores>] [-o <outdir>]
#+end_src

* [[file:beagle3-imputation.sh][genotype refinement using beagle 3]]

#+begin_src shell
Usage: beagle3-imputation.sh [options]
Pipeline of genotype refinement for median depth sequencing data using beagle3

-h,          Display help
-i,          Input VCF/BCF file
-o,          Output folder
-f,          MAF filters before imputation
#+end_src

* [[file:calc_imputed_gt_discord.py][calculate genotype discordance]]

When you run imputation analysis with =BEAGLE= (or other imputation tools), you may want to know the distribution of genotype discordance between the original vcf and imputed vcf.

#+begin_src shell
usage: calc_imputed_gt_discord.py [-h] [-chr STRING] VCF1 VCF2 OUT
#+end_src

*warning* : Before running the script, you must be sure the two vcfs have the exact same sites and samples for each chromosome.


[[file:image/calc_imputed_gt_discord.png]]

* [[file:plot-r2-vs-maf.R][R2 vs MAF plot]]
plot INFO/R2 after imputation by =BEAGLE= etc.

[[file:image/r2-vs-maf.png]]

* [[file:njtree.R][Phylogenetic Tree]]

[[file:image/njtree-circular.png]]
* [[file:fixref.py][fixref]]

Before running =bcftools merge=, you maybe need to fix the ref and alt and corresponding genotypes, otherwise =bcftools= will surprise you.

#+begin_src shell
usage: fixref.py [-h] REF_VCF IN_VCF OUT_VCF
#+end_src