search

statgen / swiss

Software to help identify overlap between association scan results and GWAS hit catalogs.
Other
14 stars 5 forks source link

Hg38 swiss-create-data issue and raising the p-value #18

Open smcamcamca opened 2 years ago

smcamcamca commented 2 years ago

Command-line used to cause the issue

swiss-create-data --genome-build GRCh38p12 --dbsnp-build b151

Problem description

Hello,

I am trying to loosen the p-value threshold for GWAS catalog variants to 5e-05 using the --gwas-cat-p option. It runs fine, but I noticed the output is identical to using the default of 5e-08 (the lowest/least significant GWAS_LOG_PVAL was 7.3 for both runs). This was the command I used: 

swiss --assoc input_file.tsv --delim tab --build hg38  --gwas-cat-p 5e-05  --out output_file

I assumed this might be related to filtering done on the data provided by swiss ---download-data, so I tried to download my own hg38 genome build using swiss-create-data but I get the error below. I have also tried GRCh38p13 and GRCh38 with the same issue.

Thank you!

Program output

Genome build: GRCh38p12 dbSNP build: b151 Downloading rsID merge history... RsMergeArch.bcp.gz: 153MB [00:32, 4.78MB/s]
Downloading rsID deletion history... SNPHistory.bcp.gz: 106MB [00:32, 3.25MB/s]
Downloading NCBI dbSNP VCF for b151 / GRCh38p12 @ ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606_b151_GRCh38p12/VCF/00-All.vcf.gz 00-All.vcf.gz: 0.00B [00:00, ?B/s] Traceback (most recent call last): File "/home/user/.conda/envs/gwas_swiss/bin/swiss-create-data", line 22, in main() File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/site-packages/swiss/create_data.py", line 671, in main vcf_path = download_dbsnp_vcf(opts.dbsnp_build,opts.genome_build) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/site-packages/swiss/create_data.py", line 129, in download_dbsnp_vcf urlretrieve(url,filename=outpath,reporthook=tqdm_hook(t),data=None) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/urllib.py", line 98, in urlretrieve return opener.retrieve(url, filename, reporthook, data) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/urllib.py", line 247, in retrieve fp = self.open(url, data) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/urllib.py", line 215, in open return getattr(self, name)(url) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/urllib.py", line 552, in open_ftp ftpwrapper(user, passwd, host, port, dirs) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/urllib.py", line 879, in init self.init() File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/urllib.py", line 891, in init self.ftp.cwd(_target) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/ftplib.py", line 574, in cwd return self.voidcmd(cmd) File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/ftplib.py", line 256, in voidcmd return self.voidresp() File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/ftplib.py", line 231, in voidresp resp = self.getresp() File "/home/user/.conda/envs/gwas_swiss/lib/python2.7/ftplib.py", line 226, in getresp raise error_perm, resp IOError: [Errno ftp error] 550 snp/organisms/human_9606_b151_GRCh38p12/VCF: No such file or directory

Swiss version

1.1.1

welchr commented 2 years ago

Looks like there's only a GRCh38p7 on their FTP. The VCF appears to be rather out of date at this point (2018).

They seem to have changed how they do dbSNP VCF releases, with the latest VCF now at:

https://ftp.ncbi.nih.gov/snp/redesign/latest_release/VCF/GCF_000001405.39.gz (GRCh38p13, dbSNP 155)

Unfortunately the new VCF is a bit different and wouldn't be parsed correctly by the program. In particular they've switched to using RefSeq sequence IDs for chromosomes.