Dinh 2011/NOTES/2011-7-20

From ZhangLabWiki
Jump to navigation Jump to search

"Population structure and eigenanalysis" [1][edit]

  • This paper was published in PloS Genetics in 2006 and it had been cited 423 times according to Google scholar
  • A smarter way of performing PCA was developed to find population structure in large genetic datasets
  • Software link - David Reich's laboratory: [2]
  • Downloaded software and put into /usr/local/bin (path)
  • Software works best with >100,000 markers
Command for running smartpca (PCA on input genotype data):
smartpca.perl -i example.geno -a example.snp -b example.ind -k 10 -o example.pca -p example.plot -e example.eval -l example.log -m 5 -t 2 -s 6.0

-i input genotype file
-a input snp file
-b input indiv file
-k number of principal components to output
-o output file of principal components
-p prefix for plot files of top 2 principal components
-e output file for eigenvalues
-l output logfile
-m maximum number of iterations for removing outliers
-t maximum number of outliers to remove at each iterations
-s number of stdev an individual must exceed along one of the top principal components to be removed as an outlier.
-q YES for QTL
*.geno file format:
1 line per SNP, and 1 character per individual:
0 means 0 copies of the reference allele (could use 0 for methylated)
1 means 1 copy of the reference allele (use 1 for >0 and <1 methylation)
2 means 2 copies of the reference allele (use 2 for unmethylated)
9 means missing data (use 9 for NA)
*.snp file format:
1 line per SNP, and 6 columns (last 2 optional)
1st is SNP name
2nd is chromosome. Encode X as 23, Y as 24, mtDNA as 90, XY as 91
3rd is genetic position (in Morgans). 0.0 if unknown
4th is physical positions (in bases)
5th is reference alleles 
6th is variant alleles
*.ind format:
1 line per individuals and 3 columns
1st is sample ID, limited to 39 characters
2nd is gender (M or F)
3rd is label as "Case" or "Control" or might be population group label. Can be set to "Ignore" to ignore this individual.

Ran this software on UPenn 48 samples, methylation matrix with minimum 42 samples covered, filtered out sites wherein technical replicates has >=0.1diff . (N = 398020 sites)

File:MethylMatrix2smartPCA.pl
File:RunSmartPCA.sh
File:SmartPCA.log

The steps below appears to not work well when we don't have case versus control stratification[edit]

Command for running smarteigenstrat (statification correction):
smarteigenstrat.perl -i example.geno -a example.snp -b example.ind -p example.pca -k 1 -o example.chisq -l example.log
-i input genotype file
-a input snp file
-b input indiv file
-k number of principal components to correct for stratification must be less than -k above
-o output file of chisq statistics
-p prefix for plot files of top 2 principal components
-l output logfile
-q YES for quantitative data
Command for running gc on quantitaive data (run genomic control):
qc.perl example.chisq example.chisq.GC