#!/bin/sh
# This example is on simulated data, indicative of the split between africa, europe and asia
# It was generated using the following commands, using the tool "msms" available from :

#> msms -N 10000 -t 5000 -r 4400 100000000 -ms 70 1 -I 3 30 20 20 -en 0.2 3 0.1 -ej 0.25 3 2 -en 0.25 2 0.1 -em 0.25 1
2 50 -ej 0.4 2 1 > example.txt
#> msms2cp.pl -n 10000000 example.txt example_cp
#> makeuniformrecfile.pl example_cp.phase example_cp.recombfile

# We have 15 African individuals (IND1-15) 10 Europeans (Inds16-25) and 10 Asians (Inds 26-35)

## This is how we process this in finestructure:
fs example_cp.cp -n -phasefiles example_cp.phase -recombfiles example_cp.recombfile -idfile example_cp.ids -s1minsnps
5000 -s3iters 10000 -s4iters 10000 -go
# Takes about 4 mins on a 2nd gen Intel i5 dual core laptop.

# See the help for those commands for details:
#> fs -h s1minsnps
#Help for Parameter s1minsnps : Minimum number of SNPs for EM estimation (for chromopainter -e, default: 10000)
#> fs -h s3iters
#Help for Parameter s3iters : Number of TOTAL iterations to use for MCMC. By default we assign half to burnin and half
to sampling. (default: 100000)
#$ fs -h s4iters
#Help for Parameter s4iters : Number of maximization steps when finding the best state from which the tree is built.
(default: 10000)

# follow the advice and load these results into the GUI, if you have installed it:
#> finegui -c example_cp_linked.chunkcounts.out -m example_cp_linked_x5000_y5000_z10_mcmc_run0.xml -t
example_cp_linked_x5000_y5000_z10_tree_run0.xml -m2 example_cp_linked_x5000_y5000_z10_mcmc_run1.xml -t2
example_cp_linked_x5000_y5000_z10_tree_run1.xml
# There isn't a lot to look at - the model is certain about the clustering, and it is right. We'll explore more in
example2.

13.2 Involved example with few SNPs and individuals

Listing 27: example2: complex use case on a larger dataset split by chromosome.

> cat examples/example2.sh

#!/bin/bash
### Example using the unlinked model on (nearly) unlinked data
## We are showing two different methods here.

###############################################
## First method: defining all parameters via a settings file, and performing the calculations in one go
# This is only recommended for small datasets
time fs example2a.cp -v -import example2.settings -go
# Takes about 8min 45s on a 2nd gen core i5 laptop (dual core, 4 cores with hyperthreading)
# You can specify the number of cores with e.g. -numthreads 4, but fs will use all of them by default.
# The command has run everything, including 2 independent finestructure runs!

# To examine the output, we will follow the suggested GUI output:
# > finegui -c example2a_unlinked.chunkcounts.out -m example2a_unlinked_mcmc_run0.xml -t
example2a_unlinked_tree_run0.xml
# Click File->Open
# Click "Read data file", "Read Pairwise coincidence", "Read Tree"
# Click "Done"
# Click File->Manage Second Dataset
# Change data filename from "run0" to "run1", click "Read data file"
# Change MCMC filename from "run0" to "run1", click "Read Pairwise coincidence"
# Change Tree filename from "run0" to "run1", click "Read Tree"
# Click "Done"
# View->Pairwise coincidence
# Second View->Enable alternative second view
# Second View->Use second dataset
# Second View->Pairwise coincidence
## Now run0 is displayed in the bottom left and run1 in the top right. They should be almost identical

## If you want to access the MCMC traces, they are stored in the file:
# example2a/example2a_linked_x30000_y30000_z60_mcmc.mcmctraces.tab
## The Gelman-Rubin convergence diagnostics are stored in the settings file:
#> grep "mcmcGR" example2a.cp
# mcmcGR:1.12997,1.0143,1.00456,0.999354,1.00513  # Gelman-Rubin diagnostics obtained from combining MCMC runs, for
log-posterior, K,log-beta,delta,f respectively

###############################################
## Second method: Using HPC mode, and defining the parameters inline. This achieves an identical result to the previous
approach, but allows explorting computation to another machine.
# NOTE: You need gnu "parallel" ("sudo apt-get install parallel" on Ubuntu) to be installed on your system ONLY for
this example. In a HPC setup you would replace those lines with calls to qsub_run.sh to submit the jobs to the cluster.
date > example2b.time
fs example2b.cp -hpc 1 -idfile Europe.small.ids -phasefiles EuropeSample.small.chrom{1..22}.phase -recombfiles
EuropeSample.small.chrom{1..22}.recombfile -s3iters 100000 -s4iters 50000 -s1minsnps 1000 -s1indfrac 0.1 -go
cat example2b_commandfile1.txt | parallel
fs example2b.cp -go
cat example2b_commandfile2.txt | parallel
fs example2b.cp -go
cat example2b_commandfile3.txt | parallel
fs example2b.cp -go
cat example2b_commandfile4.txt | parallel
fs example2b.cp -go
date >> example2b.time
# The output should be identical, to within MCMC tolerance

