有粉丝提问,他下载了 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.
十年后我环游世界各地的高校以及科研院所(当然包括中国大陆)的时候,如果有这样的情谊,我会优先见你。