VEP

cd $HPC_WORK
git clone https://github.com/Ensembl/ensembl-vep.git
cd ensembl-vep
# ---
# release/98
# module load htslib/1.4
# perl INSTALL.pl
# release/104
# long form:
# perl INSTALL.pl --DESTDIR Bio --ASSEMBLY GRCh38 --AUTO acfp --PLUGINS all --SPECIES homo_sapiens,homo_sapiens_merged --NO_TEST --CACHEDIR .vep
# ---
module load htslib-1.9-gcc-5.4.0-p2taavl
# short form:
perl INSTALL.pl -l Bio -y GRCh38 -a acfp -g all -s homo_sapiens,homo_sapiens_merged --NO_TEST -c .vep
ln -sf $HPC_WORK/ensembl-vep/.vep $HOME/.vep
# set up symbolic links to the executables
for f in convert_cache.pl filter_vep haplo variant_recoder vep;
    do ln -sf $HPC_WORK/ensembl-vep/$f $HPC_WORK/bin/$f; done
vep -i examples/homo_sapiens_GRCh37.vcf -o examples/homo_sapiens_GRCh37.txt \
    --force_overwrite --offline --pick --symbol

cd ${HPC_WORK}
git clone https://github.com/Ensembl/VEP_plugins
cd VEP_plugins
# download of database
vep --database --port 3337 --force_overwrite --dir_plugin . --plugin PolyPhen_SIFT,create_db=1
# offline annotation
vep -i homo_sapiens_GRCh37.vcf -o homo_sapiens_GRCh37 --offline --force_overwrite --format vcf --dir_plugin . --plugin PolyPhen_SIFT

export skips=$(grep '##' examples/homo_sapiens_GRCh37.txt | wc -l)
R --no-save -q <<END
  vo <- read.delim("examples/homo_sapiens_GRCh37.txt",skip=as.numeric(Sys.getenv("skips")))
  head(vo)
END

perl -MCPAN -e shell
install Set::IntervalTree

cd examples
bgzip -f homo_sapiens_GRCh37.vcf
tabix -Cf homo_sapiens_GRCh37.vcf.gz
tabix -h  homo_sapiens_GRCh37.vcf.gz 22:50616005-50616006 | \
vep --cache --fork 4 --port 3337 --format vcf -o - --tab --no_stats | \
grep -v '##'
cd -

