标签归档:Base

Get the reference allele based on genomic position

This post will show how to get the reference base of chr1 from 49999 to 500001 (Version: hg19).
Please note: different tools has different coordinate (0 start or 1 start).

1, SAMtools

Index reference sequence in the FASTA format or extract subsequence from indexed reference sequence. If no region is specified, faidx will index the file and create .fai on the disk. If regions are specified, the subsequences will be retrieved and printed to stdout in the FASTA format.

$samtools faidx ucsc.hg19.fasta chr1:49999-50001
>chr1:49999-50001
ATA

2, twoBitToFa

$wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
$wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
$chmod 755 twoBitToFa faToTwoBit 
$faToTwoBit ucsc.hg19.fasta ucsc.hg19.2bit
$twoBitToFa ucsc.hg19.2bit -seq=chr1 -start=49998 -end=50001 temp.out && cat temp.out && rm temp.out
>chr1:49998-50001
ATA

3, UCSC DAS server

$wget -qO- http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr1:49999,50001 | grep -v '<'
ata

继续阅读