Skip to content

Latest commit

 

History

History
212 lines (166 loc) · 12.9 KB

README.md

File metadata and controls

212 lines (166 loc) · 12.9 KB
Raw


#1. 下载和处理国际上公用公开的数据

这是初二生物学课本里面的一页。

middle school

#1.1 HAPMAP3 genotype 数据, 一般作为 LD 计算的 reference panel

打开 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 等软件使用。 通过LD的计算来找到GWAS数据里面的independent top hits，计算量很大。如果不考虑 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"}'

对应的R代码如下🔔🏮

dat <- read.table("GWAS.gz", header=T) %>% filter(P<=5e-08) %>% mutate(mb=ceiling(POS/1e+06))
dat %>% group_by(mb) %>% slice(which.min(P)) %>% ungroup() %>% select("SNP")

#1.2. 千人基因组项目的genotype数据[VCF格式]，一般作为 imputation 的 reference panel.

在千人基因组官网 下载 Phase 3 对应的 VCF 链接，GRCh37版本。 该数据是用GATK平台生成，用的reference genome 来自GATK resource bundle。 有一个文件罗列了每一个样本的人群（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

#1.4. 使用 deepvariant 生成VCF数据.

DeepVariant安装

  1. 远程链接服务器 172.18.42.37:33899
  2. VMware Workstation构建虚拟机，安装Linux Ubuntu 24.04.1 LTS。
  3. 根据deepvariant 官方指南安装Docker版本。
  4. 注意事项：（i）VMware选择个人版本可免费使用；（ii）由于服务器不存在GPU，无需安装CUDA；（iii）拉取docker镜像时，如遇访问缓慢，选择切换镜像源为https://docker.1panel.live。

DeepVariant运行

  1. 打开VMware，登陆Ubuntu系统上的终端进行操作。
  2. 具体步骤参考deepvariant 上的 If you're using GPUs, or want to use Singularity instead, see Quick Start for more details。
  3. 注意事项：如果生成的文件会出现无权限访问，需要在终端修改权限
sudo shown -R administrator /home/administrator/quickstart-output

#2. UKB 基因型和表型数据

UKB 基因型数据已下载到南科大的HPC上。 点击UKB RAP左边的 accessing phenotype data，下载TSV格式的表型数据。下载后，用下面代码将缺失数据替换为NA。

#3. GWAS 运行和结果 “挖掘”

GWAS


#3.1. 专人在服务器上运行

提取需要研究的表型数据和相关的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 文件。


#3.2. 公开的GWAS数据进行练手，或对比

最经典的，起源于美国NIH 的 GWAS Catalog. 这个页面也罗列了一些大型GWAS数据联盟。 欧洲版本，不需要下载就能通过 TwoSampleMR 远程读入。

UKB GWAS 完整的分析结果，网上发布


#3.3. GWAS的管理、QC、注释

对每一个GWAS，首先进行以下“三点”检查：

1. 用 zcat GWAS.gz | awk '{print NF}' | sort -nu | wc -l 文件每一行的列数目是一样的。 一定要扣好第一粒纽扣🤲。可用 sed 's/^\t/NA\t/; s/\t\t/\tNA\t/g; s/\t\t/\tNA\t/g; s/\t$/\tNA/' 解决。
2. 如果文件不是按照CHR和POS排序🏮，pheweb 会报错，可用sort -k 1,1V -k 2,2n。 -V是为了把chrX和chrY排到最后，但是需要把第一行先写到新文件里。
3. GWAS数据本身的问题：
   (1) Allele 最好是大写，awk 和 R 都有 toupper()功能。
   (2) P值最好不要小于1e-312。 要不然，awk 会把其当成0，有一些软件（比如LDSC）也会报错，这个时候要么用Z值，要么人为将这些P值设为1e-300。
   (3) BETA|SE|P出现“三缺一” 的情况，可用： b = se * qnorm(p/2); se = abs(b/qnorm(p/2)); se = (CI_upper - CI_lower)/(1.96*2); p = 2*pnorm(-abs(b/se))

如果生成的GWAS数据超过1G，导致无法上传到pheweb那样的软件，可用下面代码对数据进行瘦身。

