genetics, computing, life

Download ART http://www.niehs.nih.gov/research/resources/software/biostatistics/art/ 

linux64x: http://www.niehs.nih.gov/research/resources/assets/docs/artbinchocolatecherrycake031915linux64tgz.tgz 

 Download sequence of interest Example Human Potassium Channels: http://www.ncbi.nlm.nih.gov/gene/?term=potassium+channel+homo

let’s use KCNQ1

http://www.ncbi.nlm.nih.gov/gene/3784

Go to sequences http://www.ncbi.nlm.nih.gov/gene/3784#

http://www.ncbi.nlm.nih.gov/nuccore/NM_000218.2

Choose one of mRNAs

Get it in FASTA format

>gi|32479526|ref|NM_000218.2| Homo sapiens potassium voltage-gated channel, KQT-like subfamily, member 1 (KCNQ1), transcript variant 1, mRNA

GCGGCGGGGCTGGCAGCAGTGGCTGCCCGCACTGCGCCCGGGCGCTCGCCTTCGCTGCAGCTCCCGGTGC CGCCGCTCGGGCCGGCCCCCCGGCAGGCCCTCCTCGTTATGGCCGCGGCCTCCTCCCCGCCCAGGGCCGA GAGGAAGCGCTGGGGTTGGGGCCGCCTGCCAGGCGCCCGGCGGGGCAGCGCGGGCCTGGCCAAGAAGTGC CCCTTCTCGCTGGAGCTGGCGGAGGGCGGCCCGGCGGGCGGCGCGCTCTACGCGCCCATCGCGCCCGGCG

AND remove some seq in beginning

>gi|32479526|ref|NM_000218.2| Homo sapiens potassium voltage-gated channel, KQT-like subfamily, member 1 (KCNQ1), transcript variant 1, mRNA GCGGCGGGGCTGGCAGCAGTGGCTGCCCGCACTGCGCCCGGGCGCTCGCCTTCGCTGCAGCTCCCGGTGC CGCCGCTCGGG[CCGGCC]CCCCGGCAGGCCCTCCTCGTTATGGCCGCGGCCTCCTCCCCGCCCAGGGCCGA GAGGAAGCGCTGGGGTTGGGGCCGCCTGCCAGGCGCCCGGCGGGGCAGCGCGGGCCTGGCCAAGAAGTGC CCCTTCTCGCTGGAGCTGGCGGAGGGCGGCCCGGCGGGCGGCGCGCTCTACGCGCCCATCGCGCCCGGCG

Then let’s synthese reads

Normal

art_illumina -i KCNQ1_ts1.fasta -l 75 -ss HS25 -f 50 -o KCNQ1_ts

Mutated

art_illumina -i KCNQ1_ts1_del_CCGGCC.fasta -l 75 -ss HS25 -f 50 -o KCNQ1_ts1_del_CCGGCC

Now align it to hg19 and see in browser

bwa mem ~/ref/bwa/hg19.fa KCNQ1_ts1_del_CCGGCC.fq | awk -F'\t' '{if($4 != 0){print;}}' | samtools view -Shu - | samtools sort - KCNQ1_ts1_del_CCGGCC.fq.sorted

samtools index KCNQ1_ts1_del_CCGGCC.fq.sorted.bam

bwa mem ~/ref/bwa/hg19.fa KCNQ1_ts1.fq | awk -F'\t' '{if($4 != 0){print;}}' | samtools view -Shu - | samtools sort - KCNQ1_ts1.fq.sorted

samtools index KCNQ1_ts1.fq.sorted.bam

See results in IGV genome browser

I want to have file with all snps from snp138 database after exome variant calling - so my goal is to create a compatible bed file. It was easy to find description here: http://mendel.iontorrent.com/ion-docs/BED-File-Formats-and-Examples_64062608.html#BEDFileFormatsandExamples-PartialexampleofaHotSpotsBEDfile lala fa 1. Get all snp's plus reference alleles:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e ' select S.chrom, S.chromStart, S.chromEnd, S.name, concat("REF=",S.refNCBI,";OBS=;"), "." from snp138 as S limit 10;'

2. Select snps in exome target region

intersectBed -a snp138_all.bed -b AmpliSeqExome.20130612.Designed.bed > snp138_all_in_Ampliseq.bed