###############################################
## This is how we might define things in "best practice", i.e. by loading all non-input settings from a file:
# time fs example2c.cp -idfile Europe.small.ids -phasefiles EuropeSample.small.chrom{1..22}.phase -recombfiles
EuropeSample.small.chrom{1..22}.recombfile -import example2.settings -go

###############################################
# This is how we can continue previous runs.  We will test the tolerance of the model to deviations in the parameter
'c'.

## First, we don't know what the parameter might be called.
#> fs -h parameters | grep "'c'"
#Help for Parameter cval : 'c' as inferred using chromopainter. This is only used for sanity checking. See s34 args for
setting it manually.
#Help for Parameter s34args : Additional arguments to both finestructure mcmc and tree steps. Add "-c <val>" to
manually override 'c'.

## Now we've learned that we should set -s34args. But we've already run to stage 4, so we can't do this without
rerunning chromopainter. That would be madness... Lets instead duplicate the settings as of the beginning.  How do we
do that?
#> fs -h duplicate
#Help for Action    -duplicate <stage> <newfile.cp> : Copy the information from the current settings file, and then
-reset it to <stage>.

# What was inferred from the data?
#> grep "cval" example2a.cp
#cval:0.292178  # 'c' as inferred using chromopainter. This is only used for sanity checking. See s34 args for setting
it manually.

## OK, lets rerun with a chosen value of c.
fs example2a.cp -duplicate 3 example2a_c1.0.cp -s34args:-c\ 1.0 -go
fs example2a.cp -duplicate 3 example2a_c2.0.cp -s34args:-c\ 2.0 -go
fs example2a.cp -duplicate 3 example2a_c0.1.cp -s34args:-c\ 0.1 -go
fs example2a.cp -duplicate 3 example2a_c0.5.cp -s34args:-c\ 0.5 -go

#You can compare the pairwise coincidences with commands such as
# finegui -c example2a_linked.chunkcounts.out -m example2a_linked_x30000_y30000_z60_mcmc_run0.xml -t
example2a_linked_x30000_y30000_z60_tree_run0.xml -m2 example2a_linked_c1.0_x30000_y30000_z60_mcmc_run1.xml -t2
example2a_linked_c1.0_x30000_y30000_z60_tree_run1.xml &
# which shows that much of the substructure within the Orcadians is not clear at c=1.0 (though two pairs are
significantly different to the rest). At c=2.0, the main Orcadian group and Italian group might be merged, even though
the two pairs of Orcadians are not yet merged with the rest of them.
# At c=0.1 the Adygei are splitting up
# At c=0.5 looks mostly like c=1.0.
# Looking at the raw copy matrix, the structure claimed by the empirical value of 'c' does appear to be justified, and
that claimed by c=0.1 is not clear by eye.

# Similarly, if you wanted samples from chromopainter then you can rerun stage2 using
#fs example2a.cp -duplicate 2 example2a_samples.cp -s34args:-c\ 1.0 -go

##############################################
## Lets do an UNLINKED analysis
time fs example2c.cp -v -phasefiles EuropeSample.small.chrom{1..22}.phase -idfile Europe.small.ids -import
example2_bestpractice.settings -go
# Although the chunkcounts are on a different scale in an unlinked analysis, these SNPs have been thinned and therefore
the pairwise coincidence is almost identical (with just a little less evidence of substructure within the Orcadians,
but it is the same structure.)

##############################################
## An analysis of the linked data using R
ln -s ../../FinestructureLibrary.R
Rscript example2.R ##Take a look at how this works!

13.3 HGDP example, including downloading and processing.

This example gets the HGDP data from the website and processes it all. It is recommended only for HPC users, although with Chromosomes 21-22 it is fast enough for a single high specced machine.

Listing 28: example3: HGDP download and processing.

> cat examples/example_hgdp.sh

#!/bin/sh
## This example is intended to be run on a HPC machine. In particular, it will only work verbatim on a torque-based
qsub system.  Naturally, you can adapt it to run on other systems. The only part to change is "qsub_run.sh"

####### Download the sample ID info
wget -nc http://hgdp.uchicago.edu/Phased_data/H938_Clumps.clst.txt_check.Continent.format.order.NR.gz
zcat H938_Clumps.clst.txt_check.Continent.format.order.NR.gz  | cut -f2 -d' ' | awk 'NR%2==0{printf("%s%i %s
1\n",$1,NR/2,$1);}' > hgdp.ids