zcat $dat.gz | cut -f 2-6,10,13 | awk '{if (NR!=1) {$5=sprintf("%.4f",$5); $6=sprintf("%.4f",$6)} print $0}' | bgzip > $dat.lean.gz

经过QC后的GWAS数据，可用 tabix -f -S 1 -s 1 -b 2 -e 2 GWAS.gz 生成索引文件。 本课题组建议用如下标准的column名称：🐂SNP CHRPOS CHR POS EA NEA EAF N BETA SE Z P🐎。 可用下面的 bash 代码实现：

Arr1=("SNP" "CHR" "POS" "EA" "NEA" "EAF" "N" "BETA" "SE" "P")
Arr2=("snp|rsid|variant_id" "chr|chrom|chromosome" "pos|bp|base_pair" "ea|alt|eff.allele|effect_allele|a1|allele1" "nea|ref|allele0|a2|other_allele|" "eaf|a1freq|effect_allele_freq" "n|Neff" "beta" "se|standard_error" "p|pval|p_bolt_lmm")
dat=XYZ.gz
	head_row=`zcat $dat | head -1 | sed 's/\t/ /g'`; 
	snp=""; chr=""; pos=""; ea=""; nea=""; eaf=""; n=""; beta=""; se=""; p="" 
	for i in ${!Arr1[@]}; do; eval ${Arr1[$i]}=`echo $head_row | tr ' ' '\n' | grep -Einw ${Arr2[$i]} | sed 's/:.*//'`; done
	echo dat $dat, snp $SNP, ea $EA, nea $NEA, n $N, beta $BETA, p $P

对应的R代码如下：

replacement=c('SNP', 'CHR', 'POS','EA', 'NEA', 'EAF', 'N', 'BETA', 'SE', 'P') 
pattern=c('^snp$|^rsid$|variant_id', '^chr$|^chrom', '^bp$|^pos$|^position|^base_pair', '^ea$|^alt$|^a1$|^effect_allele$', '^nea$|^ref|^allele0$|^a2$|^other_allele', '^eaf$|a1freq|effect_allele_freq', '^n$|Neff', '^beta$|^effect$', '^se$|standard_error', '^p$|^pval$|^p_bolt_lmm')
names(dat) <- stringi::stri_replace_all_regex(toupper(names(dat)), pattern=toupper(pattern), replacement=replacement, vectorize_all=FALSE)

最后，可使用密西根大学开发的Pheweb 实现可视化🏮。日本版本pheweb.jp。中国版本的是本课题组建立的 pheweb.cn。 Pheweb有一个强大的add_rsids.py 的功能，但是存在先天缺陷。根据该聊天记录，用户可以在安装pheweb 后找到 add_rsids.py 文件（find /home/ -name "add_rsid*" 或者 pip show --files pheweb），修改一行代码（第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_rsid.py，然后运行如下示例命令。注意:--sep 后面有双引号。

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

密西根大学还开发了locuszoom 实现基因组局部地区的可视化🔍。


#4. Post-GWAS分析

#4.1 单个 GWAS 的深度功能（function）分析

可先尝试傻瓜相机式的FUMA 网上解读系统，见参考文献 FUMA


#4.2 基于单个GWAS的多基因风险评分PRS

相关的方法学，请参考经典版的 PLINK0.9 和新版的 PLINK1.9


#4.3. 两个或多个GWAS之间的 genetic correlation 分析

一般用LDSC软件。 compareB


#4.4. 两个或多个GWAS之间的 Mendelian Randomization 分析

如果有个体数据，可以用 OneSampleMR包。如果只有已发表的summary数据，就可以使用Bristol大学开发的TwoSampleMR R包或剑桥大学团队开发的MendelianRandomization R包。 如果工具变量太多，存在LD问题，TwoSampleMR包有clump()的功能。西湖大学杨剑开发的GSMR 也可以应对LD问题。GSMR的好处是在local电脑上运行，不好之处也是在local电脑上运行。 用户一般用 TwoSampleMR 在服务器上进行 clump，但是其过程和稳定性堪忧。



# 参考文献和网站

基因注释信息浏览器：

GWAS-PRS-MR ”三驾马车“ 入门：

GWAS:

PRS:

MR：


一些有用、有趣的实用工具：