Genome Analysis Wiki - User contributions [en]

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T15:13:53Z

Chaolong Wang: /* Population information */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European or East Asian or Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

18 East Asian populations in the HGDP dataset, including 224 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

For European ancestry:
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

For East Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.easia.2.geno -c $HGDP/HGDP.633K.easia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-EAsia &

For Central/South Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.csasia.2.geno -c $HGDP/HGDP.633K.csasia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-CSAsia &

View the results for European ancestry:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot for European ancestry:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.easia.2.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.csasia.2.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:26:19Z

Chaolong Wang: /* Estimate ancestry */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European or East Asian or Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

For European ancestry:
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

For East Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.easia.2.geno -c $HGDP/HGDP.633K.easia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-EAsia &

For Central/South Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.csasia.2.geno -c $HGDP/HGDP.633K.csasia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-CSAsia &

View the results for European ancestry:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot for European ancestry:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.easia.2.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.csasia.2.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:17:05Z

Chaolong Wang: /* Visualizing Ancestry */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European or East Asian or Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

For European ancestry:
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

For East Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.easia.2.geno -c $HGDP/HGDP.633K.easia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-EAsia &

For Central/South Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.csasia.2.geno -c $HGDP/HGDP.633K.csasia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-CSAsia &

View the results for European ancestry:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot for European ancestry:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.easia.2.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.csasia.2.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:16:23Z

Chaolong Wang: /* Visualizing Ancestry */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European or East Asian or Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

For European ancestry:
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

For East Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.easia.2.geno -c $HGDP/HGDP.633K.easia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-EAsia &

For Central/South Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.csasia.2.geno -c $HGDP/HGDP.633K.csasia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-CSAsia &

View the results for European ancestry:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot for European ancestry:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.easia.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
Rscript plotHGDP.r $HGDP/HGDP.633K.csasia.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:15:40Z

Chaolong Wang: /* Estimate ancestry */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European or East Asian or Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

For European ancestry:
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

For East Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.easia.2.geno -c $HGDP/HGDP.633K.easia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-EAsia &

For Central/South Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.csasia.2.geno -c $HGDP/HGDP.633K.csasia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-CSAsia &

View the results for European ancestry:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or for East Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or for Central/South Asian ancestry
less -S $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or
Rscript plotHGDP.r $HGDP/HGDP.633K.easia.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or
Rscript plotHGDP.r $HGDP/HGDP.633K.csasia.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:14:53Z

Chaolong Wang: /* Visualizing Ancestry */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European or East Asian or Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

For European ancestry:
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

For East Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.easia.2.geno -c $HGDP/HGDP.633K.easia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-EAsia &

For Central/South Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.csasia.2.geno -c $HGDP/HGDP.633K.csasia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-CSAsia &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or
less -S $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or
less -S $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or
Rscript plotHGDP.r $HGDP/HGDP.633K.easia.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or
Rscript plotHGDP.r $HGDP/HGDP.633K.csasia.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:12:25Z

Chaolong Wang: /* Interested in looking just at European/East Asian/Central South Asian populations? */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European or East Asian or Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

For European ancestry:
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

For East Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.easia.2.geno -c $HGDP/HGDP.633K.easia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-EAsia &

For Central/South Asian ancestry:
$LASER/laser -g $HGDP/HGDP.633K.csasia.2.geno -c $HGDP/HGDP.633K.csasia.2.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-CSAsia &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord
or
less -S $OUT/ancestry/${SAMPLE}.HGDP-EAsia.SeqPC.coord
or
less -S $OUT/ancestry/${SAMPLE}.HGDP-CSAsia.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:08:59Z

Chaolong Wang: /* Interested in looking just at European (or East Asian or Central South Asian) populations? */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European/East Asian/Central South Asian populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:05:15Z

Chaolong Wang: /* Population information */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European (or East Asian or Central South Asian) populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |300px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-09T06:04:45Z

Chaolong Wang: /* HGDP reference panel */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== HGDP reference panel ==
[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |400px|HGDP populations]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r $LASER/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:02:50Z

Chaolong Wang: /* Population information */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European (or East Asian or Central South Asian) populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |600px|HGDP populations]]

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:02:17Z

Chaolong Wang: /* Population information */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European (or East Asian or Central South Asian) populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.

17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).

8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T06:02:04Z

Chaolong Wang: /* Interested in looking just at European populations? */

{{SeqShopLogin}}

== Setup ==
''If you were sequenced, set these values:''
export SAMPLE=Sample*
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

''If you were not sequenced, set these values:''
export SAMPLE=NA12878
source /net/seqshop-server/home/mktrost/seqshop/setupSS.txt
source /net/seqshop-server/home/chaolong/LASER-Tutorial/setup.txt

After setting this, also do
mkdir -p $OUT/ancestry

Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam

== Run ==

=== Step 1: bam --> pileup ===

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.HGDP.pileup

This step takes ~1 hour for a genome sequenced at 17X.

=== Step 2: pileup --> seq ===

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.HGDP \
$OUT/ancestry/$SAMPLE.HGDP.pileup

This step takes just a few seconds.

=== Estimate ancestry ===
This step will take about 5-6 minutes.

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/${SAMPLE}.HGDP.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

==Interested in looking just at European (or East Asian or Central South Asian) populations?==
=== Population information ===
8 European populations in the HGDP dataset, including 156 individuals in total.
17 East Asian populations in the HGDP dataset, including 215 individuals after excluding 4 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration).
8 Central/South Asian populations in the HGDP dataset, including 195 individuals after excluding 5 outliers (defined as more than 5 standard deviations from mean in any of top 10 PCs, 1 removal iteration). Note that 7 out of these 8 population were collected in Pakistan and 1 was collected in China (Uygur).

=== Step 1: bam --> pileup ===
You can skip this, you already did it.

=== Step 2: pileup --> seq ===

You can skip this step too, because you already did it.

=== Estimate ancestry ===
This step will take a few seconds.

$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.HGDP.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.HGDP-Euro &

View the results:
less -S $OUT/ancestry/${SAMPLE}.HGDP-Euro.SeqPC.coord

== Visualizing Ancestry ==
Copy the R code to plot your ancestry
cp -r $LASER/plot/ $OUT/ancestry/.

Change to that new directory:
cd $OUT/ancestry/plot

Move your other plot so you don't over-write it
mv Results_on_HGDP.pdf Results_on_HGDP_All.pdf

Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord

Take a look:
evince Results_on_HGDP.pdf &

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-09T06:00:34Z

