"""
run smartpca
This uses galaxy code developed by Dan to deal with
arbitrary output files using an html dataset with it's own
subdirectory containing the arbitrary files
We create that html file and add all the links we need
Note that we execute the smartpca.perl program in the output subdirectory
to avoid having to clear out the job directory after running
Code to convert linkage format ped files into eigenstratgeno format is left here
in case we decide to autoconvert
Added a plot in R with better labels than the default eigensoft plot december 26 2007
DOCUMENTATION OF smartpca program:
smartpca runs Principal Components Analysis on input genotype data and
outputs principal components (eigenvectors) and eigenvalues.
The method assumes that samples are unrelated. (However, a small number
of cryptically related individuals is usually not a problem in practice
as they will typically be discarded as outliers.)
5 different input formats are supported. See ../CONVERTF/README
for documentation on using the convertf program to convert between formats.
The syntax of smartpca is "../bin/smartpca -p parfile". We illustrate
how parfile works via a toy example (see example.perl in this directory).
This example takes input in EIGENSTRAT format. The syntax of how to take input
in other formats is analogous to the convertf program, see ../CONVERTF/README.
The smartpca program prints various statistics to standard output.
To redirect this information to a file, change the above syntax to
"../bin/smartpca -p parfile >logfile". For a description of these
statistics, see the documentation file smartpca.info in this directory.
Estimated running time of the smartpca program is
2.5e-12 * nSNP * NSAMPLES^2 hours if not removing outliers.
2.5e-12 * nSNP * NSAMPLES^2 hours * (1+m) if m outlier removal iterations.
Thus, under the default of up to 5 outlier removal iterations, running time is
up to 1.5e-11 * nSNP * NSAMPLES^2 hours.
------------------------------------------------------------------------
DESCRIPTION OF EACH PARAMETER in parfile for smartpca:
genotypename: input genotype file (in any format: see ../CONVERTF/README)
snpname: input snp file (in any format: see ../CONVERTF/README)
indivname: input indiv file (in any format: see ../CONVERTF/README)
evecoutname: output file of eigenvectors. See numoutevec parameter below.
evaloutname: output file of all eigenvalues
OPTIONAL PARAMETERS:
numoutevec: number of eigenvectors to output. Default is 10.
numoutlieriter: maximum number of outlier removal iterations.
Default is 5. To turn off outlier removal, set this parameter to 0.
numoutlierevec: number of principal components along which to
remove outliers during each outlier removal iteration. Default is 10.
outliersigmathresh: number of standard deviations which an individual must
exceed, along one of the top (numoutlierevec) principal components, in
order for that individual to be removed as an outlier. Default is 6.0.
outlieroutname: output logfile of outlier individuals removed. If not specified,
smartpca will print this information to stdout, which is the default.
usenorm: Whether to normalize each SNP by a quantity related to allele freq.
Default is YES. (When analyzing microsatellite data, should be set to NO.
See Patterson et al. 2006.)
altnormstyle: Affects very subtle details in normalization formula.
Default is YES (normalization formulas of Patterson et al. 2006)
To match EIGENSTRAT (normalization formulas of Price et al. 2006), set to NO.
missingmode: If set to YES, then instead of doing PCA on # reference alleles,
do PCA on whether each data point is missing or nonmissing. Default is NO.
May be useful for detecting informative missingness (Clayton et al. 2005).
nsnpldregress: If set to a positive integer, then LD correction is turned on,
and input to PCA will be the residual of a regression involving that many
previous SNPs, according to physical location. See Patterson et al. 2006.
Default is 0 (no LD correction). If desiring LD correction, we recommend 2.
maxdistldregress: If doing LD correction, this is the maximum genetic distance
(in Morgans) for previous SNPs used in LD correction. Default is no maximum.
poplistname: If wishing to infer eigenvectors using only individuals from a
subset of populations, and then project individuals from all populations
onto those eigenvectors, this input file contains a list of population names,
one population name per line, which will be used to infer eigenvectors.
It is assumed that the population of each individual is specified in the
indiv file. Default is to use individuals from all populations.
phylipoutname: output file containing an fst matrix which can be used as input
to programs in the PHYLIP package, such as the "fitch" program for
constructing phylogenetic trees.
noxdata: if set to YES, all SNPs on X chr are excluded from the data set.
The smartpca default for this parameter is YES, since different variances
for males vs. females on X chr may confound PCA analysis.
nomalexhet: if set to YES, any het genotypes on X chr for males are changed
to missing data. The smartpca default for this parameter is YES.
badsnpname: specifies a list of SNPs which should be excluded from the data set.
Same format as example.snp. Cannot be used if input is in
PACKEDPED or PACKEDANCESTRYMAP format.
popsizelimit: If set to a positive integer, the result is that only the first
popsizelimit individuals from each population will be included in the
analysis. It is assumed that the population of each individual is specified
in the indiv file. Default is to use all individuals in the analysis.
The next 5 optional parameters allow the user to output genotype, snp and
indiv files which will be identical to the input files except that:
Any individuals set to Ignore in the input indiv file will be
removed from the data set (see ../CONVERTF/README)
Any data excluded or set to missing based on noxdata, nomalexhet and
badsnpname parameters (see above) will be removed from the data set.
The user may decide to output these files in any format.
outputformat: ANCESTRYMAP, EIGENSTRAT, PED, PACKEDPED or PACKEDANCESTRYMAP
genotypeoutname: output genotype file
snpoutname: output snp file
indivoutname: output indiv file
outputgroup: see documentation in ../CONVERTF/README
"""
import sys,os,time,subprocess,string,glob
from rgutils import RRun, galhtmlprefix, galhtmlpostfix, timenow, smartpca, rexe, plinke
verbose = False
def makePlot(eigpca='test.pca',title='test',pdfname='test.pdf',h=8,w=10,nfp=None,rexe=''):
"""
the eigenvec file has a # row with the eigenvectors, then subject ids, eigenvecs and lastly
the subject class
Rpy not being used here. Write a real R script and run it. Sadly, this means putting numbers
somewhere - like in the code as monster R vector constructor c(99.3,2.14) strings
At least you have the data and the analysis in one single place. Highly reproducible little
piece of research.
"""
debug=False
f = file(eigpca,'r')
R = []
if debug:
R.append('sessionInfo()')
R.append("print('dir()=:')")
R.append('dir()')
R.append("print('pdfname=%s')" % pdfname)
gvec = []
pca1 = []
pca2 = []
groups = {}
glist = [] # list for legend
ngroup = 1 # increment for each new group encountered for pch vector
for n,row in enumerate(f):
if n > 1:
rowlist = row.strip().split()
group = rowlist[-1]
v1 = rowlist[1]
v2 = rowlist[2]
try:
v1 = float(v1)
except:
v1 = 0.0
try:
v2 = float(v2)
except:
v2 = 0.0
if not groups.get(group,None):
groups[group] = ngroup
glist.append(group)
ngroup += 1 # for next group
gvec.append(groups[group]) # lookup group number
pca1.append('%f' % v1)
pca2.append('%f' % v2)
# now have vectors of group,pca1 and pca2
llist = [x.encode('ascii') for x in glist] # remove label unicode - eesh
llist = ['"%s"' % x for x in llist] # need to quote for R
R.append('llist=c(%s)' % ','.join(llist))
plist = range(2,len(llist)+2) # pch - avoid black circles
R.append('glist=c(%s)' % ','.join(['%d' % x for x in plist]))
pgvec = ['%d' % (plist[i-1]) for i in gvec] # plot symbol/colour for each point
R.append("par(lab=c(10,10,10))") # so our grid is denser than the default 5
R.append("par(mai=c(1,1,1,0.5))")
maint = title
R.append('pdf("%s",h=%d,w=%d)' % (pdfname,h,w))
R.append("par(lab=c(10,10,10))")
R.append('pca1 = c(%s)' % ','.join(pca1))
R.append('pca2 = c(%s)' % ','.join(pca2))
R.append('pgvec = c(%s)' % ','.join(pgvec))
s = "plot(pca1,pca2,type='p',main='%s', ylab='Second ancestry eigenvector'," % maint
s += "xlab='First ancestry eigenvector',col=pgvec,cex=0.8,pch=pgvec)"
R.append(s)
R.append('legend("top",legend=llist,pch=glist,col=glist,title="Sample")')
R.append('grid(nx = 10, ny = 10, col = "lightgray", lty = "dotted")')
R.append('dev.off()')
R.append('png("%s.png",h=%d,w=%d,units="in",res=72)' % (pdfname,h,w))
s = "plot(pca1,pca2,type='p',main='%s', ylab='Second ancestry eigenvector'," % maint
s += "xlab='First ancestry eigenvector',col=pgvec,cex=0.8,pch=pgvec)"
R.append(s)
R.append('legend("top",legend=llist,pch=glist,col=glist,title="Sample")')
R.append('grid(nx = 10, ny = 10, col = "lightgray", lty = "dotted")')
R.append('dev.off()')
rlog,flist = RRun(rcmd=R,title=title,outdir=nfp)
print >> sys.stdout, '\n'.join(R)
print >> sys.stdout, rlog
def getfSize(fpath,outpath):
"""
format a nice file size string
"""
size = ''
fp = os.path.join(outpath,fpath)
if os.path.isfile(fp):
n = float(os.path.getsize(fp))
if n > 2**20:
size = ' (%1.1f MB)' % (n/2**20)
elif n > 2**10:
size = ' (%1.1f KB)' % (n/2**10)
elif n > 0:
size = ' (%d B)' % (int(n))
return size
def runEigen():
""" run the smartpca prog - documentation follows
smartpca.perl -i fakeped_100.eigenstratgeno -a fakeped_100.map -b fakeped_100.ind -p fakeped_100 -e fakeped_100.eigenvals -l
fakeped_100.eigenlog -o fakeped_100.eigenout
DOCUMENTATION OF smartpca.perl program:
This program calls the smartpca program (see ../POPGEN/README).
For this to work, the bin directory containing smartpca MUST be in your path.
See ./example.perl for a toy example.
../bin/smartpca.perl
-i example.geno : genotype file in EIGENSTRAT format (see ../CONVERTF/README)
-a example.snp : snp file (see ../CONVERTF/README)
-b example.ind : indiv file (see ../CONVERTF/README)
-k k : (Default is 10) number of principal components to output
-o example.pca : output file of principal components. Individuals removed
as outliers will have all values set to 0.0 in this file.
-p example.plot : prefix of output plot files of top 2 principal components.
(labeling individuals according to labels in indiv file)
-e example.eval : output file of all eigenvalues
-l example.log : output logfile
-m maxiter : (Default is 5) maximum number of outlier removal iterations.
To turn off outlier removal, set -m 0.
-t topk : (Default is 10) number of principal components along which
to remove outliers during each outlier removal iteration.
-s sigma : (Default is 6.0) number of standard deviations which an
individual must exceed, along one of topk top principal
components, in order to be removed as an outlier.
now uses https://www.bx.psu.edu/cgi-bin/trac.cgi/galaxy/changeset/1832
All files can be viewed however, by making links in the primary (HTML) history item like:
Some Text
\n' % (progname,timenow())
lf.write('
\n') lf.write(' |
%s