3,045
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to searchSeqShop: Ancestry On Your Own Genome, June 2014 (view source)
Revision as of 14:32, 13 November 2014
, 14:32, 13 November 2014no edit summary
'''Note:''' the latest version of this practical is available at: [[SeqShop: Ancestry On Your Own Genome]]
* The ones here is the original one from the June workshop (updated to be run from elsewhere)
{{SeqShopLogin}}
{{SeqShopLogin}}
Set these values. If you used a different path for any of these, please update here. Also, be sure to specify your sample name instead of Sample_XXXXX
Set these values. If you used a different path for any of these, please update here. Also, be sure to specify your sample name instead of Sample_XXXXX
source /home/mktrost/seqshop-server/setup.txt
source /home/mktrost/seqshop/setup.2x.txt
export SAMPLE=Sample_XXXXX
export SAMPLE=Sample_XXXXX
OUT needs to point to where your alignment output went, so if your output is not ~/personal/output, please set OUT appropriately:
OUT needs to point to where your alignment output went, so if your output is not ~/personal/output, please set OUT appropriately:
export OUT=~/personal/YOUR_OUTPUT_DIR
export OUT=~/personal/YOUR_OUTPUT_DIR
After setting this, also do
mkdir -p $OUT/ancestry
Verify that this does not give an error:
Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam
ls $OUT/bams/${SAMPLE}.recal.bam
== Run ==
== Run ==
-m $HGDP/HGDP_938.site \
-m $HGDP/HGDP_938.site \
-o $OUT/ancestry/$SAMPLE.laser \
-o $OUT/ancestry/$SAMPLE.laser \
$OUT/ancestry/Sample_37224.recal.pileup
$OUT/ancestry/$SAMPLE.recal.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.laser.2 &
View the results:
less -S $OUT/ancestry/${SAMPLE}.laser.2.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}.laser.2.SeqPC.coord
Take a look:
evince Results_on_HGDP.pdf &
==Interested in looking just at European populations?==
=== Step 1: bam --> pileup ===
You can skip this, you already did it.
=== Step 2: pileup --> seq ===
python $LASER/pileup2seq/pileup2seq.py \
-m $HGDP/HGDP.633K.euro.site \
-o $OUT/ancestry/$SAMPLE.Euro.laser \
$OUT/ancestry/${SAMPLE}.recal.pileup
This step takes just a few seconds.
This step takes just a few seconds.
=== Estimate ancestry ===
=== Estimate ancestry ===
This step will take about 19 minutes.
This step will take a few seconds.
$LASER/laser -g $HGDP/HGDP_938.geno -c $HGDP/HGDP_938.RefPC.coord -s $OUT/ancestry/$SAMPLE.laser.seq -K 20 -k 4 -x 1 -y 3 -o $OUT/ancestry/$SAMPLE.laser.1 &
$LASER/laser -g $HGDP/HGDP.633K.euro.geno -c $HGDP/HGDP.633K.euro.RefPC.coord -s $OUT/ancestry/$SAMPLE.Euro.laser.seq -K 20 -k 4 -M 0.8 -o $OUT/ancestry/$SAMPLE.Euro.laser.1 &
View the results:
View the results:
less -S $OUT/ancestry/Sample_37224.laser.1.SeqPC.coord
less -S $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord
== Visualizing Ancestry ==
== Visualizing Ancestry ==
Change to that new directory:
Change to that new directory:
cd $OUT/ancestry/plot
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:
Generate the plot:
Rscript plotHGDP.r $HGDP/HGDP_938.RefPC.coord $OUT/ancestry/Sample_37224.laser.1.SeqPC.coord
Rscript plotHGDP.r $HGDP/HGDP.633K.euro.RefPC.coord $OUT/ancestry/${SAMPLE}.Euro.laser.1.SeqPC.coord
Take a look:
Take a look:
evince Results_on_HGDP.pdf &
evince Results_on_HGDP.pdf &
