Bash - Commands used for running the analyses in the Workshop

Exploring the Genotype Data

In [12]:
SDIR=~/share;
cd $SDIR/genotype_data;
ls
genotypes.geno  genotypes.ind  genotypes.snp  outgroupf3stats.params.txt

Exploring the files. Here are the first 20 individuals:

In [2]:
head -20 genotypes.ind
             Yuk_009 M    Yukagir
             Yuk_025 F    Yukagir
             Yuk_022 F    Yukagir
             Yuk_020 F    Yukagir
               MC_40 M    Chukchi
             Yuk_024 F    Yukagir
             Yuk_023 F    Yukagir
               MC_16 M    Chukchi
               MC_15 F    Chukchi
               MC_18 M    Chukchi
             Yuk_004 M    Yukagir
               MC_08 F    Chukchi
             Nov_005 M   Nganasan
               MC_25 F    Chukchi
             Yuk_019 F    Yukagir
             Yuk_011 M    Yukagir
             Sesk_47 M   Chukchi1
               MC_17 M    Chukchi
             Yuk_021 M    Yukagir
               MC_06 F    Chukchi

And here the first 20 SNP rows:

In [3]:
head -20 genotypes.snp
            1_752566     1        0.020130          752566 G A
            1_842013     1        0.022518          842013 T G
            1_891021     1        0.024116          891021 G A
            1_903426     1        0.024457          903426 C T
            1_949654     1        0.025727          949654 A G
           1_1018704     1        0.026288         1018704 A G
           1_1045331     1        0.026665         1045331 G A
           1_1048955     1        0.026674         1048955 A G
           1_1061166     1        0.026711         1061166 T C
           1_1108637     1        0.028311         1108637 G A
           1_1120431     1        0.028916         1120431 G A
           1_1156131     1        0.029335         1156131 T C
           1_1157547     1        0.029356         1157547 T C
           1_1158277     1        0.029367         1158277 G A
           1_1161780     1        0.029391         1161780 C T
           1_1170587     1        0.029450         1170587 C T
           1_1205155     1        0.029735         1205155 A C
           1_1211292     1        0.029785         1211292 C T
           1_1235792     1        0.030045         1235792 C T
           1_1254255     1        0.030111         1254255 G A

And here are the first 20 genotypes of the first 100 individuals:

In [4]:
head -20 genotypes.geno | cut -c1-100
0101101211102210102021200100010200000011001000200001010110001100001111101001110200110111212110000011
2012121012210011122100111202201222121102222121121012202221211212202201101201220222122012011112122022
1100112001110021001001111000011200000111100001110001110100002100110111120000102200110112011020001001
0000112210222121221121100202221222122112112211202122222221022222111221102200112222122211012121022011
0000000000000000000000000000100000000000000000100010000000000000000000000000000100000011201120010000
1012100221102201101121110120110000010012002010200100010011100100011011101110120200010100000002010111
2222222222222222222222222222222222222222222222222222222222222222222222222222222222121222222222122222
2211222002212022102001212222212212222210122212121222112222221112122111222222122021221122222222222222
2211222002212022102001212202012212212210122212121122112221221112121111222122112021211121011121222211
2222222222102222202222222222222222222222222211222212122222122122222222222222222122221212212212222222
2212222212122222222222222222221222222222222220221122222222122221212222221222222202222210012292222222
1101100001000001001000000222010021200001202110101111110122100021211110001221120002110000202110121222
1221121221222211221222222121221222212222222222222222222211121221212122221202101222212222021012222222
1221121221222211221222222121221222212222222222222222222211121221212122221202101222212222021012222222
1221121221222211221222222121221222212222222222222222222211121221212122221202101222212222021012222222
2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222021012222222
2222222222222222222222222222222222222222222222222222222222222222222222222222222222121211000001222222
1011111102100111001100200122221022211211222021212200120222112121221120012221222102020101000001222212
2222222222222222222222222222222222222222222222222222222222222222222222222222222222222211000112222222
2122121212102221202222222222222221122212222192222211122222222112222222222122222122221222222222122222

