User:Timothee Flutre/Notebook/Postdoc/2012/01/30

From OpenWetWare

< User:Timothee Flutre | Notebook | Postdoc | 2012 | 01
Revision as of 11:28, 30 January 2012 by Timothee Flutre (Talk | contribs)
(diff) ←Older revision | Current revision (diff) | Newer revision→ (diff)
Jump to: navigation, search
Project name Main project page
Previous entry      Next entry

Use BWA

  • index the whole genome:
echo "bwa index -p wholegenome_bwaidx -a bwtsw wholegenome.fa" | qsub -l h_vmem=5g -N job_wholegenome_bwaidx -j y -cwd -V
  • map short sequences on it:
fasta2fastq.py shortseq.fa shortseq.fq 1
echo "bwa aln wholegenome_bwaidx shortseq.fq > shortseq_bwa_aln.sai" | qsub -l h_vmem=5g -N job_bwa_aln -j y -cwd -V
echo "bwa samse wholegenome_bwaidx shortseq_bwa_aln.sai shortseq.fq > shortseq_bwa_samse.sam" | qsub -l h_vmem=5g -N job_bwa_samse -j y -cwd -V
grep -v "@" shortseq_bwa_samse.sam | awk -F"\t" 'BEGIN{print "flag\toccurrences"} {a[$2]++} END{for(i in a)print i"\t"a[i]}'
parse_sam.py -i shortseq_bwa_samse.sam -e 2 -o shortseq_bwa_uniq_perfect_e2
  • Example for hg19:
wget --timestamping 'ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz'
tar xzvf chromFa.tar.gz
rm -f hg19_autosomes.fa; for i in {1..22}; do cat "chr"$i".fa" >> hg19_autosomes.fa; done
rm -f hg19_autosomes_sex.fa; for i in "hg19_autosomes.fa" "chrX.fa" "chrY.fa"; do cat $i >> hg19_autosomes_sex.fa; done


Personal tools