打开 https://www.broadinstitute.org/medical-and-population-genetics/hapmap-3, 点击 How To Download This Release 下面的 A. SNP Genotype Data 段落的中间3个链接。 文件名字里面有 "b36",现在一般都用 b37(比如 UK Biobank),甚至有的用 b38, 所以下载后解压后需要将那个 .map 文件先用 liftOver 转化为 b37 格式,然后用 PLINK 生成 bed/bim/fam 文件。 这个基因数据可供 LDSC 和 GSMR 等软件使用。
在千人基因组官网 下载 Phase 3 对应的 VCF 链接, 有一个文件罗列了每一个样本的人群(pop)和人种 (super_pop),以及性别,可以用PLINK --keep 选取特定人种的样本。 下载下来的数据,有将近一个亿的SNP,每个染色体都是单独的文件。后续跑 GWAS 或提取 PRS 的时候,也是每条染色体的数据分开来跑。 PLINK的网站上也有“1000 genomes phase3” 数据。PLINK 不允许 SNP 名字有重复,可以用下面的命令来处理。
awk '{if(array[$2]=="Y") {i++; $2=$2".DUP"i}; print $0; array[$2]="Y"}' chr1.bim.COPY > chr1.bim
阅读 Data access guide 文件,里面会提到如何下载用来下载UKB数据的小软件(比如ukbunpack和unkconv)。
申请得到批准后,从最上面的 Researcher log in 登录后获取。基因型数据我已下载到南科大的HPC上,表型数据见百度网盘上的 ukb50136.enc。
ukbiobank 官网,点击 data showcase --> Essential information --> Accessing UK Biobank data。具体的流程和代码,请见 scripts 文件夹下的 phe.sh, gen.sh 以及 phe.R 和 gen.R。
通过 ukbconv 提取很多列的时候,可以先写一个 MY.fields.txt 文件,列出想提取的变量和对应的 data-field,比如第一行是sex 31,第二行是 age 21022,等。然后用 ukbconv ukb50136.enc_ukb r -iMY.fields.id -oMY 对于表型数据的提取,有一个 ukbtools R软件包。 但不是太好用,并且很慢,可供参考。 用下面的代码将 icd.tab 文件整合为两列,便于读入R。
cat icd.tab | sed -e 's/\tNA//g' -e 's/\t/,/2g' | \
awk '{ if(NR==1) print "IID icd"; else if (NF==1) print $1 " NA"; else print $0"," }' > icd.2cols
提取需要研究的表型数据和相关的covariates,比如 age, sex, PCs。 一般来说,quantitative的表型数据要 adjust for covariates 和转化成正态分布,这个可以在R里面用下面的命令来实现。
trait_res = residuals(lm(trait ~ age+sex+PC1+PC2, na.action=na.exclude)
trait_inv = qnorm((rank(trait_res,na.last="keep")-0.5) / length(na.omit(trait_res)))
对于疾病的binary 表型,只需要把需要 adjust 的covarites 和表型数据放在同一个表型数据文件里面, 然后在 GWAS里面的plink命令指明哪个是表型,哪些是 covariates。 目前GWAS 由专人负责运行,一般来说就是通过下面这样的PLINK命令来跑
for chr in {1..22}; do
plink2 --memory 12000 --threads 16 --pfile chr$chr --extract ukb.chr$chr.good.snps --pheno cvd.EUR.pheno --no-psam-pheno --pheno-name XXX --1 --glm cols=+ax,+a1freq,+a1freqcc,+a1count,+a1countcc,+beta,+orbeta,+nobs hide-covar no-x-sex --covar pheno/ukb.cov --covar-name age,sex,PC1-PC10 --out chr$chr
done
上述命令顺利跑完后,确认生成的文件没有问题后,可以把所有的染色体的数据串到一起,形成一个单一的 XXX.gwas.gz 文件。 最终合并成的 XXX.gwas.gz 文件用 TAB 分割,CHR:POS 排好序,要不然 LocusZoom 那样的软件不能处理。也可以用 tabix -f -S 1 -s 1 -b 2 -e 2 XXX.gwas.gz 对数据进行索引,便于 LocalZoom 那样的软件去处理。
最经典的,起源于美国NIH 的 GWAS Catalog. 这个页面也罗列了一些大型GWAS数据联盟。 欧洲版本,不需要下载就能通过 TwoSampleMR 远程读入。他们提倡 使用 VCF 格式的GWAS文件。
UKB GWAS 完整的分析结果,网上发布
- 美国哈佛大学:http://www.nealelab.is/uk-biobank
- 英国爱丁堡大学:geneatlas: http://geneatlas.roslin.ed.ac.uk
- 哈佛大学的 CVD knowlege portal: https://hugeamp.org/
可使用密西根大学开发的Pheweb 流水线作业。日本版本pheweb.jp。**版本的是本课题组建立的 pheweb.cn。 Pheweb有一个强大的add_rsids.py 的功能,但是存在先天缺陷。根据该聊天记录,用户可以在安装pheweb 后找到 add_rsids.py 文件(find /home/ -name "add_rsid*"),修改一行代码(第140行)。。
修改前:rsids = [rsid['rsid'] for rsid in rsid_group if cpra['ref'] == rsid['ref'] and are_match(cpra['alt'], rsid['alt'])]
修改后:rsids = [rsid['rsid'] for rsid in rsid_group if (cpra['ref'] == rsid['ref'] and are_match(cpra['alt'], rsid['alt'])) or (cpra['ref'] == rsid['alt'] and are_match(cpra['alt'], rsid['ref']))]
用户也可以在得到pheweb网站上的 rsids-v154-hgXX.tsv.gz 文件(7亿多行)后,在本Github的 scripts文件夹下载本课题组修订的 add_rsid.py,dos2unix add_rsid2.py,然后运行如下示例命令。注意:--sep 后面有双引号。 GWAS文件必须按照CHR和POS排好序,可用sort -k 1,1V -k 2,2n来实现,-V是为了把chrX和chrY排到最后,但是需要把第一行先写到新文件里。
python add_rsid.py -i test.tsv --sep "\t" --chr CHR --pos POS --ref NEA --alt EA -d files/dbsnp/rsids-v154-hg19.tsv.gz -o out.tsv
通过LD的计算来找到GWAS数据里面的independent top hits,也有一些问题。比如,g1k的LD不是金标准,r2也不是最合理的筛选办法,并且计算量很大。如果不考虑 SNP之间的LD,只考虑距离,假设GWAS的第1,2,3 列分别是 SNP, CHR, POS,最后一列是P,可以用下面这个简单的代码来寻找GWAS数据里面每1MB区间的top SNP。
zcat ABC.gwas.gz | awk 'NR==1 || $NF<5e-8 {b=sprintf("%.0f",$3/1e6); print $1,$2,$3,$NF,b}' | \
sort -k 2,2n -k 5,5n -k 4,4g | awk '{if (arr[$NF] !="Y") print $0; arr[$NF] ="Y"}'
要把上述得到的显著区域跟已发表的文章中的SNP进行比较,看是不是有重叠(1MB范围之内的重叠都算),可以用 bedtools。
如果有个体数据数据,可以采用 GSMR,参考GSMR 文章 。 如果没有个体数据,只有别人报道的 exposure 和 outcome 的 BETA 和 SE,就可以使用Bristol大学开发的TwoSampleMR R包或剑桥大学团队开发的MendelianRandomization R包。 前者有非常强大的数据调取和数据管理功能,而后者比较简单直观。
基因注释信息浏览器:
- 非常经典的 dbSNP: https://www.ncbi.nlm.nih.gov/snp/
- 非常经典的 UCSC genome browser: https://www.genome.ucsc.edu/
- 美国精准医学All of Us:https://www.researchallofus.org/ 和 https://databrowser.researchallofus.org/
- 密西根大学公卫学院 TopMed browser: https://bravo.sph.umich.edu/
- 一天发了7篇 NATURE系列文章的Gnomad项目的 browser: https://gnomad.broadinstitute.org/
- 带”Global“的 GlobalBiobankEngine:https://github.com/rivas-lab
理论学习
进化对人类疾病的影响:
- Y 2021. Nature Review Genetics. The influence of evolutionary history on human health and disease
- Y 2023. Science Bulletin. Recent positive selection signatures reveal phenotypic evolution in the Han Chinese population
GWAS-PRS-MR ”三驾马车“ 入门:
GWAS:
- Y 2006. Nature Review Genetics. A tutorial on statistical methods for population association studies
- Y 2014. Nature Protocols. Quality control and conduct of genome-wide association meta-analyses
- Y 2021. Nature Reviews Methods Primers. Genome-wide association studies
- 芬兰赫尔辛基大学 GWAS 课程
PRS:
- Y 2019. Lancet Respiratory Medicine. Identification of risk loci and a polygenic risk score for lung cancer: a large-scale prospective cohort study in Chinese populations
- Y 2020. Nature Protocols. Tutorial: a guide to performing polygenic risk score analyses
- Y 2022. EHJ. A polygenic risk score improves risk stratification of coronary artery disease: a large-scale prospective Chinese cohort study
- Y 2023. Nature Medicine. A multi-ancestry polygenic risk score improves risk prediction for coronary artery disease
MR:
- Y 2012. 经典案例 Lancet. Plasma HDL cholesterol and risk of myocardial infarction: a mendelian randomisation study
- Y 2017. 一篇解读我的PRIMe方法的文章 Ele Zeggini. Statistical methods to detect pleiotropy in human complex traits
- Y 2021. 通过MR进行中介分析 Mendelian randomisation for mediation analysis: current methods and challenges for implementation
- Y 2022. 入门必读 Nature Reviews Methods Primers. Mendelian randomization
- Y 2023. 一篇非常简单的JAMA子刊 JAMA Psychiatry. Association of Genetically Predicted Insomnia With Risk of Sepsis
一些有用、有趣的实用工具:
- The R Graph Gallery
- Doing and reporting your first mediation analysis in R
- Add P-values and Significance Levels to ggplots
- Top 100 R resources on COVID-19 Coronavirus
- 以及 scitb, CanvasXpress, modelSummary, forplo,sankey diagram, CellChat, ComplexHeatmap,等等
- X2Y: X到Y
- Y4x: Y中数据中为(for)了X的部分,里面并没有X的数据,因此小写x。
- X8M:X加上(bus)了M,的并集。
- dat0:0原始大量的数据,最好读一次,不要多次读取。
- dat1:1临时数据,在loop里面使用,要不然在loop外面和里面都用dat,会出问题。
- dat:正常分析用的数据。