Counting how many individuals and SNPs there are:

In [5]:
wc -l genotypes.ind
wc -l genotypes.snp
1351 genotypes.ind
593124 genotypes.snp

And now we check that the first row of the *.geno file indeed contains the same number of columns:

In [6]:
head -1 genotypes.geno | wc -c
1352

which is one more, including the newline character at the end of the line. Now counting the number of rows in the *.geno-file (this takes a few seconds, as the file is >2Gb):

In [7]:
wc -l genotypes.geno
593124 genotypes.geno

Great, the number of rows and columns agrees with the numbers indicated in the *.ind and *.snp file! Now we're counting how many different populations there are. Let's first see the first 10 populations in the sorted list, alongside the number of individuals in each group:

In [8]:
awk '{print $3}' genotypes.ind | sort | uniq -c | head -20
      9 Abkhasian
     16 Adygei
      6 Albanian
      7 Aleut
      4 Aleut_Tlingit
      7 Altaian
     10 Ami
     10 Armenian
      9 Atayal
     10 Balkar
     29 Basque
     25 BedouinA
     19 BedouinB
     10 Belarusian
      6 BolshoyOleniOstrov
      9 Borneo
     10 Bulgarian
      8 Cambodian
      2 Canary_Islander
      2 ChalmnyVarre

As you can see, there are a number of populations with only one sample. Those are typically ancient individuals which are considered individually in most analyses. Let's filter those out and count only populations with at least two individuals and count them:

In [9]:
awk '{print $3}' genotypes.ind | sort | uniq -c | awk '$1>1' | wc -l
114

OK, so there are 114 populations in this dataset.

Running smartPCA

We first have to create a parameter file with the input and output data.

In [73]:
WDIR=~/work/share/solutions
mkdir -p $WDIR
cd $WDIR
cat <<EOF > pca.WestEurasia.params.txt
genotypename: $SDIR/genotype_data/genotypes.geno
snpname: $SDIR/genotype_data/genotypes.snp
indivname: $SDIR/genotype_data/genotypes.ind
evecoutname: $WDIR/pca.WestEurasia.evec
evaloutname: $WDIR/pca.WestEurasia.eval
poplistname: $SDIR/WestEurasia.poplist.txt
lsqproject: YES
numoutevec: 4
numthreads: 1
EOF

Let's see whether it worked:

In [74]:
cat pca.WestEurasia.params.txt
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
evecoutname: /home/training/work/share/solutions/pca.WestEurasia.evec
evaloutname: /home/training/work/share/solutions/pca.WestEurasia.eval
poplistname: /home/training/share/WestEurasia.poplist.txt
lsqproject: YES
numoutevec: 4

Good, now we can run smartPCA:

In [ ]:
time smartpca -p pca.WestEurasia.params.txt > pca.WestEurasia.log.txt

This has run for 15 minutes. Now we'll run the AllEurasia one:

In [75]:
WDIR=~/work/share/solutions
mkdir -p $WDIR
cd $WDIR
cat <<EOF > pca.AllEurasia.params.txt
genotypename: $SDIR/genotype_data/genotypes.geno
snpname: $SDIR/genotype_data/genotypes.snp
indivname: $SDIR/genotype_data/genotypes.ind
evecoutname: $WDIR/pca.AllEurasia.evec
evaloutname: $WDIR/pca.AllEurasia.eval
poplistname: $SDIR/AllEurasia.poplist.txt
lsqproject: YES
numoutevec: 4
numthreads: 1
EOF
In [76]:
cat pca.AllEurasia.params.txt
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
evecoutname: /home/training/work/share/solutions/pca.AllEurasia.evec
evaloutname: /home/training/work/share/solutions/pca.AllEurasia.eval
poplistname: /home/training/share/AllEurasia.poplist.txt
lsqproject: YES
numoutevec: 4
In [ ]:
time smartpca -p pca.AllEurasia.params.txt

Running Admixture F3

Preparing the population file:

In [39]:
WDIR=~/work/share/solutions
mkdir -p $WDIR
cd $WDIR
cat <<EOF > $WDIR/f3.poplist.txt
Nganasan French Finnish 
Nganasan Icelandic Finnish 
Nganasan Lithuanian Finnish 
Nganasan Norwegian Finnish 
BolshoyOleniOstrov French Finnish 
BolshoyOleniOstrov Icelandic Finnish 
BolshoyOleniOstrov Lithuanian Finnish 
BolshoyOleniOstrov Norwegian Finnish
EOF
In [41]:
cat $WDIR/f3.poplist.txt
Nganasan French Finnish 
Nganasan Icelandic Finnish 
Nganasan Lithuanian Finnish 
Nganasan Norwegian Finnish 
BolshoyOleniOstrov French Finnish 
BolshoyOleniOstrov Icelandic Finnish 
BolshoyOleniOstrov Lithuanian Finnish 
BolshoyOleniOstrov Norwegian Finnish
In [42]:
WDIR=~/work/share/solutions
mkdir -p $WDIR
cd $WDIR
cat <<EOF > f3stats.params.txt
genotypename: $SDIR/genotype_data/genotypes.smaller.geno
snpname: $SDIR/genotype_data/genotypes.smaller.snp
indivname: $SDIR/genotype_data/genotypes.smaller.ind
popfilename: $WDIR/f3.poplist.txt
inbreed: YES
EOF
In [43]:
cat f3stats.params.txt
genotypename: /home/training/work/share/genotype_data/genotypes.smaller.geno
snpname: /home/training/work/share/genotype_data/genotypes.smaller.snp
indivname: /home/training/work/share/genotype_data/genotypes.smaller.ind
popfilename: /home/training/work/share/solutions/f3.poplist.txt
inbreed: YES
In [44]:
time qp3Pop -p f3stats.params.txt
qp3Pop: parameter file: f3stats.params.txt
### THE INPUT PARAMETERS
##PARAMETER NAME: VALUE
genotypename: /home/training/work/share/genotype_data/genotypes.smaller.geno
snpname: /home/training/work/share/genotype_data/genotypes.smaller.snp
indivname: /home/training/work/share/genotype_data/genotypes.smaller.ind
popfilename: /home/training/work/share/solutions/f3.poplist.txt
inbreed: YES
## qp3Pop version: 435
inbreed set
nplist: 8
number of blocks for block jackknife: 711
snps: 593124
                      Source 1             Source 2               Target           f_3       std. err           Z    SNPs
 result:              Nganasan               French              Finnish     -0.004539       0.000510      -8.894  442567
 result:              Nganasan            Icelandic              Finnish     -0.005297       0.000563      -9.404  427954
 result:              Nganasan           Lithuanian              Finnish     -0.005062       0.000590      -8.574  426231
 result:              Nganasan            Norwegian              Finnish     -0.004744       0.000569      -8.332  428161
 result:    BolshoyOleniOstrov               French              Finnish     -0.002814       0.000444      -6.341  402958
 result:    BolshoyOleniOstrov            Icelandic              Finnish     -0.002590       0.000486      -5.323  386418
 result:    BolshoyOleniOstrov           Lithuanian              Finnish     -0.001523       0.000536      -2.840  384134
 result:    BolshoyOleniOstrov            Norwegian              Finnish     -0.001553       0.000502      -3.092  386203
##end of qp3Pop:      190.789 seconds cpu      327.410 Mbytes in use

real	3m10.903s
user	3m10.156s
sys	0m0.664s

Running D stats

