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

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