6p21.31这种表示方式在人类遗传学中用得比较多。
在参考基因组出来之前，人们通过做实验染色的方法将染色体染成很多段（就是那种很经典的黑白灰分段染色体图，分辨率现在可以达到大几百上千条带），用特定的编号来表示染色体位置，在此基础上再去研究染色体结构正常或异常对表型的影响。
6p21.31
- 6表示第6号染色体
- p表示染色体短臂，与之对应的q表示长臂
- 21.31表示2区1带3亚带1次亚带
怎么查询6p21.31的物理位置？
点击go之后就能出现物理位置和区段长度了
但其实这些区段（Chromosome Bands）信息是可以下载的，链接是：ftp://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/，在这个文件夹下面搜索band关键字即可找到。
$ less cytoBand.txt | grep "chr6.*p21\.31" | cut -f1-3 > 6p21.31.bed
(base) hsy 00:51:18 ~/bioinformatics/ref/Human/ensembl
$ cat 6p21.31.bed
chr6 33500000 36600000
注意：这个文件是起始位置为0，左闭右开的bed格式，和上面那个查询结果有点区别。
6p21.31上面有多少基因？
需要用到bedtools
$ less -S -x 20 Homo_sapiens.GRCh38.98.chr.gtf.gz | awk -F "\t" '{if($3=="gene") print "chr"$1"\t"$4-1"\t"$5}' > tmp1
$ less -S -x 20 Homo_sapiens.GRCh38.98.chr.gtf.gz | awk -F "\t" '{if($3=="gene") print $9}' > tmp2
$ cat test.pl
#! /usr/bin/perl
use warnings;
use strict;
my $file_name = $ARGV[0];
open my $in_fh, "<", "$file_name";
while (<$in_fh>) {
chomp $_;
my @one_line = (split("\t", $_));
my @old_name = ( $one_line[0] =~ /gene_name \".*?\"/g );
print "@old_name\n";
}
close $in_fh;
$ perl test.pl tmp2 | cut -d "\"" -f2 > tmp3
$ paste tmp1 tmp3 > gene.bed
$ rm -f tmp*
$ head -n 4 gene.bed
chr1 11868 14409 DDX11L1
chr1 14403 29570 WASH7P
chr1 17368 17436 MIR6859-1
chr1 29553 31109 MIR1302-2HG
$ bedtools intersect -a 6p21.31.bed -b gene.bed -wb | sort -k1,1 -k2,2n | cut -f7 | head
GGNBP1
RN7SL26P
Metazoa_SRP
BAK1
LINC00336
ITPR3
UQCC2
MIR3934
IP6K3
LEMD2
这个区段有92个基因。
最后需要注意的是，Chromosome Bands文件也是需要区分参考基因组版本的。
网友评论