####### Download the genetic map
wget -nc http://hapmap.ncbi.nlm.nih.gov/downloads/recombination/2011-01_phaseII_B37/genetic_map_HapMapII_GRCh37.tar.gz
tar -xzvf genetic_map_HapMapII_GRCh37.tar.gz

####### Create scripts that download the raw data and convert it into chromopainter input format
mkdir -p getdata
cmdfile="getdata_raw.sh"
rm -f $cmdfile
for chr in ‘seq 21 22‘; do
    cmdf="getdata/getdata_chr$chr.sh"
    echo "#!/bin/sh" > $cmdf
    echo "wget -nc http://hgdp.uchicago.edu/Phased_data/chrom${chr}_hapguess_switch.out.gz" >> $cmdf
    echo "wget -nc http://hgdp.uchicago.edu/Phased_data/chrom${chr}.final.gz" >> $cmdf
    echo "fastphase2cp.sh chrom${chr}_hapguess_switch.out.gz chrom${chr}.final.gz chrom${chr}.phase" >> $cmdf
    echo "convertrecfile.pl -M hapmap chrom${chr}.phase genetic_map_GRCh37_chr${chr}.txt chrom${chr}.recombfile" >>
$cmdf
    chmod +x $cmdf
    echo "$cmdf" >> $cmdfile
done

## If you want the list of RSIDs
zcat chrom⋆.final.gz | tail -n +2 | cut -f1 > hgdprsids.txt

######## Submit these scripts to the qsub system
qsub_run.sh -n 8 -f getdata_raw.sh -w 1
## (nb: you could use cat getdata_raw.sh | parallel if you are using a single machine with many cores)

## Now we start the finestucture run proper.  First we do the painting.
fs hgdp_2chr.cp -n -phasefiles chrom21.phase chrom22.phase -recombfiles chrom21.recombfile chrom22.recombfile -idfile
hgdp.ids -hpc 1 -s1indfrac 0.1 -go # We are just using chr21-22.
qsub_run.sh -f hgdp_2chr/commandfiles/commandfile1.txt -n 8 -m 16 # choose -m to get the number of jobs we want. There
are 186 commands to run, so this gives us 12 jobs. Don't go above a couple of hundred.
fs hgdp_2chr.cp -indsperproc 100 -go # Set indsperproc again to control the amount of jobs. This is a short-term hack
because the default (1) produces too many files - 22 ⋆ 938 ⋆ 9 = 185,724 which breaks filesystem lookups
qsub_run.sh -f hgdp_2chr/commandfiles/commandfile2.txt -n 8 # We're already doing 100 inds per node requested. So there
are only 20 commands this time, and we don't need -m to run them serially
## Finestructure section, no more painting.  So things scale as N^3 with no dependence on L.
fs hgdp_2chr.cp -s3iters 1000000 -go
qsub_run.sh -f hgdp_2chr/commandfiles/commandfile3.txt -n 2
## But we received a the following messages:
#WARNING: Failed Gelman-Rubin MCMC diagnostic check with maximum potential scale reduction factor 1.39013 (threshold
1.3)
#WARNING: Stage 3 convergence criterion failed! Use "-ignoreGR" to continue regardless. (Set the parameter
"-threshGR:-1" to ignore the GR statistic in future.)  Re-running MCMC for longer: this is attempt 1 of 5.
## So we rerun, with automatically doubled duration:
fs hgdp_2chr.cp -go
qsub_run.sh -f hgdp_2chr/commandfiles/commandfile3.txt -n 2
## ... And again...
fs hgdp_2chr.cp -go
qsub_run.sh -f hgdp_2chr/commandfiles/commandfile3.txt -n 2
## ... And again...
fs hgdp_2chr.cp -go
qsub_run.sh -f hgdp_2chr/commandfiles/commandfile3.txt -n 2
### Finally, convergence is achieved. We can run the tree. This is quicker
fs hgdp_2chr.cp -go
qsub_run.sh -f hgdp_2chr/commandfiles/commandfile4.txt -n 2
## And we are done:
fs hgdp_2chr.cp -go

[next] [prev] [prev-tail] [front] [up]

PaintMyChromosomes.comfineSTRUCTURE v2 & GLOBETROTTER

fineSTRUCTURE

Data

Methodology

Other

Authors

13 Examples

13.1 Simple and quick example

13.2 Involved example with few SNPs and individuals

13.3 HGDP example, including downloading and processing.

PaintMyChromosomes.com
fineSTRUCTURE v2 & GLOBETROTTER