BGI-shenzhen / LDBlockShow

LDBlockShow: a fast and convenient tool for visualizing linkage disequilibrium and haplotype blocks based on VCF files
MIT License
134 stars 40 forks source link

how to solve this error #37

Closed wangshuang2024 closed 2 weeks ago

wangshuang2024 commented 8 months ago

@hewm2008 hello, I have some trouble with this tool. when I use the following code, it always return to a error. Thank you for your help.

$ LDBlockShow -InVCF thal_qc_4.vcf.gz -Region 1:28400000:29400000 -OutPdf -OutPut chr1:28968982_C/A

open OUT File error: chr1:28968982_C/A.region.vcf

less -S chr1:28968982_C/A.region.vcf

fileformat=VCFv4.2

FILTER=

fileDate=20231221

source=PLINKv1.90

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

contig=

INFO=

FORMAT=

bcftools_viewVersion=1.17+htslib-1.17

bcftools_viewCommand=view -r chr1:28400000-29400000 -o chr1_28968982_C_A.vcf thal_qc_4.vcf.gz; Date=Sat Dec 23 21:47:30 2023

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ....

hewm2008 commented 8 months ago

you shoud remove the "/" in your output file. so you can try as :

chr1:28968982_C/A ---> chr1_28968982_C_A

$ LDBlockShow -InVCF thal_qc_4.vcf.gz -Region 1:28400000:29400000 -OutPdf -OutPut chr1_28968982_C_A

wangshuang2024 commented 8 months ago

you shoud remove the "/" in your output file. so you can try as :

chr1:28968982_C/A ---> chr1_28968982_C_A

$ LDBlockShow -InVCF thal_qc_4.vcf.gz -Region 1:28400000:29400000 -OutPdf -OutPut chr1_28968982_C_A

Thank you so much. I solve this problem.

wangshuang2024 commented 8 months ago

@hewm2008 hello, Weiming, I have another problem with my plot. I use "plink --ld "to calculate LD value (r2) between chr1:28968982_C/A and chr1:28968075_T/A, it showed r2=0.662337

while, I use LDBlockShow to plot LD plot. the chr1:28968075_T/A point showed as grey colour. Could you tell me the problem?how can I solve this problem? this site is not a multiple allele loci. Thanks for your warm heart.

1704548413169

hewm2008 commented 8 months ago

the default is : the R^2 of this size with the peak site (the red site ) .
chr1:28968982_C/A with [peak site] is grey chr1:28968075_T/A with [peak site ] is grey

you want chr1:28968982_C/A with chr1:28968075_T/A : red
you can try add para ShowLDSVG -TopSite chr1:28968982 ... ...

wangshuang2024 commented 8 months ago

the default is : the R^2 of this size with the peak site (the red site ) . chr1:28968982_C/A with [peak site] is grey chr1:28968075_T/A with [peak site ] is grey

you want chr1:28968982_C/A with chr1:28968075_T/A : red you can try add para ShowLDSVG -TopSite chr1:28968982 ... ...

I use -TopSite chr1:28968982 to plot the chr1:28968982_C/A site, the red point in the top is this site. What I wonder is, one of the grey site is chr1:28968075_T/A, I use plink --ld to calculate the LD relationship between these two site, r2=0.66, but in my plot it is grey.

hewm2008 commented 8 months ago

the R2 is the same with haplopview . see this https://github.com/BGI-shenzhen/LDBlockShow/issues/18 to see the D' R^2 with these two sites.

wangshuang2024 commented 8 months ago

the R2 is the same with haplopview . see this #18 to see the D' R^2 with these two sites.

老师,我还是不太理解。我用plink --ld算的r2是0.66,D‘=0.9。但是图中红色方形的chr1:28968982_C/A位点与灰色的chr1:28968075_T/A位点没显示有连锁,chr1:28968075_T/A位点(另一个灰色的位点是chr1:28968073位点)呈现灰色的状态,这是表示什么呢?下面LDblock里面也没有这个点chr1:28968075_T/A的信息。

hewm2008 commented 8 months ago

我猜可能是 你的gwas的位点和下面的vcf的位点不一致, 确保过滤条件都是一样的,即最终都是用到相同的位点.
解决方案1 修改过滤条件
-MAF Min minor allele frequency filter [0.05] -Miss Max ratio of miss allele filter [0.25] -HWE Exact test of Hardy-Weinberg Equilibrium for SNP Pvalue[0] -Het Max ratio of het allele filter [1.0]

解决方案2 从提取gwas里面对应的位点出来

wangshuang2024 commented 8 months ago

我猜可能是 你的gwas的位点和下面的vcf的位点不一致, 确保过滤条件都是一样的,即最终都是用到相同的位点. 解决方案1 修改过滤条件 -MAF Min minor allele frequency filter [0.05] -Miss Max ratio of miss allele filter [0.25] -HWE Exact test of Hardy-Weinberg Equilibrium for SNP Pvalue[0] -Het Max ratio of het allele filter [1.0]

解决方案2 从提取gwas里面对应的位点出来

我检查了vcf和gwas_result的结果。我这个vcf是用的qc之后的,也就是上gwas之前的质控数据,总位点的个数vcf和gwas_result是对的上的。我搞不太懂为什么有的会变成灰色,而且那些灰色的位点都不是多等位位点。我再看看吧。谢谢老师解答。

hewm2008 commented 8 months ago

把gwas只留 xx.site.gz 这个文件位点出来再重新试看