bedtools查找基因组位置的信息

从gencode下载注释数据gff3

我们需要的是第一列，第三列，第四列，第五列
去掉注释部分起名tt
然后使用awk 获取这些列，分隔符为空格

awk -F " " '{print $1,$4,$5,$3}' tt > GR37.annotation

由于bed文件是制表符分隔的，所以我们需要把整理好的annotation文件变为制表符分隔

awk 'BEGIN{ FS=" ";OFS="\t" }{ print $1,$2,$3,$4 }' GR37.annotation > GR37.annotation1

这里的GR37.annotation1就是我们需要的文件了
由于基因位置在这个文件中包含10种，有很多是重复的，所以把CDS，gene，transcript去掉

然后把gross deletion的数据整理成bed文件

gross是逗号分隔的文件，把染色体和起始位置提取出来

awk -F "," '{print $12,$14,$15}'  GROSS_DELETIONS_hg19_ALL_TAGS_HGMD_PRO_19.3.csv >  temp

离我们的bed文件还需要把表头删掉，在染色体数字前加chr，然后改成制表符分隔
这个就留给你当作业吧
在每行前面加chr的代码是

sed 's/^/chr&/g' file > file1

最后成果

最后用bedtools

bedtools intersect -a ~/work/HGMD/data/gross.bed  -b  GR37.bed -wa -wb > t1

结果：前三列是gross文件里的信息，后四列是reference文件里的信息
比对回gross文件，发现有的位置没有，就说明其在intron区域，补齐就行了

这种方法可以获得断裂片段的位置信息，若只是想知道断点的信息，更改bed文件的染色体位置就行了

bedtools intersect -a ~/work/HGMD/data/gross.bed  -b  GR37.bed -wa -wb |bedtools groupby -i - -g 1-3 -c 7 -o collapse> t2

可以这样写整理一下