export ref=~/rds/post_qc_data/interval/reference_files/genetic/reference_files_genotyped_imputed/
export chunk_size=10000
seq 22 | \
parallel -j1 --env ref -C' ' '
  export n=$(awk "END{print NR-1}" ${ref}/impute_{}_interval.snpstats)
  export g=$(expr ${n} / ${chunk_size})
  export s=$(expr $n - $(( $g * $chunk_size)))
  (
    for i in $(seq ${g}); do
      (
        awk "BEGIN{print \"##fileformat=VCFv4.0\"}"
        awk -vOFS="\t" "BEGIN{print \"#CHROM\",\"POS\",\"ID\",\"REF\",\"ALT\",\"QUAL\",\"FILTER\",\"INFO\"}"
        (
          sed "1d" ${ref}/impute_{}_interval.snpstats | \
          awk -v i=${i} -v chunk_size=${chunk_size} -v OFS="\t" "NR==(i-1)*chunk_size+1,NR==i*chunk_size {
              if(\$1==\".\") \$1=\$3+0 \":\" \$4 \"_\" \$5 \"/\" \$6; print \$3+0,\$4,\$1,\$5,\$6,\".\",\".\",\$19}"
          if [ ${s} -gt 0 ] &&  [ ${i} -eq ${g} ]; then
             sed "1d" ${ref}/impute_{}_interval.snpstats | \
             awk -v i=${i} -v chunk_size=${chunk_size} -v OFS="\t" -v n=${n} "NR==i*chunk_size+1,NR==n {
                 if(\$1==\".\") \$1=\$3+0 \":\" \$4 \"_\" \$5 \"/\" \$6; print \$3+0,\$4,\$1,\$5,\$6,\".\",\".\",\$19}"
          fi
        )
      ) | \
      vep  --cache --offline --format vcf -o - --tab --pick --no_stats  \
           --species homo_sapiens --assembly GRCh37 --port 3337 | \
      (
        if [ ${i} -eq 1 ]; then cat; else grep -v "#"; fi
      )
    done
  ) | \
  gzip -f > work/INTERVAL-{}.vep.gz
'

wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/ensemblToGeneName.txt.gz
zgrep -e ENST00000368485 -e ENST00000515190 ensemblToGeneName.txt.gz

R -q -e "grex::grex(\"ENSG00000160712\")"

mysql  --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 <<END
select distinct G.gene,N.value from ensGtp as G, ensemblToGeneName as N where
   G.transcript=N.name and
   G.gene in ("ENSG00000160712");
END

mysql -h ensembldb.ensembl.org --port 5306  -u anonymous -D homo_sapiens_core_64_37 -A <<END
select distinct
   G.stable_id,
   S.synonym
from
  gene_stable_id as G,
  object_xref as OX,
  external_synonym as S,
  xref as X ,
  external_db as D
where
  D.external_db_id=X.external_db_id and
  X.xref_id=S.xref_id and
  OX.xref_id=X.xref_id and
  OX.ensembl_object_type="Gene" and
  G.gene_id=OX.ensembl_id and
  G.stable_id in ("ENSG00000160712");
END

# Compressed VCF file/Index file
# GCRh37
curl ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar.vcf.gz -o clinvar_GRCh37.vcf.gz
curl ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar.vcf.gz.tbi -o clinvar_GRCh37.vcf.gz.tbi
# GCRh38
curl ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz -o clinvar_GRCh38.vcf.gz
curl ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz.tbi -o clinvar_GRCh38.vcf.gz.tbi

vep --id "1 154426970 154426970 A/C 1" --species homo_sapiens -o rs2228145 --cache --offline --force_overwrite \
    --assembly GRCh37 --custom clinvar_GRCh37.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN

vep --id "1 154426970 154426970 A/C 1" --species homo_sapiens -o rs2228145 --cache --offline --force_overwrite \
    --custom clinvar_GRCh37.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN --tab \
    --fields Uploaded_variation,Gene,Consequence,ClinVar_CLNSIG,ClinVar_CLNREVSTAT,ClinVar_CLNDN

wget ftp://dbnsfp:dbnsfp@dbnsfp.softgenetics.com/dbNSFP4.1a.zip
unzip dbNSFP4.1a.zip -d dbNSFP4.1a
cd dbNSFP4.1a
# code for 4.0a
# zcat dbNSFP4.0a_variant.chr1.gz | head -n1 > h
# zgrep -h -v ^#chr dbNSFP4.0a_variant.chr* | sort -k1,1 -k2,2n - | cat h - | bgzip -c > dbNSFP4.0a.gz
# efficient version
(
  zcat dbNSFP4.1a_variant.chr1.gz | head -n1
  export prefix=dbNSFP4.1a_variant.chr
  echo ${prefix}{1..22}.gz ${prefix}M.gz ${prefix}X.gz ${prefix}Y.gz | \
  xargs -I {} bash -c "zcat {} | sed '1d'"
) | bgzip -c > dbNSFP4.1a.gz
tabix -s 1 -b 2 -e 2 dbNSFP4.1a.gz
cd -
vep -i examples/homo_sapiens_GRCh37.vcf -o test --cache --force_overwrite --offline \
    --plugin dbNSFP,dbNSFP4.1a/dbNSFP4.1a.gz,LRT_score,FATHMM_score,MutationTaster_score

java -jar search_dbNSFP41a.jar -i tryhg19.in -o tryhg19.out -v hg19
java -jar search_dbNSFP41a.jar -i tryhg38.in -o tryhg38.out

wget -qO- ftp://ftp.ccb.jhu.edu/pub/software/genesplicer/GeneSplicer.tar.gz | \
tar xvfz -
mv GeneSplicer.pm ~/.vep/Plugins
./vep -i variants.vcf --plugin GeneSplicer,[path_to_genesplicer_bin],[path_to_training_dir],[option1=value],[option2=value]

cd loftee
# human_ancestor_fa
## samtools --version gives 0.1.19
wget http://www.broadinstitute.org/~konradk/loftee/human_ancestor.fa.rz
wget http://www.broadinstitute.org/~konradk/loftee/human_ancestor.fa.rz.fai
## samtools --version gives 1.x
wget https://s3.amazonaws.com/bcbio_nextgen/human_ancestor.fa.gz
wget https://s3.amazonaws.com/bcbio_nextgen/human_ancestor.fa.gz.fai
wget https://s3.amazonaws.com/bcbio_nextgen/human_ancestor.fa.gz.gzi
# conservation_file -- note the data, schema and GERP files are required only to build the sql file
wget -qO- https://personal.broadinstitute.org/konradk/loftee_data/GRCh37/phylocsf_gerp.sql.gz | \
gunzip -c > phylocsf_gerp.sql
# wget https://www.broadinstitute.org/~konradk/loftee/phylocsf_data.tsv.gz
# wget https://www.broadinstitute.org/~konradk/loftee/phylocsf_data_schema.sql
# wget https://personal.broadinstitute.org/konradk/loftee_data/GRCh37/GERP_scores.final.sorted.txt.gz
# wget https://personal.broadinstitute.org/konradk/loftee_data/GRCh37/GERP_scores.exons.txt.gz
# annotation
vep --id "1 154426970 154426970 A/C 1" --species homo_sapiens -o rs2228145 --cache --offline --force_overwrite \
    --assembly GRCh37 --plugin LoF,loftee_path:.,human_ancestor_fa:human_ancestor.fa.gz,conservation_file:phylocsf_gerp.sql.gz
# VEP documentation example
vep --input_file homo_sapiens_GRCh37.vcf --output_file test --cache --dir_cache ${HPC_WORK}/ensembl-vep/.vep --dir_plugins ${HPC_WORK}/loftee --offline \
    --pick --force_overwrite --species homo_sapiens --assembly GRCh37 \
    --plugin LoF,loftee_path:.,human_ancestor_fa:human_ancestor.fa.gz,conservation_file:phylocsf_gerp.sql

# error
samtools faidx human_ancestor.fa.rz
# rework
gunzip -c human_ancestor.fa.rz | bgzip -c > output.fa.gz
# check and put it back
ll human_ancestor.fa.rz output.fa.gz
mv output.fa.gz human_ancestor.fa.rz
# now possible and check
samtools faidx human_ancestor.fa.rz
ll human_ancestor.fa.*

# 1. Download and extract source code
cd $HOME
wget -qO- https://github.com/ucscGenomeBrowser/kent/archive/v335_base.tar.gz | \
tar xzf -
# 2. Set up compiling flags
export KENT_SRC=$HOME/kent-335_base/src
export MACHTYPE=$(uname -m)
export CFLAGS="-fPIC"
export MYSQLINC=`mysql_config --include | sed -e 's/^-I//g'`
export MYSQLLIBS=`mysql_config --libs`
# 3. Amend parameters
cd $KENT_SRC/lib
echo 'CFLAGS="-fPIC"' > ../inc/localEnvironment.mk
# 4. Build library
make clean && make
cd ../jkOwnLib
make clean && make
# 5. On Mac OSX
ln -s $KENT_SRC/lib/x86_64/* $KENT_SRC/lib/
# 6. Install Perl modules
cd ${HPC_WORK}/ensembl-vep
cpan JSON
cpan Bio::DB::BigFile
# 7. Test
perl -Imodules t/AnnotationSource_File_BigWig.t

# GRCh38 files
# wget https://personal.broadinstitute.org/konradk/loftee_data/GRCh38/gerp_conservation_scores.homo_sapiens.GRCh38.bw
# wget https://personal.broadinstitute.org/konradk/loftee_data/GRCh38/human_ancestor.fa.gz
# wget https://personal.broadinstitute.org/konradk/loftee_data/GRCh38/human_ancestor.fa.gz.fai
# wget https://personal.broadinstitute.org/konradk/loftee_data/GRCh38/human_ancestor.fa.gz.gzi
# wget https://personal.broadinstitute.org/konradk/loftee_data/GRCh38/loftee.sql.gz

# Annotation
export ENSEMBL=~/rds/rds-jmmh2-public_databases/ensembl-vep
export LOFTEE38=${ENSEMBL}/loftee/loftee_data/GRCh38
export LOFTEE38GERP=${LOFTEE38}/gerp_conservation_scores.homo_sapiens.GRCh38.bw
export LOFTEE38HA=${LOFTEE38}/human_ancestor.fa.gz
export LOFTEE38SQL=${LOFTEE38}/loftee.sql
export rds=..  # ~/rds/rds-jmmh2-public_databases/ensembl-vep will be user-specific
vep --input_file for_VEP.txt --format ensembl --output_file ${rds}/for_VEP_output2.txt --force_overwrite --offline --symbol --merged \
    --fasta Homo_sapiens.GRCh38.dna.toplevel.fa --dir_cache ${rds}/.vep --dir_plugins . \
    --protein --symbol --tsl --canonical --mane_select --biotype --check_existing --sift b --polyphen b \
    --plugin LoF,loftee_path:.,gerp_bigwig:${LOFTEE38GERP},human_ancestor_fa:${LOFTEE38HA},conservation_file:${LOFTEE38SQL}

mkdir REVEL
cd REVEL
wget https://rothsj06.u.hpc.mssm.edu/revel-v1.3_all_chromosomes.zip
unzip revel-v1.3_all_chromosomes.zip
cat revel_with_transcript_ids | tr "," "\t" > tabbed_revel.tsv
sed '1s/.*/#&/' tabbed_revel.tsv > new_tabbed_revel.tsv
bgzip new_tabbed_revel.tsv
# GRCh37:
tabix -f -s 1 -b 2 -e 2 new_tabbed_revel.tsv.gz
# GRCh38:
zcat new_tabbed_revel.tsv.gz | head -n1 > h
zgrep -h -v ^#chr new_tabbed_revel.tsv.gz | awk '$3 != "." ' | sort -k1,1 -k3,3n - | cat h - | bgzip -c > new_tabbed_revel_grch38.tsv.gz
tabix -f -s 1 -b 3 -e 3 new_tabbed_revel_grch38.tsv.gz

vep --input_file <input> --plugin REVEL,${ENSEMBL}/.vep/Plugins/new_tabbed_revel.tsv.gz