下载了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文件，非常的复杂，列比较多，如下所示：

$ 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脚本处理即可，如下所示：

$ cat gwas_catalog_v1.0.2-associations_e105_r2021-12-21.tsv |awk -F"\t" '{print $12,$13,$13+1,$21}'|headCHR_ID CHR_POS 1 STRONGEST SNP-RISK ALLELE9 82904977 82904978 rs7045089-A20 49084487 49084488 rs2075677-A3 43812828 43812829 rs6772840-A5 152808169 152808170 rs4958581-T3 158470022 158470023 rs6787172-T20 32643466 32643467 rs6141769-C10 3668695 3668696 rs10795061-T1 20847936 20847937 rs17449243-T5 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的，进行如下所示的代码：

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 > tmpawk '{print "chr"$0}' tmp |uniq > gwas_with_rs.bed#一般来说，bed文件是不足够的，需要排序，并且加上chr的标签，上面的示例代码需要一个软件。

然后，对那些不含dbsnp数据库的ID的，进行如下所示代码：

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.txtchr1:158585415-Tchr7:135628370-Cchr9:135864436-Gchr14:66113725-?chr2:48696432-?chr5:75996909-Achr20:42658274-Cchr22:37462936-Gchr10:94450233-Tchr9:6282511-A

前面的gwas_with_rs.bed我测试了，没有问题，是一个正常的bed文件，后面的gwas_without_rs.txt里面的信息也足矣做出bed格式，就交给大家作为练习题吧。

记住：bed格式最重要的是前面的3列，也就是染色体编号，起始终止坐标即可，剩余的3列或者6列都是可以选择的。

转载自《生信技能树》，如有侵权，请联系删除。