把gwas信息转为bed格式

生信技能树

发布于 2022-03-03 14:07:17

8020

发布于 2022-03-03 14:07:17

文章被收录于专栏：生信技能树

有粉丝提问，他下载了 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}'|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的，进行如下所示的代码：

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的，进行如下所示代码：


 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列都是可以选择的。

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-02-06，如有侵权请联系 cloudcommunity@tencent.com 删除

typescript