In [17]:
WDIR=~/work/share/solutions
cat <<EOF > $WDIR/dstat.poplist.txt
Mbuti Nganasan French Finnish 
Mbuti Nganasan Icelandic Finnish 
Mbuti Nganasan Lithuanian Finnish 
Mbuti Nganasan Norwegian Finnish 
Mbuti BolshoyOleniOstrov French Finnish 
Mbuti BolshoyOleniOstrov Icelandic Finnish 
Mbuti BolshoyOleniOstrov Lithuanian Finnish 
Mbuti BolshoyOleniOstrov Norwegian Finnish
EOF
In [18]:
cat $WDIR/dstat.poplist.txt
Mbuti Nganasan French Finnish 
Mbuti Nganasan Icelandic Finnish 
Mbuti Nganasan Lithuanian Finnish 
Mbuti Nganasan Norwegian Finnish 
Mbuti BolshoyOleniOstrov French Finnish 
Mbuti BolshoyOleniOstrov Icelandic Finnish 
Mbuti BolshoyOleniOstrov Lithuanian Finnish 
Mbuti BolshoyOleniOstrov Norwegian Finnish
In [19]:
cat <<EOF > $WDIR/dstats.params.txt
genotypename: $SDIR/genotype_data/genotypes.geno
snpname: $SDIR/genotype_data/genotypes.snp
indivname: $SDIR/genotype_data/genotypes.ind
popfilename: $WDIR/dstat.poplist.txt
f4mode: YES
EOF
In [20]:
cat $WDIR/dstats.params.txt
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
popfilename: /home/training/work/share/solutions/dstat.poplist.txt
f4mode: YES
In [21]:
time qpDstat -p $WDIR/dstats.params.txt
qpDstat: parameter file: /home/training/work/share/solutions/dstats.params.txt
### THE INPUT PARAMETERS
##PARAMETER NAME: VALUE
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
popfilename: /home/training/work/share/solutions/dstat.poplist.txt
f4mode: YES
## qpDstat version: 755
number of quadruples 8
  0                Mbuti   10
  1             Nganasan   11
  2   BolshoyOleniOstrov    6
  3               French   32
  4            Icelandic   12
  5           Lithuanian   10
  6            Norwegian   11
  7              Finnish    8
jackknife block size:     0.050
snps: 593124  indivs: 100
number of blocks for jackknife: 711
nrows, ncols: 100 593124
result:      Mbuti   Nganasan     French    Finnish      0.002363     19.016   29254  27852 593124 
result:      Mbuti   Nganasan  Icelandic    Finnish      0.001721     11.926   28915  27894 593124 
result:      Mbuti   Nganasan Lithuanian    Finnish      0.001368      9.664   28745  27933 593124 
result:      Mbuti   Nganasan  Norwegian    Finnish      0.001685     11.663   28933  27934 593124 
result:      Mbuti BolshoyOleniOstrov     French    Finnish      0.001962     16.737   27249  26175 547486 
result:      Mbuti BolshoyOleniOstrov  Icelandic    Finnish      0.001084      7.776   26876  26282 547486 
result:      Mbuti BolshoyOleniOstrov Lithuanian    Finnish      0.000554      3.942   26683  26380 547486 
result:      Mbuti BolshoyOleniOstrov  Norwegian    Finnish      0.000952      6.707   26873  26351 547486 
##end of qpDstat:      172.851 seconds cpu      731.099 Mbytes in use

real	2m52.891s
user	2m52.136s
sys	0m0.744s

Outgroup F3 stats

