Error: No variants remaining after --extract

[Ana Marija]

Hello Chris,

I downloaded 1000G from here:

wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz

converted it to Plink format and did all required QC on that data and got 1kG_MDS5 files. Later I extracted the variants present in my dataset (output12) from the 1000 genomes dataset.(1kG_MDS5)

awk '{print$2}' output12.bim > HapMap_SNPs.txt
/projects/com_grassim/anamaria/anamaria/plink/plink --bfile 1kG_MDS5 --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6

but I got this error:

PLINK v1.90b6.12 64-bit (28 Oct 2019)
Options in effect:
  --bfile 1kG_MDS4
  --extract HapMap_SNPs.txt
  --make-bed
  --out 1kG_MDS6

Random number seed: 1585071880
128528 MB RAM detected; reserving 64264 MB for main workspace.
8240745 variants loaded from .bim file.
629 people (0 males, 0 females, 629 ambiguous) loaded from .fam.
Ambiguous sex IDs written to 1kG_MDS6.nosex .
Error: No variants remaining after --extract.

Looking at 1kG_MDS5.bim

22    rs62240043    0 51233312 G A
22    22:51235447[b37]A,G   0 51235447 A G
22    22:51235496[b37]A,G   0 51235496 G A
22    22:51235498[b37]A,G   0 51235498 G A
22    22:51235547[b37]A,C   0 51235547 C A
22    22:51235959[b37]C,T   0 51235959 C T

it seems to be about formatting of this 2nd column because

head  HapMap_SNPs.txt

1:785989:T:C
1:1130727:A:C
1:1156131:T:C
1:1158631:A:G

...

Is the solution for this to change format of the 2nd column in 1kG_MDS5.bim (remove [b37] and replace , with : ) so that it matches HapMap_SNPs.txt or I have to download a different 1000G dataset (and which one) or something else?

Thanks

Ana

[Ana Marija]

I tried reformatting 1kG_MDS5.bim and it looks like

1 rs58108140 0 10583 A G
1 1:11508:A:G 0 11508 A G
1 1:15820:G:T 0 15820 T G
1 1:16257:C:G 0 16257 C G

but still I got the same error after:

plink --bfile 1kG_MDS5 --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6

128528 MB RAM detected; reserving 64264 MB for main workspace.

5808310 variants loaded from .bim file.

629 people (0 males, 0 females, 629 ambiguous) loaded from .fam.
Ambiguous sex IDs written to 1kG_MDS6.nosex .
Error: No variants remaining after --extract.

Can you please advise?

[Christopher Chang]

https://www.biostars.org/p/306616/

On Tuesday, March 24, 2020 at 12:23:05 PM UTC-7, Ana Marija wrote:

I tried reformatting 1kG_MDS5.bim and it looks like

1 rs58108140 0 10583 A G
1 1:11508:A:G 0 11508 A G
1 1:15820:G:T 0 15820 T G
1 1:16257:C:G 0 16257 C G

but still I got the same error after:

plink --bfile 1kG_MDS5 --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6

128528 MB RAM detected; reserving 64264 MB for main workspace.
5808310 variants loaded from .bim file.
629 people (0 males, 0 females, 629 ambiguous) loaded from .fam.
Ambiguous sex IDs written to 1kG_MDS6.nosex .
Error: No variants remaining after --extract.

Can you please advise?

[Ana Marija]

just to confirm you are advising to use hapmap3_r1_b37_fwd_consensus.qc.poly.recode files downloaded from there instead of my 1000G downloaded and QCd data?

so do this:

awk '{print$2}' output12.bim > HapMap_SNPs.txt

/projects/com_grassim/anamaria/anamaria/plink/plink --bfile hapmap3_r1_b37_fwd_consensus.qc.poly.recode --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6

--
You received this message because you are subscribed to the Google Groups "plink2-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/d80f2d11-c26a-4751-a38e-ab02c3225d7d%40googlegroups.com.

[Christopher Chang]

No, you need to convert in the other direction.  The original HapMap files are very old and use hg18 coordinates.  I provided a link with HapMap SNPs converted to hg19 coordinates, which matches 1000G.

On Tuesday, March 24, 2020 at 12:41:18 PM UTC-7, Ana Marija wrote:

just to confirm you are advising to use hapmap3_r1_b37_fwd_consensus.qc.poly.recode files downloaded from there instead of my 1000G downloaded and QCd data?

so do this:

awk '{print$2}' output12.bim > HapMap_SNPs.txt

/projects/com_grassim/anamaria/anamaria/plink/plink --bfile hapmap3_r1_b37_fwd_consensus.qc.poly.recode --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users+unsubscribe@googlegroups.com.

[Ana Marija]

I am a bit confused about what I have to do. Let me just put this in a quick context:

I converted my Affymetrix output (CEL files) to VCF files and then converted that to Plink format using:

https://github.com/freeseek/gtc2vcf

http://apol1.blogspot.com/2014/11/best-practice-for-converting-vcf-files.html

I got my output files. Then I followed GWAS tutorial here: https://www.ncbi.nlm.nih.gov/pubmed/29484742

