Question: Question about ROH analysis by Plink 1.9
1
gravatar for 60343011s
5 months ago by
60343011s10
60343011s10 wrote:

Hi all,

I have recently tried to estimate runs of homozygosity (ROH) from my vcf file by using plink 1.9.

I ran following code to generate binary files that plink required:

plink --vcf myfile.vcf --make-bed --out out_name --no-sex --no-parents --no-fid --no-pheno --allow-extra-chr

This vcf file only contains one individual and about 3 million SNPs.

I used --allow-extra-chr here because I mapped my sequences to a drift genome.

Then, I used following code (with default parameters), trying to estimate ROH of my sample:

plink -bfile out_name --homozyg --allow-extra-chr

The result gave me 0 ROH, and only header produced in .hom file.

I also tried different parameters with different SNP windows and criterions, such like:

plink -bfile out_name --homozyg --homozyg-window-snp 50 --homozyg-snp 50 --homozyg-window-missing 3 --homozyg-kb 100 --homozyg-density 1000 --allow-extra-chr

However, all the results were the same, that showed :

PLINK v1.90b6.12 64-bit (28 Oct 2019)          www.cog-genomics.org/plink/1.9/
(C) 2005-2019 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to plink.log.
Options in effect:

  --allow-extra-chr
  --bfile out_name
  --homozyg
  --homozyg-density 1000
  --homozyg-kb 100
  --homozyg-snp 50
  --homozyg-window-missing 3
  --homozyg-window-snp 50

515905 MB RAM detected; reserving 257952 MB for main workspace.
3708761 variants loaded from .bim file.
1 person (0 males, 0 females, 1 ambiguous) loaded from .fam.
Ambiguous sex ID written to plink.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 1 founder and 0 nonfounders present.
Calculating allele frequencies... done.
3708761 variants and 1 person pass filters and QC.
Note: No phenotypes present.
--homozyg: Scan complete, found 0 ROH.

Results saved to plink.hom + plink.hom.indiv + plink.hom.summary .

Does anyone have idea why this happened to my files?

Will be grateful for any suggestions.

snp plink next-gen • 492 views
ADD COMMENTlink modified 5 months ago by chrchang5236.9k • written 5 months ago by 60343011s10
3
gravatar for chrchang523
5 months ago by
chrchang5236.9k
United States
chrchang5236.9k wrote:

This is unsurprising if your VCF contains only positions where your sample differs from the reference genome; practically all the ROHs will span regions which are excluded by your VCF. You may need to reconstruct your VCF in a way that includes homozygous-REF calls in an unbiased manner.

ADD COMMENTlink written 5 months ago by chrchang5236.9k

Yes, my VCF file contain only genotype 0/1 sites. I called those sites by

bcftools mpileup -a AD,DP -q 20 -f ref.fa input.bam | bcftools call -V indels -m -v -O v -o myfile.vcf

So I need to recall all the sites (0/0,0/1,1/1) by function like --gvcf under bcftools, right?

ADD REPLYlink modified 5 months ago • written 5 months ago by 60343011s10

Just to update the results, I've successfully used a VCF file that contain all sites for ROH analysis. Thanks!

ADD REPLYlink written 5 months ago by 60343011s10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1612 users visited in the last hour