Help needed: smartpca filters all SNPs

[CPerscheid]

Hello,

I need some help in understanding smartpca's output. Currently, it filters out all SNPs and I do not know why/what I need to change in the input data to avoid that.

These are my parameters:

genotypename: /path/to/pedfile.ped
snpname: /path/to/snpfile.pedsnp
indivname: /path/to/pedind.pedind
evecoutname: /path/to/test.pca.evec
evaloutname: /path/to/test.eval
altnormstyle: NO
numoutevec: 10
numoutlieriter: 0
numoutlierevec: 10
outliersigmathresh: 6
qtmode: 0

This is what is contained in the .log file:

altnormstyle: NO
numoutevec: 10
numoutlieriter: 0
numoutlierevec: 10
outliersigmathresh: 6
qtmode: 0
smartpca version: 16000
norm used

genetic distance set from physical distance
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 2300 230038
zzz 8617 39740
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 2624 262414
zzz 8618 39741
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 2862 286186
zzz 8619 39742
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 6698 669836
zzz 8620 39743
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 6756 675561
zzz 8621 39744
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 9096 909648
zzz 8622 39745
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 10835 1083519
zzz 8623 39746
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 10984 1098352
zzz 8624 39747
snp order check fail; snp list not ordered: /path/to/pedsnp.pedsnp (processing continues) 2 12092 1209219
zzz 8625 39748
genotype file processed
number of samples used: 974 number of snps used: 74808
Using 8 threads, and partial sum lookup algorithm.
snp snp1 ignored . allelecnt: 0 missing: 0
snp snp2 ignored . allelecnt: 0 missing: 0
snp snp3 ignored . allelecnt: 0 missing: 0
snp snp4 ignored . allelecnt: 0 missing: 0
snp snp5 ignored . allelecnt: 0 missing: 0
snp snp6 ignored . allelecnt: 0 missing: 0
snp snp7 ignored . allelecnt: 0 missing: 0
snp snp8 ignored . allelecnt: 0 missing: 0
snp snp9 ignored . allelecnt: 0 missing: 0
snp snp10 ignored . allelecnt: 0 missing: 0
total number of snps killed in pass: 74808 used: 0

This is how my pedsnp file looks like (snps are ordered by position):

1 snp1 0 494589 G A
1 snp2 0 633971 T C
1 snp3 0 633990 T C
1 snp4 0 956911 A G
1 snp5 0 956912 T A
1 snp6 0 961658 G A
1 snp7 0 962441 C T
1 snp8 0 964004 G A
1 snp9 0 966556 C T
1 snp10 0 970728 G T
1 snp11 0 979128 C T
1 snp12 0 979362 G A
1 snp13 0 999551 C G
...
X snp79094 0 155900547 C T
X snp79095 0 155939723 G T
X snp79096 0 155943159 G A
X snp79097 0 155997784 G T
X snp79098 0 156003498 C G
X snp79099 0 156007531 G C
X snp79100 0 156009838 C T
X snp79101 0 156010199 G A
X snp79102 0 156025453 C T

Any help appreciated. Thanks!

[kristjanmoore]

(Disclaimer: I'm not a maintainer of EIGENSOFT, I just use it a lot!)

smartpca expects that all chromosomes are coded as (small) integers, so the "X" chromosome here is likely messing some things up!

Additionally, you might want to just double-check that the variants really are sorted in numeric order by chromosome and position -- so chromosomes go "1, 2, 3..." rather than alphabetically-sorted "1, 10, 11...". (and make sure that you don't unsynchronise the .bed/.ped and .bim/.pedsnp files!)

The allelecnt: 0 lines are strange, but they could be downstream from one of the above issues. If not, try running something simple like plink --bfile yourData --missing --freq and verifying the output there looks alright (i.e. you don't have loads of SNPs with 0 MAF or 100% missingness).

[CPerscheid]

Oh somewhere read to encode chr23 this way...
However, you are right with your guess of alphabetical order. I will try this one first and then check on the other issue. It could actually be that I have many SNPs that occur only for one particular sample - these would be filtered out then?

[kristjanmoore]

As a rule of thumb you usually don't want to include variants with a really low minor allele frequency / count in a PCA, although if you want to look at very recent population history that can be useful. I think smartpca ignores variants that have a MAC of zero -- I don't think it "kills" variants if MAC is 1!

[CPerscheid]

I have run your plink command and it turns out frequency counts for all snps are zero... so they are probably filtered out (see output below)? Plink sets MAF to zero even if there is one occurrence of a SNP?
The data I am using are TCGA mutation annotation format (MAF) files (somatic variants), i.e. annotated VCF that I bring into the corresponding format required by smartpca. The samples there do not have so many mutations in general, so it is kind of expected that SNPs have a very low minor allele count.

CHR SNP A1 A2 MAF NCHROBS
1 snp1 0 G 0 1948
1 snp2 0 T 0 1948
1 snp3 0 T 0 1948
1 snp4 0 A 0 1948
1 snp5 0 T 0 1948
1 snp6 0 G 0 1948
1 snp7 0 C 0 1948
1 snp8 0 G 0 1948

[bumblenick]

A good way to diagnose what is going on is to run convertf with outputformat: eigenstrat This way you can see the genotypes and look if they are all zero. Also I recommend making allele counts/SNP. I'm sure PLINK will do that for you. SNPs that are not polymorphic are useless for PCA (and for most other purposes). Nick

…

On Thu, Oct 1, 2020 at 8:24 AM Cindy Perscheid ***@***.***> wrote: I have run your plink command and it turns out frequency counts for all snps are zero... so they are probably filtered out (see output below)? Plink sets MAF to zero even if there is one occurrence of a SNP? The data I am using are TCGA mutation annotation format (MAF) files (somatic variants), i.e. annotated VCF that I bring into the corresponding format required by smartpca. The samples there do not have so many mutations in general, so it is kind of expected that SNPs have a very low minor allele count. CHR SNP A1 A2 MAF NCHROBS 1 snp1 0 G 0 1948 1 snp2 0 T 0 1948 1 snp3 0 T 0 1948 1 snp4 0 A 0 1948 1 snp5 0 T 0 1948 1 snp6 0 G 0 1948 1 snp7 0 C 0 1948 1 snp8 0 G 0 1948 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#53 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEE77B4VN6KK5ESSBWQBA5TSIRYHLANCNFSM4R4TYUBA> .

[kristjanmoore]

@CPerscheid Your use case seems pretty unusual! It seems pretty likely to me that smartpca might be in some ways overkill or not well optimised for what you're trying to do.

In PLINK output, MAF zero really does mean zero! If these mutations are so rare that there aren't many non-zero datapoints, then you might consider writing the genotype files out to some kind of text format, and then reading it into R and running prcomp(), or something like that. Or filtering the VCF directly with some VCF-focused tool and reading it into some statistical computing software to run PCA on, although I can't offer much advice on working with VCFs.

If you really want to go down the VCF -> ped -> smartpca route, at the very least I think you should pass your fileset through plink --mac 1. (Or 2.) You might also want to check that the conversion from VCF to ped worked properly... There's a lot of moving parts here. Anyway, it does seem like the problem is with your data rather than with smartpca!

[CPerscheid]

@kristjanmoore you are probably right; the data just does not fit to the scope covered by smartpca. I wanted to compute genotype covariates (to then use these for eQTL computation) and stuck to a reference workflow that used smartpca. I will give prcomp a try and see if this helps.

Many thanks for your help here!
Best,
Cindy