最近在做遺傳統(tǒng)計學的作業(yè)，卻慘遭ensembl背刺
具體表現(xiàn)為線上工具連接1kGenome數(shù)據(jù)庫失效
Error with message:
cannot open https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel No such file or directory at /net/isilonP/public/ro/ensweb-software/sharedsw/e108/1000G-tools/vcftoped/vcftoped.pl line 268.
不過在工具的index頁面還看到了一個pl版本,并且提供一些可以設置的參數(shù)和一個例子
想想要不就把vcf和panel文件下載到本地，用.pl手動轉(zhuǎn)換吧
此處以RET(10:43077069-43130351)為例
數(shù)據(jù)
vcf(chr10)下載地址: https://ftp.ensembl.org/pub/data_files/homo_sapiens/GRCh38/variation_genotype/ALL.chr10_GRCh38.genotypes.20170504.vcf.gz
tbi(chr10)下載地址: https://ftp.ensembl.org/pub/data_files/homo_sapiens/GRCh38/variation_genotype/ALL.chr10_GRCh38.genotypes.20170504.vcf.gz.tbi
panel下載地址：https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel
tbi是索引文件，如果工作路徑?jīng)]有與.vcf.gz相應的.tbi，這個perl腳本會要求用-tabix提供索引文件路徑，從而調(diào)用tabix，進一步調(diào)用bcftools來創(chuàng)建.vcf.gz相應的.tbi。費時費力 而且中間裝軟件設置環(huán)境又是100年(別問我怎么知道的)
運行
注意：腳本要求perl環(huán)境5以上 我使用的是conda-forge里能安裝的最新的5.34
perl vcf_to_ped_convert.pl \
-vcf ALL.chr10_GRCh38.genotypes.20170504.vcf.gz \
-sample_panel_file integrated_call_samples_v3.20130502.ALL.panel \
-region 10:43077000-43131000 \
-population CHS \
-output_info RET.info \
-output_ped RET.ped
補充說明
參數(shù)解釋
-vcf [.vcf.gz文件]
-sample_panel_file [.panel文件]
-region [染色體:1234-12345]
-population [人群 CHS為中國南方]
可選參數(shù)
問題
.pl版本工具相比網(wǎng)頁那樣方便，網(wǎng)頁工具可以只保留Biallelic
如果選擇BRCA2這種基因，后續(xù)haploview會出現(xiàn)問題(haploview連鎖分析時默認為只有兩個等位基因)
不愧是遺傳腫瘤易感明星基因
