sh refseq chroms dict convert

Posted 2021-05-07

curl ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/H_sapiens/ARCHIVE/BUILD.37.3/Assembled_chromosomes/chr_NC_gi | \
awk 'BEGIN { OFS="\t"; FS="\t"; } NR > 1 && $2 ~ "NC_*" { print $2 " " $1 }' > acc_nums.txt

curl ftp://ftp.ncbi.nlm.nih.gov/genomes/H_sapiens/GRCh37.p13_interim_annotation/interim_GRCh37.p13_top_level_2017-01-13.gff3.gz > interim_GRCh37.p13_top_level_2017-01-13.gff3.gz
extract interim_GRCh37.p13_top_level_2017-01-13.gff3.gz

sed -f <(sed '/./!d;s/\([^ ]*\) *\(.*\)/\\|\1|s||\2|g/' acc_nums.txt ) interim_GRCh37.p13_top_level_2017-01-13.gff3 | grep 'exon' | sort -k1,1 -k2,2n  > hg19.gff3
bedtools sort -i hg19.gff3 | bgzip > hg19.gff3.gz
tabix -p gff hg19.gff3.gz

# Extract genes
sed -f <(sed '/./!d;s/\([^ ]*\) *\(.*\)/\\|\1|s||\2|g/' acc_nums.txt ) interim_GRCh37.p13_top_level_2017-01-13.gff3 | \
awk '$3 == "gene" { match($9, /gene=([^;])+;/, a); ; FS="\t"; OFS="\t"; print $1 "\t" $4 "\t" $5 "\t" substr($9, RSTART+5, RLENGTH-6);}' > gene.bed