Get information about direction of agreement for any two lists of SNPs, based on LDlink and plink
Compares every SNP in one list to each SNP in another:
Pairs are the filtered by distance on the genome and linkage disequilibrium. The filtered list includes phasing information, allowing instant lookup of the allele in one SNP given an allele in the other. Lookup for a single pair is also possible.
The software works best used as a python library, but two scripts are provided also:
ldlists
: allows comparison of two SNP lists as described aboveldpair
: allows comparison of just one pair with a more detailed report than given byldlists
To run this code, you will need a copy of a dbSNP sqlite database and a copy of all the 1000genomes VCF files.
Note: Currently under development and unstable/unusable.
The code can be used as
To run ldlists
you will need a funtional version of PLINK 1.9: https://www.cog-genomics.org/plink, to run ldpair
you will need a functional version of tabix htslib
I recommend installing both of these tools.
All python requirements will be automatically downloaded when installing using the setup.py file.
Essential requirements:
Required for table handling in ldlists:
Recommended for progress bars:
Recommended if you have access to a torque or slurm cluster:
All of the above except fyrd will be installed automatically, if you want to parallelize on a cluster, install and configure fyrd on your system.
To run this code you will need the 1000genomes data in plink format and a copy of dbSNP. It could hypothetically be generalized to run on any plink file and without dbSNP.
To get the 1000genomes data, download all the ALL.chr*vcf.gz files on ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502.
To convert to plink, you will need to run the following commands on every file:
i=ALL.chr6.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz;
j=$(echo $i | sed 's/.vcf.gz//');
vftools --gzvcf $i --plink-tped --out $j &&
plink --tped $j.tped --tfam $j.tfam --out $j
To run these on every file:
for i in ALL.chr*.vcf.gz; do
j=$(echo $i | sed 's/.vcf.gz//');
vftools --gzvcf $i --plink-tped --out $j &&
plink --tped $j.tped --tfam $j.tfam --out $j
done
We use a dbSNP database interface that I wrote for position lookup: https://github.com/MikeDacre/dbSNP You can download any version of dbSNP you want and build a database using that code. However, this is a very slow process. You can find pre-built databases at sftp [email protected]/dbSNP, password is 'givemethedb'. Note that we support sftp only, not regular FTP, and this is just my personal machine at Stanford, if there are too many downloads I will have to stop providing data.