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:14, 19 June 2014
, 14:14, 19 June 2014no edit summary
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/chaolong/LASER-Tutorial/setup.txt
source /home/mktrost/seqshop-server/setup.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:
export OUT=~/personal/YOUR_OUTPUT_DIR
Verify that this does not give an error:
ls $OUT/bams/${SAMPLE}.recal.bam
$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/hs37d5.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.recal.pileup
$GC/bin/samtools mpileup -q 30 -Q 20 -f $REF/hs37d5.fa -l $HGDP/HGDP_938.bed $OUT/bams/${SAMPLE}.recal.bam > $OUT/ancestry/${SAMPLE}.recal.pileup
This step takes 5-6 minutes.
=== Step 2: pileup --> seq ===
=== Step 2: pileup --> seq ===
python $LASER/pileup2seq/pileup2seq.py \
python $LASER/pileup2seq/pileup2seq.py \
-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_37224.recal.pileup
This step takes just a few seconds.
=== Estimate ancestry ===
=== Estimate ancestry ===
This step will take about 19 minutes.
$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_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 &
View the results:
less -S $OUT/ancestry/Sample_37224.laser.1.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_37224.laser.1.SeqPC.coord
Take a look:
evince Results_on_HGDP.pdf &
