有粉丝提问，他下载了 gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv 文件，希望我可以帮忙看看他自己的一些表观调控区域里面是否有这些gwas位点信息，也就是说两个坐标文件的交集。
这个时候需要把两个文件都弄成为bed格式， 然后使用 bedtools "intersect" 命令即可。
看了看他下载的 gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv 文件，非常的复杂， 列比较多，如下所示：
代码语言：javascript复制$ cat gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv | head -1 |tr 't' 'n' |cat -n
1 DATE ADDED TO CATALOG
2 PUBMEDID
3 FIRST AUTHOR
4 DATE
5 JOURNAL
6 LINK
7 STUDY
8 DISEASE/TRAIT
9 INITIAL SAMPLE SIZE
10 REPLICATION SAMPLE SIZE
11 REGION
12 CHR_ID
13 CHR_POS
14 REPORTED GENE(S)
15 MAPPED_GENE
16 UPSTREAM_GENE_ID
17 DOWNSTREAM_GENE_ID
18 SNP_GENE_IDS
19 UPSTREAM_GENE_DISTANCE
20 DOWNSTREAM_GENE_DISTANCE
21 STRONGEST SNP-RISK ALLELE
22 SNPS
23 MERGED
24 SNP_ID_CURRENT
25 CONTEXT
26 INTERGENIC
27 RISK ALLELE FREQUENCY
28 P-VALUE
29 PVALUE_MLOG
30 P-VALUE (TEXT)
31 OR or BETA
32 95% CI (TEXT)
33 PLATFORM [SNPS PASSING QC]
34 CNV
35 MAPPED_TRAIT
36 MAPPED_TRAIT_URI
37 STUDY ACCESSION
38 GENOTYPING TECHNOLOGY
bed格式最重要的是前面的3列，也就是 染色体编号，起始终止坐标即可，剩余的3列或者6列都是可以选择的。所以我只需要简单的awk脚本处理即可， 如下所示：
代码语言：javascript复制$ cat gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv |awk -F"t" '{print $12,$13,$13 1,$21}'|head
CHR_ID CHR_POS 1 STRONGEST SNP-RISK ALLELE
9 82904977 82904978 rs7045089-A
20 49084487 49084488 rs2075677-A
3 43812828 43812829 rs6772840-A
5 152808169 152808170 rs4958581-T
3 158470022 158470023 rs6787172-T
20 32643466 32643467 rs6141769-C
10 3668695 3668696 rs10795061-T
1 20847936 20847937 rs17449243-T
5 104356621 104356622 rs10067176-A
# 命令如下所示：
cat gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv |awk -F"t" '{print $12,$13,$13 1,$21}' > gwas.bed
但是这个bed文件里面很多非法字符，非常的不干净。我上面的命令输出的bed文件在后续分析，总是被bedtools等工具给出报错。
我们可以尝试根据里面的信息，是否含有dbsnp数据库的ID来进行分类讨论，分别输出不同bed文件，再合并。
首先，对gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv 文件里面含有 dbsnp数据库的ID的，进行如下所示的代码：
代码语言：javascript复制grep rs gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv |perl -F"t" -alne '{print if $F[11] }'|awk -F"t" '{print $12,$13,$13 1 }'| grep -v ';'| grep -v 'x'|uniq > gwas_with_rs.bed
# 这个代码又臭又长， 就是因为这个gwas_catalog文件实在是太恶心了，里面很多非法字符。
sort-bed gwas_with_rs.bed > tmp
awk '{print "chr"$0}' tmp |uniq > gwas_with_rs.bed
# 一般来说， bed文件是不足够的，需要排序，并且加上 chr的标签，上面的示例代码需要一个软件。
然后，对那些不含dbsnp数据库的ID的，进行如下所示代码 ：
代码语言：javascript复制
grep -v rs gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv |awk -F"t" '{print $21}'|grep chr > gwas_without_rs.txt
(m6a) jmzeng 15:57:34 ~/m6a
$ head gwas_without_rs.txt
chr1:158585415-T
chr7:135628370-C
chr9:135864436-G
chr14:66113725-?
chr2:48696432-?
chr5:75996909-A
chr20:42658274-C
chr22:37462936-G
chr10:94450233-T
chr9:6282511-A
前面的 gwas_with_rs.bed 我测试了，没有问题， 是一个正常的bed文件，后面的 gwas_without_rs.txt 里面的信息也足矣做出bed格式。
记住：bed格式最重要的是前面的3列，也就是 染色体编号，起始终止坐标即可，剩余的3列或者6列都是可以选择的。
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
代码语言：javascript复制We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。