In [27]:
cat <<EOF > $WDIR/outgroupf3stats.params.txt
genotypename: $SDIR/genotype_data/genotypes.geno
snpname: $SDIR/genotype_data/genotypes.snp
indivname: $SDIR/genotype_data/genotypes.ind
popfilename: $WDIR/outgroupF3pops_Han.txt
EOF
In [28]:
cat $WDIR/outgroupf3stats.params.txt
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
popfilename: /home/training/work/share/solutions/outgroupF3pops_Han.txt
In [29]:
time qp3Pop -p $WDIR/outgroupf3stats.params.txt
qp3Pop: parameter file: /home/training/work/share/solutions/outgroupf3stats.params.txt
### THE INPUT PARAMETERS
##PARAMETER NAME: VALUE
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
popfilename: /home/training/work/share/solutions/outgroupF3pops_Han.txt
## qp3Pop version: 435
nplist: 32
number of blocks for block jackknife: 711
snps: 593124
                      Source 1             Source 2               Target           f_3       std. err           Z    SNPs
 result:                   Han              Chuvash                Mbuti      0.233652       0.002072     112.782  502678
 result:                   Han             Albanian                Mbuti      0.215629       0.002029     106.291  501734
 result:                   Han             Armenian                Mbuti      0.213724       0.001963     108.882  504370
 result:                   Han            Bulgarian                Mbuti      0.216193       0.001979     109.266  504310
 result:                   Han                Czech                Mbuti      0.218060       0.002002     108.939  504089
 result:                   Han                Druze                Mbuti      0.209551       0.001919     109.205  510853
 result:                   Han              English                Mbuti      0.216959       0.001973     109.954  504161
 result:                   Han             Estonian                Mbuti      0.220730       0.002019     109.332  503503
 result:                   Han              Finnish                Mbuti      0.223447       0.002044     109.345  502217
 result:                   Han               French                Mbuti      0.216623       0.001969     110.012  509613
 result:                   Han             Georgian                Mbuti      0.214295       0.001935     110.721  503598
 result:                   Han                Greek                Mbuti      0.215203       0.001984     108.465  507475
 result:                   Han            Hungarian                Mbuti      0.217894       0.001999     109.004  507409
 result:                   Han            Icelandic                Mbuti      0.218683       0.002015     108.553  504655
 result:                   Han        Italian_North                Mbuti      0.215332       0.001978     108.854  507589
 result:                   Han        Italian_South                Mbuti      0.211787       0.002271      93.265  492400
 result:                   Han           Lithuanian                Mbuti      0.219615       0.002032     108.098  503681
 result:                   Han              Maltese                Mbuti      0.210359       0.001956     107.542  503985
 result:                   Han            Mordovian                Mbuti      0.223469       0.002008     111.296  503441
 result:                   Han            Norwegian                Mbuti      0.218873       0.002023     108.197  504621
 result:                   Han             Orcadian                Mbuti      0.217773       0.002014     108.115  504993
 result:                   Han              Russian                Mbuti      0.223993       0.001995     112.274  506525
 result:                   Han            Sardinian                Mbuti      0.213230       0.001980     107.711  508413
 result:                   Han             Scottish                Mbuti      0.218489       0.002039     107.145  499784
 result:                   Han             Sicilian                Mbuti      0.212272       0.001975     107.486  505477
 result:                   Han        Spanish_North                Mbuti      0.215885       0.002029     106.383  500853
 result:                   Han              Spanish                Mbuti      0.213869       0.001975     108.297  513648
 result:                   Han            Ukrainian                Mbuti      0.218716       0.002007     108.950  503981
 result:                   Han           Levanluhta                Mbuti      0.236252       0.002383      99.123  263049
 result:                   Han   BolshoyOleniOstrov                Mbuti      0.247814       0.002177     113.849  457102
 result:                   Han         ChalmnyVarre                Mbuti      0.233499       0.002304     101.345  366220
 result:                   Han             Saami.DG                Mbuti      0.236198       0.002274     103.852  489038
##end of qp3Pop:      296.879 seconds cpu      327.410 Mbytes in use

