CookHLA

SNP2HLA

Web: SNP2HLA v1.0.3 (BEAGLE utitlities)

Installation

It turns out to be difficult to download beagle 3.0.4 as indicated so it is included in the files/ directory (beagle 3.0.4, documentation, example).

wget -qO- https://software.broadinstitute.org/mpg/snp2hla/data/SNP2HLA_package_v1.0.3.tar.gz | \
tar xvfz -
cd SNP2HLA_package_v1.0.3/SNP2HLA
# add beagle2linkage.jar as above
# test.sh is adapted from SNP2HLA.csh by removing argument checking and as an executable.
# The .dos file described in README.txt is actually the .dosage file generated internally from test.sh
# The association statistics will be in .assoc.assoc.dosage; the double .assoc guarantees an .assoc.log file.
csh SNP2HLA.csh 1958BC HM_CEU_REF 1958BC_IMPUTED plink
csh ParseDosage.csh 1958BC_IMPUTED.bgl.gprobs > 1958BC_IMPUTED.bgl.dos
plink --noweb --dosage 1958BC_IMPUTED.dosage noheader format=1 --fam 1958BC_IMPUTED.fam --logistic --out 1958BC_IMPUTED.assoc

As indicated, we call the .csh scripts with csh which is actually tsch on csd3.

Example

The screen output is as follows,

  SNP2HLA: Performing HLA imputation for dataset 1958BC
- Java memory = 2000Mb
- Beagle window size = 1000 markers
[1] Extracting SNPs from the MHC.
[2] Performing SNP quality control.
[3] Convering data to beagle format.
[4] Performing HLA imputation (see 1958BC_IMPUTED.bgl.log for progress).
[5] Converting posterior probabilities to PLINK dosage format.
[6] Converting imputation genotypes to PLINK .ped format.
DONE!

The output is 1958BC_IMPUTED in PLINK binary format.

MakeReference

We first create MakeReference.tcsh such that calls to .pl scripts are prefixed with perl, so that

diff MakeReference.csh MakeReference.tcsh

77c77
<     ./HLAtoSequences.pl $HLA_DATA HLA_DICTIONARY_AA.txt AA > $OUTPUT.AA.ped
---
>     perl HLAtoSequences.pl $HLA_DATA HLA_DICTIONARY_AA.txt AA > $OUTPUT.AA.ped
81c81
<     ./encodeVariants.pl $OUTPUT.AA.ped $OUTPUT.AA.map $OUTPUT.AA.CODED
---
>     perl encodeVariants.pl $OUTPUT.AA.ped $OUTPUT.AA.map $OUTPUT.AA.CODED
93c93
<     ./encodeHLA.pl $HLA_DATA $OUTPUT.HLA.map > $OUTPUT.HLA.ped
---
>     perl encodeHLA.pl $HLA_DATA $OUTPUT.HLA.map > $OUTPUT.HLA.ped
100c100
<     ./HLAtoSequences.pl $HLA_DATA HLA_DICTIONARY_SNPS.txt SNPS > $OUTPUT.SNPS.ped
---
>     perl HLAtoSequences.pl $HLA_DATA HLA_DICTIONARY_SNPS.txt SNPS > $OUTPUT.SNPS.ped
104c104
<     ./encodeVariants.pl $OUTPUT.SNPS.ped $OUTPUT.SNPS.map $OUTPUT.SNPS.CODED
---

tcsh MakeReference.tcsh HAPMAP_CEU HAPMAP_CEU_HLA.ped HM_CEU_REF plink

tcsh MakeReference.tcsh HAPMAP_CEU HAPMAP_CEU_HLA.ped HM_CEU_REF plink
Creating reference panel: HM_CEU_REF
[1] Generating amino acid sequences from HLA types.
[2] Encoding amino acids positions.
[3] Encoding HLA alleles.
[4] Generating DNA sequences from HLA types.
[5] Encoding SNP positions.
[6] Extracting founders.
[7] Merging SNP, HLA, and amino acid datasets.
Warning: Variants 'SNP_A_30018350' and 'AA_A_-11_30018350' have the same
position.
Warning: Variants 'SNP_A_30018461_G' and 'SNP_A_30018461_A' have the same
position.
Warning: Variants 'SNP_A_30018461_T' and 'SNP_A_30018461_G' have the same
position.
987 more same-position warnings: see log file.
[8] Performing quality control.
[9] Preparing files for Beagle.
[10] Converting to beagle format.
[11] Phasing reference using Beagle (see progress in HM_CEU_REF.bgl.log).
[12] Done.

↩

HATK

Web: https://github.com/WansonChoi/HATK

Installation

git clone https://github.com/WansonChoi/HATK

Example

source ~/COVID-19/py37/bin/activate
python3 HATK.py \
        --variants example/wtccc_filtered_58C_RA.hatk.300+300.chr6.hg18 \
        --hped example/wtccc_filtered_58C_RA.hatk.300+300.hped \
        --2field \
        --pheno example/wtccc_filtered_58C_RA.hatk.300+300.phe \
        --pheno-name RA \
        --out work/RESULT_EXAMPLE_wtccc_filtered_58C_RA.hatk.300+300.chr6.hg18 \
        --imgt 3320 \
        --hg 18 \
        --imgt-dir example/IMGTHLA3320 \
        --multiprocess 2

This is from the documentation, where --variants reads in the genotype files and --hped the .hped file to be followed by specification of the RA phenotype in a logistic regression. Note that the example is more desirable compared to the toy data in SNP2HLA given its 600 samples and 29,373 variants; we proceed with the imputation with the 1000Genomes panel provided with CookHLA.

export hatk=${HPC_WORK}/HATK
export cookhla=${HPC_WORK}/CookHLA

cd ${cookhla}
python CookHLA.py \
       -i ${hatk}/example/wtccc_filtered_58C_RA.hatk.300+300.chr6.hg18 \
       -hg 18 \
       -o ${hatk}/work/hla_CookHLA \
       -ref ${cookhla}/1000G_REF/1000G_REF.EUR.chr6.hg18.29mb-34mb.inT1DGC \
       -gm ${caprion}/hla_IMPUTED.mach_step.avg.clpsB \
       -ae ${caprion}/hla_IMPUTED.aver.erate \
       -mem 20g \
       -mp 8
cd -

↩

HLA-TAPAS

Web: https://github.com/immunogenomics/HLA-TAPAS

Installation

git clone https://github.com/immunogenomics/HLA-TAPAS

Example

python HLA-TAPAS.py \
       --target example/Case+Control.300+300.chr6.hg18 \
       --reference example/1000G.EUR.chr6.hg18.28mb-35mb \
       --hped-Ggroup example/1000G.EUR.Ggroup.hped \
       --pheno example/Case+Control.300+300.phe \
       --hg 18 \
       --out MyHLA-TAPAS/Case+Control+1000G_EUR_REF \
       --mem 4g \
       --nthreads 4

↩

HIBAG

Web: https://hibag.s3.amazonaws.com/index.html (Bioconductor)

Installation

Rscript -e 'BiocManager:install("HIBAG")'

Examples

See Amazon but most are available from GitHub.

↩