There are several lines in Manhattan plot.

Question

[kgwkk2]

I generated Manhattan plot with multi-calling vcf file. But when I got the result, there were several weird lines that are not usually observed in other Manhattan plots. I want to know whether my code is wrong or my vcf file is wrong.

My code was like this.

plink --vcf Asp3.vcf.gz --allow-no-sex --recode --out myplink
plink --vcf Asp3.vcf.gz --allow-no-sex --recode --pheno MANN_pheno.txt --make-bed --out myplink
plink --bfile myplink --allow-no-sex --make-bed --mind 0.05 --maf 0.05 --geno 0.1 --hwe 1e-6 --recode --out myplink.QC
plink --bfile myplink.QC --allow-no-sex --linear --pheno MANN_pheno.txt --out myplink_result

With the result named 'myplink_result.assoc.linear', I made plot by R.

gwas_results <- read.table("myplink_result.assoc.linear", header=T, sep="\t")
gwas_results$P = as.numeric(as.character(gwas_results$P))
manhattan(gwas_results, chr="CHR", bp="BP", p="P", main="GWAS Manhattan Plot", ylim=c(0, -log10(5e-8)))

[Ram]

I assume you're using the qqman package - you should have mentioned which package you're using.

The vignette shows the horizontal lines as well as the vertical chromosome bands: https://cran.r-project.org/web/packages/qqman/vignettes/qqman.html

Can you show us an example of a plot where you don't see these lines?

Answer 1

[bk11]

Checking your lines of codes-

plink --vcf Asp3.vcf.gz --allow-no-sex --recode --out myplink

You don't have to run the above code as you are converting vcf into plink .bed/.bim/.fam file and updating phenotypes at the same time in the following code-

plink --vcf Asp3.vcf.gz --allow-no-sex --recode --pheno MANN_pheno.txt --make-bed --out myplink

In the following QC line of code, --recode is unnecessary.

plink --bfile myplink --allow-no-sex --make-bed --mind 0.05 --maf 0.05 --geno 0.1 --hwe 1e-6 --recode --out myplink.QC

Can you tell me why you have --pheno here again? If you have any covariates in your MANN_pheno.txt, you can adjust for covariates as follows-

plink --bfile myplink.QC --linear --covar MANN_pheno.txt --out myplink_result

instead of doing this-

plink --bfile myplink.QC --allow-no-sex --linear --pheno MANN_pheno.txt --out myplink_result

Can you also answer the following questions?

1. Do you know what type of trait are you testing in these data?
2. Were multiple datasets merged before you started with this vcf (Asp3.vcf.gz) data? Alleles flips/mismatch could have happen if the data were not correctly merged.
3. Are there any genetically related individuals in your data?

I am suspecting that you are getting the flat lines of p-values (similar p-values) all over the genome because of related individuals. You could run --genome and check if there are genetically related individuals in your data. I would recommend not to include related individuals in the association test.

[kgwkk2]

Thank you for your answer.

1. I wanted to visualize about 80 fungal genome's SNPs via 8 chromosomes and detect highest and various SNPs regions.
2. Yes. I made vcf file by merging multiple vcf files of each strain. I used GATK4 for vcf calling and merged with bcf
3. I couldn't get your word 'individuals'. I collected same species but different strain's genomes. So they are genetically related.

I want to know more about you said 'related individuals'. And '--genome' function must be talking about plink, right?

[bk11]

Yes --genome function is from Plink. Looks like you have genetically related data. That could be one of the reason for the several lines seen in your manhattan plot. I would also suggest you to check if alleles flip happened while merging your data.