real	4m56.909s
user	4m56.140s
sys	0m0.756s
In [30]:
cat <<EOF > $WDIR/outgroupf3stats.params.txt
genotypename: $SDIR/genotype_data/genotypes.geno
snpname: $SDIR/genotype_data/genotypes.snp
indivname: $SDIR/genotype_data/genotypes.ind
popfilename: $WDIR/outgroupF3pops_MA1.txt
EOF
In [31]:
cat $WDIR/outgroupf3stats.params.txt
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
popfilename: /home/training/work/share/solutions/outgroupF3pops_MA1.txt
In [32]:
time qp3Pop -p $WDIR/outgroupf3stats.params.txt
qp3Pop: parameter file: /home/training/work/share/solutions/outgroupf3stats.params.txt
### THE INPUT PARAMETERS
##PARAMETER NAME: VALUE
genotypename: /home/training/share/genotype_data/genotypes.geno
snpname: /home/training/share/genotype_data/genotypes.snp
indivname: /home/training/share/genotype_data/genotypes.ind
popfilename: /home/training/work/share/solutions/outgroupF3pops_MA1.txt
## qp3Pop version: 435
nplist: 32
number of blocks for block jackknife: 711
snps: 593124
                      Source 1             Source 2               Target           f_3       std. err           Z    SNPs
 result:             MA1_HG.SG              Chuvash                Mbuti      0.243818       0.002349     103.781  350484
 result:             MA1_HG.SG             Albanian                Mbuti      0.236494       0.002296     103.008  344332
 result:             MA1_HG.SG             Armenian                Mbuti      0.231399       0.002264     102.229  349612
 result:             MA1_HG.SG            Bulgarian                Mbuti      0.237498       0.002281     104.103  349800
 result:             MA1_HG.SG                Czech                Mbuti      0.243224       0.002328     104.457  349553
 result:             MA1_HG.SG                Druze                Mbuti      0.226740       0.002197     103.193  359004
 result:             MA1_HG.SG              English                Mbuti      0.243135       0.002317     104.941  349321
 result:             MA1_HG.SG             Estonian                Mbuti      0.247065       0.002362     104.619  348861
 result:             MA1_HG.SG              Finnish                Mbuti      0.245684       0.002379     103.266  347208
 result:             MA1_HG.SG               French                Mbuti      0.240235       0.002269     105.886  357842
 result:             MA1_HG.SG             Georgian                Mbuti      0.232645       0.002253     103.243  349082
 result:             MA1_HG.SG                Greek                Mbuti      0.236566       0.002280     103.757  355261
 result:             MA1_HG.SG            Hungarian                Mbuti      0.241720       0.002313     104.483  355340
 result:             MA1_HG.SG            Icelandic                Mbuti      0.244488       0.002386     102.481  350287
 result:             MA1_HG.SG        Italian_North                Mbuti      0.236407       0.002273     104.002  354999
 result:             MA1_HG.SG        Italian_South                Mbuti      0.230839       0.002767      83.427  321217
 result:             MA1_HG.SG           Lithuanian                Mbuti      0.246864       0.002403     102.718  348656
 result:             MA1_HG.SG              Maltese                Mbuti      0.230200       0.002259     101.903  347725
 result:             MA1_HG.SG            Mordovian                Mbuti      0.245284       0.002346     104.571  350058
 result:             MA1_HG.SG            Norwegian                Mbuti      0.243930       0.002301     106.031  350182
 result:             MA1_HG.SG             Orcadian                Mbuti      0.243614       0.002320     105.008  351053
 result:             MA1_HG.SG              Russian                Mbuti      0.245212       0.002298     106.698  355953
 result:             MA1_HG.SG            Sardinian                Mbuti      0.231967       0.002264     102.449  355548
 result:             MA1_HG.SG             Scottish                Mbuti      0.244598       0.002434     100.512  339441
 result:             MA1_HG.SG             Sicilian                Mbuti      0.231141       0.002260     102.297  351028
 result:             MA1_HG.SG        Spanish_North                Mbuti      0.238479       0.002426      98.319  341661
 result:             MA1_HG.SG              Spanish                Mbuti      0.235386       0.002257     104.293  361951
 result:             MA1_HG.SG            Ukrainian                Mbuti      0.243551       0.002345     103.881  348948
 result:             MA1_HG.SG           Levanluhta                Mbuti      0.247704       0.003055      81.090  172055
 result:             MA1_HG.SG   BolshoyOleniOstrov                Mbuti      0.256041       0.002624      97.561  305851
 result:             MA1_HG.SG         ChalmnyVarre                Mbuti      0.249619       0.002862      87.212  239594
 result:             MA1_HG.SG             Saami.DG                Mbuti      0.251530       0.002622      95.922  326072
##end of qp3Pop:      219.087 seconds cpu      327.410 Mbytes in use

real	3m39.156s
user	3m38.432s
sys	0m0.672s