[Bug]: 03g issue with cg07959070

[drsora]

Contact Details

sonja.rajic@tuni.fi

Scripts

03g-perform_positive_control.sh

What happened?

03g doesn't run fully, output is below I can't upload 03g because some files are missing, I uploaded a-d

How can the bug be reproduced?

No response

R version

4.4.0 (April, 2024)

Python version

None

Relevant log output

make manhattan and qq plots
Reading in /data/sonja/godmc/results/03/positive_control_untransformed_cg07959070.PHENO1.glm.linear.gz GWAS results

Expecting a large meQTL near 22:50053871
Lowest p-value within 1e+05 base pairs: 0.00127855
WARNING!
There doesn't appear to be a QTL for this positive control
Please upload this section and contact GoDMC analysts before continuing.

Generating QQ-plot without cis chromosome for/data/sonja/godmc/results/03/positive_control_untransformed_cg07959070 with lambda 0.943745114347533
Generating manhantten plot without cis chromosome /data/sonja/godmc/results/03/positive_control_untransformed_cg07959070
Error in manhattan(man_data, bp = names(man_data)[pos_column], chr = names(man_data)[chr_column],  : 
  #CHROM column should be numeric. Do you have 'X', 'Y', 'MT', etc? If so change to numbers and try again.
Calls: main -> manhattan
Execution halted

[drsora]

I ran this in my other cohort and get a different error message: make manhattan and qq plots Reading in /data/sonja/godmc/luric/results/03/positive_control_untransformed_cg07959070.PHENO1.glm.linear.gz GWAS results Error in main() : Wrong column specified for p-values Execution halted

[SiyiSEA]

Hi @drsora ,

I’m just back from holiday, so apologies for the slight delay.

Regarding the first error, it seems that the lowest p-value(0.00127855) within the control window is greater than the positive_control_threshold (0.001). The positive_control_threshold can be set in the parameter file depend on the sample size of your cohort. In addition, the Chromosome column may contain characters other than numbers from 1 to 23.

The second error suggests that the P-value is outside the expected range of (0,1).

Would you mind checking the columns of the /data/sonja/godmc/luric/results/03/positive_control_untransformed_cg07959070.PHENO1.glm.linear.gz file from your two corhorts?

Is column 12 the P-value column? If the P-value is out of the range (0, 1), this could be causing the error. Is column 1 the Chromosome column, and does it only contain values from 1-23? Is column 2 the Position column? Is column 3 the SNP column?

Thanks

[drsora]

Dear @SiyiSEA

I'm sorry for the very (very) late answer. I checked the data and it seems like I do indeed have an X in my chromosome column. I'm not exactly sure why though because it wasn't there in the original data. Do you know how to fix this? Thanks!

[ejh243]

This should be fixed in an update to 02a address ed PR #30 . But to avoid you have to go back and rerun those scripts, it should be an easy addition to the R script to convert X chr back to 23.

[ejh243]

@epzjlm Also wondering if we should update the version of the pipeline and ensure 2a is run on 1.1.0?

[epzjlm]

@ejh243 for script 04 chrX should be coded as 23

[ejh243]

@epzjlm which 04 script? Also can you clarify whether it needs to be added to a plink command or an Rscript?

[epzjlm]

It is 04c script where it maps against hrc ref panel.

You can use --output-chr 26 in plink to recode chromosome as numeric. This can be added in 02b.

Where does it recode to X?

[ejh243]

04c is based on hase? I am not familiar enough with Hase to make any changes to that.

I do see in 02b a plink command and yes I can add that flag --output-chr 26. In 02a we added --allow-extra-chr --chr-set 23 to prevent a loading error, which I think we might need again here as well?

[epzjlm]

What we can do is add also a command in 04a eg.

nX=grep ^X ${bfile}.bim | wc -l if [ "$nX" -gt "0" ] then ${plink2} --bfile ${bfile} --make-bed --output-chr 26 --out ${hase_dir_in} else cp ${bfile}.bim ${hase_dir_in} cp ${bfile}.fam ${hase_dir_in} cp ${bfile}.bed ${hase_dir_in} fi