In the meanwhile, I managed to track down one of the problems that leads to a difference between the missingness filter in Plink and GWASTools. There is quite a number of SNPs on the X-chromosome in my data set, which have theta values for males that are close to 0.5 except of 1.0, which implies heterozygosity where it shouldn't exist.

Example:

       Theta_1 Theta_4 Theta_8 Theta_9 Theta_10 Theta_11 Theta_12 Theta_13 Theta_14 Theta_15 Theta_16 Theta_19 Theta_20 Theta_22 Theta_25 Theta_26 Theta_27 Theta_28

0.08    0.06   0.065   0.701    0.059    0.694    0.065    0.703    0.068    0.059    0.698    0.684    0.666    0.692    0.069    0.681    0.665    0.677

Of course, the clustering leads to the following genotypes: AA, AA, AA, AB, AA and so on, where it should be AA, AA, AA, BB, AA ...

Plink automatically sets heterozygous genotypes to missing for male x-chromosomal data, which increases the number of filtered SNPs for a specific missing rate considerably. GWASTools seems to just keep the obviously wrong genotypes, or do I miss anything? This leads a the difference in filtered SNPs between the two programs. It would be very helpful if these erroneous SNPs could not only be tracked down, but could also be corrected.

Concerning the association results, I am sure that the used samples are the same and even if the number of SNPs would be different (due to different filters maybe), this shouldn't make a difference for the result of a specific SNP, since association tests are independent for each SNP, are they? I'd rather think the problem are different default options or something like that. In both cases I used very simple commands:

GWASTools

assoc.list <- lapply(unique(chr), function(x) {

  covar <- NULL
  if(x != 25) covar <- c("sex")
  start <- which(chr == x)[1]
  end <-which(chr==x)[length(which(chr==x))]
  assocRegression(genoData, outcome="cd34", covar=covar, model.type="linear",
                  snpStart=start, snpEnd=end)
})

Plink:

plink --bfile q_all --geno 0.15 --maf 0.06 --make-bed --nonfounders --out q_filt

plink --adjust --assoc --bfile q_filt --out quant.filt --qt-means

q_all are the files I exported from GWASTools, including my quantitative trait variable.

Now I found another problem, which causes the differences between the missing rate in Plink and GWASTools. I generated the input ped file for Plink with Bioconductor:

plinkWrite(genoData, pedFile="Plink", individual.col="scanID",
           phenotype.col="cd34", rs.col="rsID", mapdist.col=NULL)

For some SNPs, single data points are missing in the resulting ped-file.

Here is an example:

The genotypes of the SNP I obtain with getGenotype are: 0 0 1 0 0 0 1 0 1 0 0 0 0 2 0 1 1 1 0 0 0 0 0 0 0 0 1 0

which should come from the raw data:

B,B
B,B
A,B
B,B
B,B
B,B
A,B
B,B
A,B
B,B
B,B
B,B
B,B
A,A
B,B
A,B
A,B
A,B
B,B
B,B
B,B
B,B
B,B
B,B
B,B
B,B
A,B
B,B

When I check the ped file, I get:

A A
A A
A B
0 0
A A
A A
A B
A A
A B
A A
A A
A A
A A
B B
A A
A B
A B
A B
A A
A A
A A
A A
A A
A A
A A
A A
A B
A A

The data is suddenly missing in the forth line. I found SNPs for which half of the data was missing, which explains the higher missing rate in PLINK.