Chaolong Wang: /* HGDP reference panel */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== HGDP reference panel ==
[[File:HGDP Popualtions.png|thumb|center|alt=HGDP populations |600px|HGDP populations]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r $LASER/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-09T06:00:07Z

Chaolong Wang: /* HGDP reference panel */

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-09T05:59:55Z

Chaolong Wang: /* HGDP reference panel */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== HGDP reference panel ==
[[File:HGDP Popualtions.png.png|thumb|center|alt=HGDP populations |400px|HGDP populations]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r $LASER/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-09T05:59:34Z

Chaolong Wang: /* LASER workflow */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== HGDP reference panel ==
[[File:File:HGDP Popualtions.png.png|thumb|center|alt=HGDP populations |400px|HGDP populations]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r $LASER/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

File:HGDP Popualtions.png

2014-12-09T05:57:44Z

Chaolong Wang:

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:34:08Z

Chaolong Wang: /* Estimate ancestry */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:32:51Z

Chaolong Wang: /* Step 2: pileup --> seq */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:32:33Z

Chaolong Wang: /* Step 2: pileup --> seq */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:32:02Z

Chaolong Wang: /* Step 2: pileup --> seq */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:30:26Z

Chaolong Wang: /* Visualizing Ancestry */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:29:30Z

Chaolong Wang: /* Estimate ancestry */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:28:05Z

Chaolong Wang: /* Step 2: pileup --> seq */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-09T05:27:28Z

Chaolong Wang: /* Step 1: bam --> pileup */

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-08T05:51:20Z

Chaolong Wang: /* Step 1: Estimate ancestry for sequenced samples */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r $LASER/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-08T05:46:17Z

Chaolong Wang: /* Step 1: bam --> pileup */

SeqShop: Ancestry On Your Own Genome, December 2014

2014-12-08T05:45:00Z

Chaolong Wang: /* Setup */

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-08T04:59:46Z

Chaolong Wang: /* Visualizing results */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

mosrun -e $LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
mosrun -e $LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r $LASER/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-08T04:56:00Z

Chaolong Wang: /* Step 1: Estimate ancestry for sequenced samples */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

mosrun -e $LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
mosrun -e $LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r ./LASER-2.01/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-08T04:54:07Z

Chaolong Wang: /* Estimating ancestry coordinates */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# $LASER/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r ./LASER-2.01/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-08T04:49:22Z

Chaolong Wang: /* Step 2: pileup --> seq */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# ./LASER-2.01/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

./LASER-2.01/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
./LASER-2.01/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r ./LASER-2.01/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-08T04:48:47Z

Chaolong Wang: /* Step 2: pileup --> seq */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python $LASER/LASER-2.01/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# ./LASER-2.01/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

./LASER-2.01/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
./LASER-2.01/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r ./LASER-2.01/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

SeqShop: Estimates of Genetic Ancestry Practical, December 2014

2014-12-07T23:45:38Z

Chaolong Wang: /* Setup your run environment */

== Introduction ==
See the [[Media:LASER-tutorial.pdf|tutorial slides]] for an introduction of the LASER analysis workflow, input/output file formats, and usage of the LASER software.

The main purpose of this page is to provide step-by-step command lines for using LASER to estimate ancestry of 6 targeted sequenced samples (2 HapMap trios) in a principal component space generated using genome-wide SNP data from the Human Genome Diversity Project (HGDP). The HGDP reference panel contains genotype data across 632,958 autosomal loci for 938 individuals from 53 populations worldwide.

For more details about the options and usage of LASER, please read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual].

== LASER workflow ==
[[File:LASER-workflow.png|thumb|center|alt=LASER workflow|400px|LASER workflow]]

== Setup in person at the SeqShop Workshop ==
''This section is specifically for the SeqShop Workshop computers.''
<div class="mw-collapsible mw-collapsed" style="width:600px">
''If you are not running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

{{SeqShopLogin}}

=== Setup your run environment===
This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

This will setup some environment variables to point you to
* Tutorial input files
source /home/chaolong/LASER-Tutorial/setup.txt
* You won't see any output after running <code>source</code>
** It silently sets up your environment

<div class="mw-collapsible mw-collapsed" style="width:400px">
What is in the setup.txt file:
<div class="mw-collapsible-content">
export GC=/home/mktrost/seqshop/gotcloud
export REF=/home/mktrost/seqshop/gotcloud/gotcloud.ref
export HGDP=/home/chaolong/LASER-Tutorial/HGDP
export BAM=/home/chaolong/LASER-Tutorial/BAM
export LASER=/home/chaolong/LASER-Tutorial/LASER-2.01
</div>
</div>

Point to where you want the output to go replacing the path with where you would like your output to go
<pre>export OUT=~/seqshop_output/</pre>
</div>
</div>

== Setup when running on your own outside of the SeqShop Workshop ==
''This section is specifically for running on your own outside of the SeqShop Workshop.''
<div class="mw-collapsible" style="width:600px">
''If you are running during the SeqShop Workshop, please skip this section.''
<div class="mw-collapsible-content">

This tutorial uses samtools from GotCloud, as well as example data downloaded in the Sequence Mapping & Assembly tutorial, so if you have not already installed GotCloud and the tutorial data in a previous tutorial, please do so now: [[SeqShop:_Sequence_Mapping_and_Assembly_Practical#Setup_when_running_on_your_own_outside_of_the_SeqShop_Workshop|Tutorial Setup]]

{{SeqShopRemoteEnv}}

<ul>
<li> Additional variables for Ancestry:</li>
<ul>
<div class="mw-collapsible" style="width:500px">
<li>Using bash (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>export REF=$SS/ancestry/ref
export HGDP=$SS/ancestry/HGDP
export BAM=$SS/ancestry/bams</pre>
</div>
</div>
<div class="mw-collapsible mw-collapsed" style="width:500px">
<li>Using tcsh (replace the paths below with the appropriate paths):</li>
<div class="mw-collapsible-content">
:<pre>setenv REF $SS/ancestry/ref
setenv HGDP $SS/ancestry/HGDP
setenv BAM $SS/ancestry/bams</pre>
</div>
</div>
</ul>
</ul>
</div>
</div>

== Getting started ==
Create a working directory:

mkdir $OUT/ancestry
cd $OUT/ancestry

Download and decompress software package:

wget http://www.sph.umich.edu/csg/chaolong/LASER/LASER-2.01.tar.gz
tar xzvf LASER-2.01.tar.gz

== Preparing input files for LASER ==
=== Step 0: vcf --> geno ===

This step prepares the reference panel by converting a VCF genotype file to a GENO file. We will skip this step and use a ready-to-use HGDP reference panel. A typical command to run the vcf2geno tool is given in the file "./LASER-2.01/vcf2geno/cmd.sh":
# cd ./LASER-2.01/vcf2geno/
# ./vcf2geno --inVcf exampleVCF/example.vcf.gz --updateID test.updateId --out test

=== Step 1: bam --> pileup ===

This step uses samtools to generate pileup files from bam files.

<div class="mw-collapsible mw-collapsed" style="width:500px">
In person at workshop notes:
<div class="mw-collapsible-content">
*Please only try one sample so that we won't overload the sever with everyone running 6 jobs at the same time. Pileup files for these 6 samples have been prepared for later steps.
*It takes about 2 mins for each pileup job.
</div>
</div>

<div class="mw-collapsible" style="width:500px">
Outside of the workshop notes:
<div class="mw-collapsible-content">
*The BAMs provided as part of the download are chr22 only BAMs. They are used to demonstrate how to run this step.
*Pileup files for the whole genome BAMs are provided with the download and will be used in the next step.
* You only need to try one of these.
</div>
</div>

$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101035.recal.bam > 121101035.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101043.recal.bam > 121101043.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101050.recal.bam > 121101050.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101052.recal.bam > 121101052.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101415.recal.bam > 121101415.recal.pileup
# $GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/human.g1k.v37.fa -l $HGDP/HGDP_938.bed $BAM/121101861.recal.bam > 121101861.recal.pileup

We use -q 30 and -Q 20 to exclude reads that have mapping quality score lower than 30 or base quality score lower than 20.

=== Step 2: pileup --> seq ===

In this step, we will generate a file called "hapmap_trios.seq", containing the information of 6 samples. It takes about 30 seconds to run.
We will use the pre-generated pileup files in the $BAM folder.
* These pre-generated pileup files are for the whole genome of all 6 samples

python ./LASER-2.01/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP_938.site \
-b $BAM/AMD_roi_1-based.bed \
-i $BAM/AMD_hapmap_trios_id.txt \
-o hapmap_trios \
$BAM/121101035.recal.pileup \
$BAM/121101043.recal.pileup \
$BAM/121101050.recal.pileup \
$BAM/121101052.recal.pileup \
$BAM/121101415.recal.pileup \
$BAM/121101861.recal.pileup

In the above command, -b provides the targeted regions to exclude and -i specifies alternative IDs for the BAM files to be used in the .seq file (including popID and indivID).
-b and -i are optional.

== Estimating ancestry coordinates ==

=== Step 0: Generate the reference ancestry space ===

LASER can perform principal components analysis (PCA) on genotype data of the reference panel to generate a reference ancestry space.

# ./LASER-2.01/laser -g $HGDP/HGDP_938.geno -pca 1 -k 30 -o HGDP_938

The above command takes ~20 minutes to finish.
We will skip this step, and use a set of reference ancestry coordinates that have been generated in the file $HGDP/HGDP_938.RefPC.coord.
View the reference coordinates:

less -S $HGDP/HGDP_938.RefPC.coord

=== Step 1: Estimate ancestry for sequenced samples ===

Submit two jobs to place sequenced samples into the reference ancestry space:

./LASER-2.01/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 1 -y 3 -o hapmap_trios.1-3 &
./LASER-2.01/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s hapmap_trios.seq -K 20 -k 4 -x 4 -y 6 -o hapmap_trios.4-6 &

The first job will process samples 1 to 3 and the second job will processed samples 4 to 6.
Each sequenced sample will be projected from a 20-dimensional PCA space onto a 4-dimensional reference ancestry space.
The running time is ~10 minutes for processing 3 samples in each job.

=== Step 2: Combine results ===

Results from previous step will be output to two files "hapmap_trios.1-3.SeqPC.coord" and "hapmap_trios.4-6.SeqPC.coord".
Here we simply concatenate the two files while skipping the header line of the second file.

cp hapmap_trios.1-3.SeqPC.coord hapmap_trios.SeqPC.coord
more +2 hapmap_trios.4-6.SeqPC.coord >> hapmap_trios.SeqPC.coord

View the results:
less -S hapmap_trios.SeqPC.coord

The results should look like below (results will vary slightly):

popID indivID L1 Ci K t PC1 PC2 PC3 PC4
YRI NA19238 78386 0.170864 20 0.999688 467.989 -210.294 -14.1729 -14.4204
CEU NA12892 85486 0.185973 20 0.999723 10.796 199.095 -9.90387 -21.4534
CEU NA12891 87588 0.190442 20 0.99973 2.04224 196.07 -19.5705 -12.8022
CEU NA12878 83213 0.181748 20 0.999711 4.34591 199.861 -12.4825 -22.6281
YRI NA19239 87564 0.193424 20 0.999734 474.464 -215.96 -9.02921 -19.7372
YRI NA19240 95866 0.213874 20 0.999748 469.914 -214.94 -14.9923 -13.6559

== Visualizing results ==

Example R codes are available in ./LASER-2.01/plot/. Let's copy the folder to current working directory:

cp -r ./LASER-2.01/plot/ ./

Go to the plot folder and run the script to plot results:
cd plot
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord ../hapmap_trios.SeqPC.coord

A figure named "Results_on_HGDP.pdf" will be generated.
Visualize the figure:

evince Results_on_HGDP.pdf &

We expect to see the following figure, in which 3 CEU samples cluster with HGDP Europeans and 3 YRI samples cluster with HGDP Africans:
[[File:Results_on_HGDP.png|thumb|center|alt=LASER results|400px|LASER results]]

LASER

2014-07-30T15:22:43Z

Chaolong Wang: /* Introduction */

= Introduction =

LASER, which stands for Locating Ancestry using SEquencing Reads, is a C++ software package that can estimate individual ancestry directly from genome-wide shortgun sequencing reads without calling genotypes. The method relies on the availability of a set of reference individuals whose genome-wide SNP genotypes and ancestral information are known. We first construct a reference coordinate system by applying principal components analysis (PCA) to the genotype data of the reference individuals. Then, for each sequencing sample, use the genome-wide sequencing reads to place the sample into the reference PCA space. With an appropriate reference panel, the estimated coordinates of the sequencing samples identify their ancestral background and can be directly used to correct for population structure in association studies or to ensure adequate matching of cases and controls.

Note:
The goal of this wiki page is to help you get start using LASER.
This page was created for LASER 1.0. Some of the information might be outdated for LASER 2.0.
A more updated wiki page can be found at [http://genome.sph.umich.edu/wiki/SeqShop:_Estimates_of_Genetic_Ancestry_Practical 2014 UM Sequencing Workshop].
We also encourage you to read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual] for more details of the software.

= Download =

To get a copy of the software and manual, go to the [http://www.sph.umich.edu/csg/chaolong/LASER/ LASER Download] page.

= Workflow =

LASER generates the coordinates from both reference individuals and sequence samples. It requires essentially two input files:

[[File:LASER-Workflow.png|thumb|center|alt=LASER workflow|400px|LASER Workflow]]

*Seq file: a text file processed from BAM (alignment) files. (See [[#Process sequencing file (BAM)|Processing sequencing file]] for how to prepare seq file)
*Geno file: genotypes of reference individuals. (See [[#Geno file|Geno file]] to understand geno file format)

LASER typically outputs two coord files: (1) in reference individuals' coord file(Reference.coord), LASER outputs the reference coordinates in the PCA space; (2) in sequence samples' coord files(AllSamples.coord), LASER infers their ancestries by placing their ancestry coordinates onto reference samples' PCA space.

An example result of the coord file of sequence samples is shown below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

In the header line, popID means "population ID", indivID means "individual ID", L1 means number of loci that has been covered by at least one read, Ci means "average coverage", t means Procrustes similarity. PC1 and PC2 mean coordinates of the first and second principal components.

= Tutorial =

In this tutorial, we will show you how to prepare data and run LASER.

== Process sequencing file (BAM) ==

We illustrate how to obtain .seq file from BAM files in this section.
In this example, we use HGDP data set as a reference, which contains 938 individuals and 632,958 markers.
[[File:LASER-DataProcessing.png|thumb|center|alt=LASER workflow|400px|LASER Data Processing Procedure]]

1. Obtain pileup files from BAM files

The first step is to generate a BED file:

cat ../resource/HGDP/HGDP_938.site |awk '{if (NR > 1) {print $1, $2-1, $2;}}' > HGDP_938.bed

This BED file contains the positions of all the reference markers.

Then we use ''samtools'' to extract the sequence bases overlapping these 632,958 reference markers.
Assuming your BAM file name is ''NA12878.chrom22.recal.bam'' (our example BAM file), you can use this:

samtools mpileup -q 30 -Q 20 -f ../../LASER-resource/reference/hs37d5.fa -l HGDP_938.bed exampleBAM/NA12878.chrom22.recal.bam > NA12878.chrom22.pileup

to obtain a pileup file named ''NA12878.chrom22.pileup''. It is required to keep the ''.pileup' suffix.

2. Obtain a seq file from pileup files.

After obtaining pileup files from each BAM file, you can convert them into a single seq file before running LASER.
Use the same site file and all generated pileup files from step 1 to generate a seq file:

python pileup2seq.py -m ../resource/HGDP/HGDP_938.site -o test NA12878.chrom22.pileup

You should obtain test.seq file after this step.

== Estimate ancestries of sequence samples ==

The easiest way to perform LASER using its exemplar data is:

./laser -s pileup2seq/test.seq -g resource/HGDP/HGDP_938.geno -c resource/HGDP/HGDP_938.RefPC.coord -o test -k 2

Upon successful calculation, you will find a result file "test.SeqPC.coord".

<br>

== Interpret LASER outputs ==

Upon successfully launching LASER command line as above, the output messages should be similar to below:

===================================================================
==== LASER: Locating Ancestry from SEquencing Reads ====
==== Version 1.0 | (c) Chaolong Wang 2013 ====
====================================================================
Started at: Fri Nov 15 01:05:48 2013

938 individuals are detected in the GENO_FILE.
632958 loci are detected in the GENO_FILE.
1 individuals are detected in the SEQ_FILE.
632958 loci are detected in the SEQ_FILE.
938 individuals are detected in the COORD_FILE.
100 PCs are detected in the COORD_FILE.

Parameter values used in execution:
-------------------------------------------------
GENO_FILE (-g)resource/HGDP/HGDP_938.geno
SEQ_FILE (-s)pileup2seq/test.seq
COORD_FILE (-c)resource/HGDP/HGDP_938.RefPC.coord
OUT_PREFIX (-o)test
DIM (-k)2
MIN_LOCI (-l)100
SEQ_ERR (-e)0.01
FIRST_IND (-x)1
LAST_IND (-y)1
REPS (-r)1
OUTPUT_REPS (-R)0
CHECK_FORMAT (-fmt)10
CHECK_COVERAGE (-cov)0
PCA_MODE (-pca)0
-------------------------------------------------

Fri Nov 15 01:05:50 2013
Checking data format ...
GENO_FILE: OK.
SEQ_FILE: OK.
COORD_FILE: OK.

Fri Nov 15 01:06:01 2013
Reading reference genotypes ...

Fri Nov 15 01:09:15 2013
Reading reference PCA coordinates ...

Fri Nov 15 01:09:15 2013
Analyzing sequence samples ...
Results for the sequence samples are output to 'test.SeqPC.coord'.

Finished at: Fri Nov 15 01:09:21 2013
====================================================================

The ancestry of input samples are store in the file '''test.SeqPC.coord''', which content is shown below:

popID indivID L1 Ci t PC1 PC2
NA12878.chrom22 NA12878.chrom22 1601 0.00858193 0.977243 31.522 224.098

The ancestry coordinates for NA12878 samples are given in PC1 (31.522) and PC2 (224.098).

It is recommended to visualize this results with HGDP reference samples whose coordinates are given in file: resource/HGDP/HGDP_938.RefPC.coord

In our manuscript, an example figure is shown:

[[File:LASER paper Figure 2.png|thumb|center|alt=LASER example outputs as in Figure 2|400px|LASER Outputs]]

In this figure, 238 individuals were randomly selected from the total 938 HGDP samples as the testing set (colored symbols),
and the remaining 700 HGDP individuals were used as the reference panel (gray symbols).

= File format =

== Geno file ==

Geno file are from reference samples. LASER use genotype of these samples as a reference panel. You can obtain geno file from VCF files using [https://github.com/zhanxw/vcf2geno vcf2geno].

In our resource folder, we provide an example geno file for the HGDP data set (resource/HGDP/HGDP_938.geno):

Brahui HGDP00001 1 2 1 1 0 2 0 2 1 2 2 2 1 1 2 1 0
Brahui HGDP00003 0 0 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00005 0 2 2 0 0 1 0 2 1 2 2 2 2 1 2 2 1
Brahui HGDP00007 0 2 2 0 0 2 0 2 0 2 2 2 1 1 2 2 1
Brahui HGDP00009 0 1 0 1 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00011 1 1 2 1 1 2 1 1 1 2 2 2 1 1 2 2 0
Brahui HGDP00013 1 2 2 1 1 2 1 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00015 1 1 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00017 1 1 2 0 0 1 0 0 0 2 0 1 1 2 2 2 0
Brahui HGDP00019 0 2 2 0 0 1 0 1 0 2 1 2 2 1 2 2 0

The first and second columns represent the population id and individual id.
From the third column, each number represents a genotype.
To be consistent with the sequence data, genotypes should be given on the '''forward strand'''. Genotypes are coded by 0, 1, or 2, representing copies of the
reference allele at a locus in one individual.

In this geno file, we have 632,960 columns which contains 632,958 markers from column 3 to the last column.

== Seq file ==
Seq file is generated from pileup files. It contains sequencing information and organize it in a LASER readable format.
The first two columns represent population id and individual id.
Subsequent columns are total read depths and reference base counts.
For example, column 3 and 4 are 0, 0 in the following example. That means at first marker, the sequence read depth is 0 and thus none of the reads has reference base.
We enforce tab delimiters between markers and space delimiters between each read depths and reference base counts.
On line of seq file looks like below:

NA12878.chrom22 NA12878.chrom22 0 0 0 0 0 0 0 0 0

== Pileup file ==

Pileup file are generated using samtools. An example pileup file is shown below:

22 17094749 A 1 c D
22 17202602 T 1 . D
22 17411899 A 1 . C
22 17450515 G 2 ., 9<
22 17452966 T 1 c 5
22 17470779 C 1 , A
22 17492203 G 1 , B
22 17504945 C 3 ,.. BCA
22 17529814 T 3 .., CCC

The columns are chromosome, position (1-based), reference base, depth, bases and base qualities.

== BED file ==
BED file represents genomic regions and it follows [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 UCSC conventions]:

1 752565 752566
1 768447 768448
1 1005805 1005806
1 1018703 1018704
1 1021414 1021415

The columns are: chromosome, start position (0-based) and end position (1-based).

== Coord file ==
Coord files represent the ancestries of both reference samples and sequence samples.
An example coord file looks like below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

The columns are: popID means "population ID", indivID means "individual ID", L1 means number of loci has been covered, Ci means "average coverage", t means Procrustes similarity.
PC1, PC2 means coordinates of first and second principal components. You may notice L1, Ci, and t are omitted in the coord files of reference samples. The reason is that reference samples use genotypes and do not have coverage information.

== Site file ==
Site file is equivalent to BED file and it is used here to represent marker positions. An example site file looks like below:
CHR POS ID REF ALT
1 752566 rs3094315 G A
1 768448 rs12562034 G A
1 1005806 rs3934834 C T
1 1018704 rs9442372 A G
1 1021415 rs3737728 A G

The site file has header line, and it contains chromosome, position(1-based), id (usually marker name), ref (reference allele) and alt (alternative allele).

= Advanced options =

LASER has advanced options including (1) parallel computing; (2) increase ancestry inference accuracy using repeated runs; (3) generate PCA coordiates using genotypes.
See [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf LASER Manual] for detailed information.

= Contact =
Comments on this wiki page or questions related to preparing input files for LASER can be sent to [mailto:zhanxw@umich.edu Xiaowei Zhan].
Comments on the LASER software or the user's manual can be sent to [mailto:chaolong@umich.edu Chaolong Wang].
This project was directed by Gonçalo Abecasis and Sebastian Zöllner at the University of Michigan.

LASER

2014-07-30T15:22:28Z

Chaolong Wang: /* Introduction */

= Introduction =

Note:
The goal of this wiki page is to help you get start using LASER.
This page was created for LASER 1.0. Some of the information might be outdated for LASER 2.0.
A more updated wiki page can be found at [http://genome.sph.umich.edu/wiki/SeqShop:_Estimates_of_Genetic_Ancestry_Practical 2014 UM Sequencing Workshop].
We also encourage you to read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual] for more details of the software.

LASER, which stands for Locating Ancestry using SEquencing Reads, is a C++ software package that can estimate individual ancestry directly from genome-wide shortgun sequencing reads without calling genotypes. The method relies on the availability of a set of reference individuals whose genome-wide SNP genotypes and ancestral information are known. We first construct a reference coordinate system by applying principal components analysis (PCA) to the genotype data of the reference individuals. Then, for each sequencing sample, use the genome-wide sequencing reads to place the sample into the reference PCA space. With an appropriate reference panel, the estimated coordinates of the sequencing samples identify their ancestral background and can be directly used to correct for population structure in association studies or to ensure adequate matching of cases and controls.

= Download =

To get a copy of the software and manual, go to the [http://www.sph.umich.edu/csg/chaolong/LASER/ LASER Download] page.

= Workflow =

LASER generates the coordinates from both reference individuals and sequence samples. It requires essentially two input files:

[[File:LASER-Workflow.png|thumb|center|alt=LASER workflow|400px|LASER Workflow]]

*Seq file: a text file processed from BAM (alignment) files. (See [[#Process sequencing file (BAM)|Processing sequencing file]] for how to prepare seq file)
*Geno file: genotypes of reference individuals. (See [[#Geno file|Geno file]] to understand geno file format)

LASER typically outputs two coord files: (1) in reference individuals' coord file(Reference.coord), LASER outputs the reference coordinates in the PCA space; (2) in sequence samples' coord files(AllSamples.coord), LASER infers their ancestries by placing their ancestry coordinates onto reference samples' PCA space.

An example result of the coord file of sequence samples is shown below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

In the header line, popID means "population ID", indivID means "individual ID", L1 means number of loci that has been covered by at least one read, Ci means "average coverage", t means Procrustes similarity. PC1 and PC2 mean coordinates of the first and second principal components.

= Tutorial =

In this tutorial, we will show you how to prepare data and run LASER.

== Process sequencing file (BAM) ==

We illustrate how to obtain .seq file from BAM files in this section.
In this example, we use HGDP data set as a reference, which contains 938 individuals and 632,958 markers.
[[File:LASER-DataProcessing.png|thumb|center|alt=LASER workflow|400px|LASER Data Processing Procedure]]

1. Obtain pileup files from BAM files

The first step is to generate a BED file:

cat ../resource/HGDP/HGDP_938.site |awk '{if (NR > 1) {print $1, $2-1, $2;}}' > HGDP_938.bed

This BED file contains the positions of all the reference markers.

Then we use ''samtools'' to extract the sequence bases overlapping these 632,958 reference markers.
Assuming your BAM file name is ''NA12878.chrom22.recal.bam'' (our example BAM file), you can use this:

samtools mpileup -q 30 -Q 20 -f ../../LASER-resource/reference/hs37d5.fa -l HGDP_938.bed exampleBAM/NA12878.chrom22.recal.bam > NA12878.chrom22.pileup

to obtain a pileup file named ''NA12878.chrom22.pileup''. It is required to keep the ''.pileup' suffix.

2. Obtain a seq file from pileup files.

After obtaining pileup files from each BAM file, you can convert them into a single seq file before running LASER.
Use the same site file and all generated pileup files from step 1 to generate a seq file:

python pileup2seq.py -m ../resource/HGDP/HGDP_938.site -o test NA12878.chrom22.pileup

You should obtain test.seq file after this step.

== Estimate ancestries of sequence samples ==

The easiest way to perform LASER using its exemplar data is:

./laser -s pileup2seq/test.seq -g resource/HGDP/HGDP_938.geno -c resource/HGDP/HGDP_938.RefPC.coord -o test -k 2

Upon successful calculation, you will find a result file "test.SeqPC.coord".

<br>

== Interpret LASER outputs ==

Upon successfully launching LASER command line as above, the output messages should be similar to below:

===================================================================
==== LASER: Locating Ancestry from SEquencing Reads ====
==== Version 1.0 | (c) Chaolong Wang 2013 ====
====================================================================
Started at: Fri Nov 15 01:05:48 2013

938 individuals are detected in the GENO_FILE.
632958 loci are detected in the GENO_FILE.
1 individuals are detected in the SEQ_FILE.
632958 loci are detected in the SEQ_FILE.
938 individuals are detected in the COORD_FILE.
100 PCs are detected in the COORD_FILE.

Parameter values used in execution:
-------------------------------------------------
GENO_FILE (-g)resource/HGDP/HGDP_938.geno
SEQ_FILE (-s)pileup2seq/test.seq
COORD_FILE (-c)resource/HGDP/HGDP_938.RefPC.coord
OUT_PREFIX (-o)test
DIM (-k)2
MIN_LOCI (-l)100
SEQ_ERR (-e)0.01
FIRST_IND (-x)1
LAST_IND (-y)1
REPS (-r)1
OUTPUT_REPS (-R)0
CHECK_FORMAT (-fmt)10
CHECK_COVERAGE (-cov)0
PCA_MODE (-pca)0
-------------------------------------------------

Fri Nov 15 01:05:50 2013
Checking data format ...
GENO_FILE: OK.
SEQ_FILE: OK.
COORD_FILE: OK.

Fri Nov 15 01:06:01 2013
Reading reference genotypes ...

Fri Nov 15 01:09:15 2013
Reading reference PCA coordinates ...

Fri Nov 15 01:09:15 2013
Analyzing sequence samples ...
Results for the sequence samples are output to 'test.SeqPC.coord'.

Finished at: Fri Nov 15 01:09:21 2013
====================================================================

The ancestry of input samples are store in the file '''test.SeqPC.coord''', which content is shown below:

popID indivID L1 Ci t PC1 PC2
NA12878.chrom22 NA12878.chrom22 1601 0.00858193 0.977243 31.522 224.098

The ancestry coordinates for NA12878 samples are given in PC1 (31.522) and PC2 (224.098).

It is recommended to visualize this results with HGDP reference samples whose coordinates are given in file: resource/HGDP/HGDP_938.RefPC.coord

In our manuscript, an example figure is shown:

[[File:LASER paper Figure 2.png|thumb|center|alt=LASER example outputs as in Figure 2|400px|LASER Outputs]]

In this figure, 238 individuals were randomly selected from the total 938 HGDP samples as the testing set (colored symbols),
and the remaining 700 HGDP individuals were used as the reference panel (gray symbols).

= File format =

== Geno file ==

Geno file are from reference samples. LASER use genotype of these samples as a reference panel. You can obtain geno file from VCF files using [https://github.com/zhanxw/vcf2geno vcf2geno].

In our resource folder, we provide an example geno file for the HGDP data set (resource/HGDP/HGDP_938.geno):

Brahui HGDP00001 1 2 1 1 0 2 0 2 1 2 2 2 1 1 2 1 0
Brahui HGDP00003 0 0 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00005 0 2 2 0 0 1 0 2 1 2 2 2 2 1 2 2 1
Brahui HGDP00007 0 2 2 0 0 2 0 2 0 2 2 2 1 1 2 2 1
Brahui HGDP00009 0 1 0 1 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00011 1 1 2 1 1 2 1 1 1 2 2 2 1 1 2 2 0
Brahui HGDP00013 1 2 2 1 1 2 1 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00015 1 1 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00017 1 1 2 0 0 1 0 0 0 2 0 1 1 2 2 2 0
Brahui HGDP00019 0 2 2 0 0 1 0 1 0 2 1 2 2 1 2 2 0

The first and second columns represent the population id and individual id.
From the third column, each number represents a genotype.
To be consistent with the sequence data, genotypes should be given on the '''forward strand'''. Genotypes are coded by 0, 1, or 2, representing copies of the
reference allele at a locus in one individual.

In this geno file, we have 632,960 columns which contains 632,958 markers from column 3 to the last column.

== Seq file ==
Seq file is generated from pileup files. It contains sequencing information and organize it in a LASER readable format.
The first two columns represent population id and individual id.
Subsequent columns are total read depths and reference base counts.
For example, column 3 and 4 are 0, 0 in the following example. That means at first marker, the sequence read depth is 0 and thus none of the reads has reference base.
We enforce tab delimiters between markers and space delimiters between each read depths and reference base counts.
On line of seq file looks like below:

NA12878.chrom22 NA12878.chrom22 0 0 0 0 0 0 0 0 0

== Pileup file ==

Pileup file are generated using samtools. An example pileup file is shown below:

22 17094749 A 1 c D
22 17202602 T 1 . D
22 17411899 A 1 . C
22 17450515 G 2 ., 9<
22 17452966 T 1 c 5
22 17470779 C 1 , A
22 17492203 G 1 , B
22 17504945 C 3 ,.. BCA
22 17529814 T 3 .., CCC

The columns are chromosome, position (1-based), reference base, depth, bases and base qualities.

== BED file ==
BED file represents genomic regions and it follows [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 UCSC conventions]:

1 752565 752566
1 768447 768448
1 1005805 1005806
1 1018703 1018704
1 1021414 1021415

The columns are: chromosome, start position (0-based) and end position (1-based).

== Coord file ==
Coord files represent the ancestries of both reference samples and sequence samples.
An example coord file looks like below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

The columns are: popID means "population ID", indivID means "individual ID", L1 means number of loci has been covered, Ci means "average coverage", t means Procrustes similarity.
PC1, PC2 means coordinates of first and second principal components. You may notice L1, Ci, and t are omitted in the coord files of reference samples. The reason is that reference samples use genotypes and do not have coverage information.

== Site file ==
Site file is equivalent to BED file and it is used here to represent marker positions. An example site file looks like below:
CHR POS ID REF ALT
1 752566 rs3094315 G A
1 768448 rs12562034 G A
1 1005806 rs3934834 C T
1 1018704 rs9442372 A G
1 1021415 rs3737728 A G

The site file has header line, and it contains chromosome, position(1-based), id (usually marker name), ref (reference allele) and alt (alternative allele).

= Advanced options =

LASER has advanced options including (1) parallel computing; (2) increase ancestry inference accuracy using repeated runs; (3) generate PCA coordiates using genotypes.
See [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf LASER Manual] for detailed information.

= Contact =
Comments on this wiki page or questions related to preparing input files for LASER can be sent to [mailto:zhanxw@umich.edu Xiaowei Zhan].
Comments on the LASER software or the user's manual can be sent to [mailto:chaolong@umich.edu Chaolong Wang].
This project was directed by Gonçalo Abecasis and Sebastian Zöllner at the University of Michigan.

LASER

2014-07-30T15:21:20Z

Chaolong Wang: /* Introduction */

LASER

2014-07-30T15:20:51Z

Chaolong Wang: /* Introduction */

= Introduction =

LASER, which stands for Locating Ancestry using SEquencing Reads, is a C++ software package that can estimate individual ancestry directly from genome-wide shortgun sequencing reads without calling genotypes. The method relies on the availability of a set of reference individuals whose genome-wide SNP genotypes and ancestral information are known. We first construct a reference coordinate system by applying principal components analysis (PCA) to the genotype data of the reference individuals. Then, for each sequencing sample, use the genome-wide sequencing reads to place the sample into the reference PCA space. With an appropriate reference panel, the estimated coordinates of the sequencing samples identify their ancestral background and can be directly used to correct for population structure in association studies or to ensure adequate matching of cases and controls.

Note:
The goal of this wiki page is to help you get start using LASER.
This page was created for LASER version 1.0. Some of the information might be outdated for LASER version 2.0.
A more updated wiki page can be found at [http://genome.sph.umich.edu/wiki/SeqShop:_Estimates_of_Genetic_Ancestry_Practical 2014 UM Sequencing Workshop].
We also encourage you to read the [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf manual] for more details of the software.

= Download =

To get a copy of the software and manual, go to the [http://www.sph.umich.edu/csg/chaolong/LASER/ LASER Download] page.

= Workflow =

LASER generates the coordinates from both reference individuals and sequence samples. It requires essentially two input files:

[[File:LASER-Workflow.png|thumb|center|alt=LASER workflow|400px|LASER Workflow]]

*Seq file: a text file processed from BAM (alignment) files. (See [[#Process sequencing file (BAM)|Processing sequencing file]] for how to prepare seq file)
*Geno file: genotypes of reference individuals. (See [[#Geno file|Geno file]] to understand geno file format)

LASER typically outputs two coord files: (1) in reference individuals' coord file(Reference.coord), LASER outputs the reference coordinates in the PCA space; (2) in sequence samples' coord files(AllSamples.coord), LASER infers their ancestries by placing their ancestry coordinates onto reference samples' PCA space.

An example result of the coord file of sequence samples is shown below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

In the header line, popID means "population ID", indivID means "individual ID", L1 means number of loci that has been covered by at least one read, Ci means "average coverage", t means Procrustes similarity. PC1 and PC2 mean coordinates of the first and second principal components.

= Tutorial =

In this tutorial, we will show you how to prepare data and run LASER.

== Process sequencing file (BAM) ==

We illustrate how to obtain .seq file from BAM files in this section.
In this example, we use HGDP data set as a reference, which contains 938 individuals and 632,958 markers.
[[File:LASER-DataProcessing.png|thumb|center|alt=LASER workflow|400px|LASER Data Processing Procedure]]

1. Obtain pileup files from BAM files

The first step is to generate a BED file:

cat ../resource/HGDP/HGDP_938.site |awk '{if (NR > 1) {print $1, $2-1, $2;}}' > HGDP_938.bed

This BED file contains the positions of all the reference markers.

Then we use ''samtools'' to extract the sequence bases overlapping these 632,958 reference markers.
Assuming your BAM file name is ''NA12878.chrom22.recal.bam'' (our example BAM file), you can use this:

samtools mpileup -q 30 -Q 20 -f ../../LASER-resource/reference/hs37d5.fa -l HGDP_938.bed exampleBAM/NA12878.chrom22.recal.bam > NA12878.chrom22.pileup

to obtain a pileup file named ''NA12878.chrom22.pileup''. It is required to keep the ''.pileup' suffix.

2. Obtain a seq file from pileup files.

After obtaining pileup files from each BAM file, you can convert them into a single seq file before running LASER.
Use the same site file and all generated pileup files from step 1 to generate a seq file:

python pileup2seq.py -m ../resource/HGDP/HGDP_938.site -o test NA12878.chrom22.pileup

You should obtain test.seq file after this step.

== Estimate ancestries of sequence samples ==

The easiest way to perform LASER using its exemplar data is:

./laser -s pileup2seq/test.seq -g resource/HGDP/HGDP_938.geno -c resource/HGDP/HGDP_938.RefPC.coord -o test -k 2

Upon successful calculation, you will find a result file "test.SeqPC.coord".

<br>

== Interpret LASER outputs ==

Upon successfully launching LASER command line as above, the output messages should be similar to below:

===================================================================
==== LASER: Locating Ancestry from SEquencing Reads ====
==== Version 1.0 | (c) Chaolong Wang 2013 ====
====================================================================
Started at: Fri Nov 15 01:05:48 2013

938 individuals are detected in the GENO_FILE.
632958 loci are detected in the GENO_FILE.
1 individuals are detected in the SEQ_FILE.
632958 loci are detected in the SEQ_FILE.
938 individuals are detected in the COORD_FILE.
100 PCs are detected in the COORD_FILE.

Parameter values used in execution:
-------------------------------------------------
GENO_FILE (-g)resource/HGDP/HGDP_938.geno
SEQ_FILE (-s)pileup2seq/test.seq
COORD_FILE (-c)resource/HGDP/HGDP_938.RefPC.coord
OUT_PREFIX (-o)test
DIM (-k)2
MIN_LOCI (-l)100
SEQ_ERR (-e)0.01
FIRST_IND (-x)1
LAST_IND (-y)1
REPS (-r)1
OUTPUT_REPS (-R)0
CHECK_FORMAT (-fmt)10
CHECK_COVERAGE (-cov)0
PCA_MODE (-pca)0
-------------------------------------------------

Fri Nov 15 01:05:50 2013
Checking data format ...
GENO_FILE: OK.
SEQ_FILE: OK.
COORD_FILE: OK.

Fri Nov 15 01:06:01 2013
Reading reference genotypes ...

Fri Nov 15 01:09:15 2013
Reading reference PCA coordinates ...

Fri Nov 15 01:09:15 2013
Analyzing sequence samples ...
Results for the sequence samples are output to 'test.SeqPC.coord'.

Finished at: Fri Nov 15 01:09:21 2013
====================================================================

The ancestry of input samples are store in the file '''test.SeqPC.coord''', which content is shown below:

popID indivID L1 Ci t PC1 PC2
NA12878.chrom22 NA12878.chrom22 1601 0.00858193 0.977243 31.522 224.098

The ancestry coordinates for NA12878 samples are given in PC1 (31.522) and PC2 (224.098).

It is recommended to visualize this results with HGDP reference samples whose coordinates are given in file: resource/HGDP/HGDP_938.RefPC.coord

In our manuscript, an example figure is shown:

[[File:LASER paper Figure 2.png|thumb|center|alt=LASER example outputs as in Figure 2|400px|LASER Outputs]]

In this figure, 238 individuals were randomly selected from the total 938 HGDP samples as the testing set (colored symbols),
and the remaining 700 HGDP individuals were used as the reference panel (gray symbols).

= File format =

== Geno file ==

Geno file are from reference samples. LASER use genotype of these samples as a reference panel. You can obtain geno file from VCF files using [https://github.com/zhanxw/vcf2geno vcf2geno].

In our resource folder, we provide an example geno file for the HGDP data set (resource/HGDP/HGDP_938.geno):

Brahui HGDP00001 1 2 1 1 0 2 0 2 1 2 2 2 1 1 2 1 0
Brahui HGDP00003 0 0 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00005 0 2 2 0 0 1 0 2 1 2 2 2 2 1 2 2 1
Brahui HGDP00007 0 2 2 0 0 2 0 2 0 2 2 2 1 1 2 2 1
Brahui HGDP00009 0 1 0 1 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00011 1 1 2 1 1 2 1 1 1 2 2 2 1 1 2 2 0
Brahui HGDP00013 1 2 2 1 1 2 1 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00015 1 1 2 0 0 2 0 2 0 2 2 2 2 0 2 2 0
Brahui HGDP00017 1 1 2 0 0 1 0 0 0 2 0 1 1 2 2 2 0
Brahui HGDP00019 0 2 2 0 0 1 0 1 0 2 1 2 2 1 2 2 0

The first and second columns represent the population id and individual id.
From the third column, each number represents a genotype.
To be consistent with the sequence data, genotypes should be given on the '''forward strand'''. Genotypes are coded by 0, 1, or 2, representing copies of the
reference allele at a locus in one individual.

In this geno file, we have 632,960 columns which contains 632,958 markers from column 3 to the last column.

== Seq file ==
Seq file is generated from pileup files. It contains sequencing information and organize it in a LASER readable format.
The first two columns represent population id and individual id.
Subsequent columns are total read depths and reference base counts.
For example, column 3 and 4 are 0, 0 in the following example. That means at first marker, the sequence read depth is 0 and thus none of the reads has reference base.
We enforce tab delimiters between markers and space delimiters between each read depths and reference base counts.
On line of seq file looks like below:

NA12878.chrom22 NA12878.chrom22 0 0 0 0 0 0 0 0 0

== Pileup file ==

Pileup file are generated using samtools. An example pileup file is shown below:

22 17094749 A 1 c D
22 17202602 T 1 . D
22 17411899 A 1 . C
22 17450515 G 2 ., 9<
22 17452966 T 1 c 5
22 17470779 C 1 , A
22 17492203 G 1 , B
22 17504945 C 3 ,.. BCA
22 17529814 T 3 .., CCC

The columns are chromosome, position (1-based), reference base, depth, bases and base qualities.

== BED file ==
BED file represents genomic regions and it follows [http://genome.ucsc.edu/FAQ/FAQformat.html#format1 UCSC conventions]:

1 752565 752566
1 768447 768448
1 1005805 1005806
1 1018703 1018704
1 1021414 1021415

The columns are: chromosome, start position (0-based) and end position (1-based).

== Coord file ==
Coord files represent the ancestries of both reference samples and sequence samples.
An example coord file looks like below:

popID indivID L1 Ci t PC1 PC2
YRI NA19238 1409 0.304122 0.98933 52.7634 -39.7924
CEU NA12892 1552 0.330037 0.989709 9.82674 25.2898
CEU NA12891 1609 0.362198 0.988082 0.439573 26.8872
CEU NA12878 1579 0.334825 0.988677 8.83775 28.1342
YRI NA19239 1558 0.34898 0.988302 53.9104 -39.1727
YRI NA19240 1735 0.404142 0.990264 59.8379 -45.2765

The columns are: popID means "population ID", indivID means "individual ID", L1 means number of loci has been covered, Ci means "average coverage", t means Procrustes similarity.
PC1, PC2 means coordinates of first and second principal components. You may notice L1, Ci, and t are omitted in the coord files of reference samples. The reason is that reference samples use genotypes and do not have coverage information.

== Site file ==
Site file is equivalent to BED file and it is used here to represent marker positions. An example site file looks like below:
CHR POS ID REF ALT
1 752566 rs3094315 G A
1 768448 rs12562034 G A
1 1005806 rs3934834 C T
1 1018704 rs9442372 A G
1 1021415 rs3737728 A G

The site file has header line, and it contains chromosome, position(1-based), id (usually marker name), ref (reference allele) and alt (alternative allele).

= Advanced options =

LASER has advanced options including (1) parallel computing; (2) increase ancestry inference accuracy using repeated runs; (3) generate PCA coordiates using genotypes.
See [http://www.sph.umich.edu/csg/chaolong/LASER/LASER_Manual.pdf LASER Manual] for detailed information.

= Contact =
Comments on this wiki page or questions related to preparing input files for LASER can be sent to [mailto:zhanxw@umich.edu Xiaowei Zhan].
Comments on the LASER software or the user's manual can be sent to [mailto:chaolong@umich.edu Chaolong Wang].
This project was directed by Gonçalo Abecasis and Sebastian Zöllner at the University of Michigan.

Template:SeqShopLogin

2014-06-26T15:52:35Z

Chaolong Wang: /* Login to the seqshop-server Linux Machine */

== Login to the seqshop-server Linux Machine ==
'' This section will appear redundantly in each session. If you already know how to log in to the server, please skip this section''

<ol>
<li> Login to the windows machine</li>
* The username/password for the Windows machine should be written on the right-hand monitor
<li> Start xming so you can open external windows on our Linux machine</li>
<ul>
<li> Start->Enter "Xming" in the search and select "Xming" from the program list</li>
<li> Nothing will happen, but Xming was started.</li>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Screenshot</li>
<div class="mw-collapsible-content">
[[File:xming.png]]
</div>
</div>
</ul>
<li> Open putty</li>
<ul>
<li> Start->Enter "putty" in the search and select "PuTTY" from the program list</li>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Screenshot</li>
<div class="mw-collapsible-content">
[[File:PuttyS.png]]
</div>
</div>
</ul>
<li> Configure PuTTY in the PuTTY Configuration window</li>
<ul>
<li> Host Name: <code>seqshop-server.sph.umich.edu</code></li>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Screenshot</li>
<div class="mw-collapsible-content">
[[File:seqshop.png]]
</div>
</div>
<li> Setup to allow you to open external windows:</li>
<ul>
<li> In the left pannel: Connection->SSH->X11</li>
<ul>
<li> Add a check mark in the box next to <code>Enable X11 forwarding</code></li>
<div class="mw-collapsible mw-collapsed" style="width:200px">
<li>View Screenshot</li>
<div class="mw-collapsible-content">
[[File:seqshopX11.png]]
</div>
</div>
</ul>
</ul>
<li> Click <code>Open</code></li>
<li> If it prompts about a key, click <code>OK</code></li>
</ul>
<li> Enter your provided username & password as provided</li>
</ol>

You should now be logged into a terminal on the seqshop-server and be able to access the test files.
* If you need another terminal, repeat from step 3.

File:LASER-tutorial.pdf

2014-06-19T18:30:02Z

Chaolong Wang: Chaolong Wang uploaded a new version of "File:LASER-tutorial.pdf"

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T03:21:01Z

Chaolong Wang: /* Step 1: bam --> pileup */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T03:20:38Z

Chaolong Wang: /* Step 1: bam --> pileup */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T03:00:32Z

Chaolong Wang: /* Visualizing results */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T02:42:04Z

Chaolong Wang: /* Visualizing results */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T02:37:46Z

Chaolong Wang: /* Visualizing results */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T02:36:22Z

Chaolong Wang: /* Step 1: bam --> pileup */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T02:34:17Z

Chaolong Wang: /* Step 1: bam --> pileup */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T02:33:00Z

Chaolong Wang: /* Getting started */

SeqShop: Estimates of Genetic Ancestry Practical, June 2014

2014-06-19T01:47:58Z

Chaolong Wang: /* Visualizing results */