The last step in tutorial, which performed without issues, was removing related individuals:

plink --bfile output11 --filter-founders --make-bed --out output12

Now I am at the step of doing population stratification. At that was the point when I downloaded those 1000G...

Is there is a better tutorial on population stratification or can you please tell me exactly what do I do with the files you just shared with me?

what I was doing for preparing files for population stratification was this:

## Download 1000 Genomes data ##
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20100804/ALL.2of4intersection.20100804.genotypes.vcf.gz

# Convert vcf to Plink format.
plink --vcf ALL.2of4intersection.20100804.genotypes.vcf.gz --make-bed --out ALL.2of4intersection.20100804.genotypes

## QC on 1000 Genomes data.
# Remove variants based on missing genotype data.
plink --bfile ALL.2of4intersection.20100804.genotypes_no_missing_IDs --geno 0.2 --allow-no-sex --make-bed --out 1kG_MDS

# Remove individuals based on missing genotype data.
plink --bfile 1kG_MDS --mind 0.2 --allow-no-sex --make-bed --out 1kG_MDS2

# Remove variants based on missing genotype data.
plink --bfile 1kG_MDS2 --geno 0.02 --allow-no-sex --make-bed --out 1kG_MDS3

# Remove individuals based on missing genotype data.
plink --bfile 1kG_MDS3 --mind 0.02 --allow-no-sex --make-bed --out 1kG_MDS4

# Remove variants based on MAF.
plink --bfile 1kG_MDS4 --maf 0.05 --allow-no-sex --make-bed --out 1kG_MDS5
sed -i 's/\[b37\]\(.\),/:\1:/' 1kG_MDS5.bim

# Extract the variants present in my dataset from the 1000 genomes dataset. output12.bim is my output when related individuals are removed.
awk '{print$2}' output12.bim > HapMap_SNPs.txt

plink --bfile 1kG_MDS5 --extract HapMap_SNPs.txt --make-bed --out 1kG_MDS6

Error is in this last step.

Please advise

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/d80f2d11-c26a-4751-a38e-ab02c3225d7d%40googlegroups.com.

--

You received this message because you are subscribed to the Google Groups "plink2-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/8345b509-b20d-4ff8-97b6-ae5e42dc4e4b%40googlegroups.com.

[Christopher Chang]

Sorry, I was confused by the HapMap_SNPs.txt filename.  Ignore the first link I posted, then.

The main point is that your two datasets appear to be using different reference genome versions.  If this is true, you need to change the reference genome version of one of the datasets.

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users+unsubscribe@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/d80f2d11-c26a-4751-a38e-ab02c3225d7d%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "plink2-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users+unsubscribe@googlegroups.com.

[Christopher Chang]

The other possibility is that your REF/ALT order is exactly backwards.  Actually, you did get REF/ALT order wrong, but probably not in a way that would get *every* variant wrong since plink 1.x tries to change allele order to major/minor every time you run it, and the first plink command line you posted doesn't prevent it from doing so.  See https://www.cog-genomics.org/plink/1.9/data#ax_allele .  So you'll need to fix that too once you've gotten the reference genome versions synced up.  The simplest way is to use only plink 2.0 --pfile/--make-pgen instead of 1.x --bfile/--make-bed, at least until you've given all your variants their final names (which depend on REF/ALT being in the correct order).

[Ana Marija]

Thank you so much I will try that. I also mixed Plink2 and Plink in my workflow...I don't know if that would make any difference? (with respect to flipping order of alleles)

And to process my CEL files to Plink I was using human_g1k_v37.fasta as reference

and from 1000G which one of these you would advise using?

ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/

The last one I used was ALL.2of4intersection.20100804.genotypes.vcf.gz

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/d80f2d11-c26a-4751-a38e-ab02c3225d7d%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "plink2-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/8345b509-b20d-4ff8-97b6-ae5e42dc4e4b%40googlegroups.com.

--

You received this message because you are subscribed to the Google Groups "plink2-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/f0765a54-a3b8-42d5-b3b6-38421b36b992%40googlegroups.com.

[Ana Marija]

Hi Christopher,

can you please give me some thoughts on this last email.

Thanks

Ana

[Ana Marija]

Hi Chris,

can you please tell me how do I convert output.bed, output.bim, output.fam to --pgen/--pvar/--psam

Thanks

Ana

[Ana Marija]

to convert to plink format after Affymetix analysis I was using this:

bcftools norm -Ou -m -any gokind.bcf |
  bcftools norm -Ou -f human_g1k_v37.fasta |
  bcftools annotate -Ob -x ID \
    -I +'%CHROM:%POS:%REF:%ALT' |
  plink --bcf /dev/stdin \
    --keep-allele-order \
    --vcf-idspace-to _ \
    --const-fid \
    --allow-extra-chr 0 \
    --split-x b37 no-fail \
    --make-bed \
    --out output

is there is a way to rewrite this for Plink2 if not how do I change my output format to Plink2, pgen, psam...to proceed the analysis in Plink2

Thanks

Ana