跟着Nature Genetics学GWAS分析：emmax软件gwas分析/qqman包展示结果

用户7010445

发布于 2023-08-23 10:39:58

5480

发布于 2023-08-23 10:39:58

文章被收录于专栏：小明的数据分析笔记本

论文

Super-pangenome analyses highlight genomic diversity and structural variation across wild and cultivated tomato species

https://www.nature.com/articles/s41588-023-01340-y

西红柿NG_superPan正文.pdf

数据分析的代码

https://github.com/HongboDoll/TomatoSuperPanGenome

论文里提供了绝大部分的数据处理代码，很好的学习材料，今天的推文我们学习一下论文中GWAS分析的相关代码

vcf文件用到的是拟南芥的数据，下载链接

http://1001genomes.org/data/GMI-MPI/releases/v3.1/1001genomes_snp-short-indel_only_ACGTN.vcf.gz

这个数据之前的推文也用过

文献笔记五十四：全基因组关联分析鉴定拟南芥中控制种子大小的调节因子

但是想不起来表型数据是在哪里下载的了

对vcf文件进行过滤

关于vcf文件的操作参考这个链接 https://www.jianshu.com/p/d46d3682637d

vcftools --gzvcf 1001genomes_snp-short-indel_only_ACGTN.vcf.gz --remove-indels --recode --recode-INFO-all --min-alleles 2 --max-alleles 2 --max-missing 0.9 --out snpOnly

删除indel位点只保留二等位变异位点能分型的样本占总样本的比例至少为0.9

~/biotools/plink19/plink --vcf snpOnly.recode.vcf --recode12 --allow-extra-chr --geno 0.1 --mac 5 --allow-no-sex --out at.snp

~/biotools/plink19/plink --allow-extra-chr --file at.snp --indep-pairwise 50 5 0.2 --recode vcf-iid --out at.snp.LDpruned

~/biotools/plink19/plink --allow-extra-chr --file at.snp --recode vcf-iid --extract at.snp.LDpruned.prune.in --out at.snp.LDpruned

~/biotools/plink19/plink --allow-extra-chr --vcf at.snp.LDpruned.vcf --make-bed --out at.snp.LDpruned.bfile

~/biotools/plink19/plink --bfile at.snp.LDpruned.bfile --allow-extra-chr --allow-no-sex --pca 10 --out PCA

cat PCA.eigenvec | awk '{print $1,$2,"1",$3,$4,$5,$6,$7}' > PCA_snp.txt

~/biotools/plink19/plink --file at.snp --allow-extra-chr --recode 12 --output-missing-genotype 0 --transpose --out at_snp

~/biotools/emmax/emmax-kin-intel64 -v -d 10 at_snp

这里表型数据忘记在哪里下载了，自己随便构造吧

head -n 20 snpOnly.recode.vcf | grep '#CHROM' | awk 'gsub("\t","\n")' | tail -n 1135 > pheno.txt

library(tidyverse)

read_delim("D:/R_4_1_0_working_directory/env001/data/20230510/pheno.txt",
           delim = "\t",
           col_names = FALSE) %>%
  mutate(X2=X1,
         X3=sample(seq(0,100,by=0.05),1135)) %>% 
  write_delim("D:/R_4_1_0_working_directory/env001/data/20230510/pheno01.txt",
              delim = "\t",
              col_names = FALSE)

最后格式前两列是样本名，最后一列是表型数据，如果有缺失可以用NA代替分隔符是制表符

image.png

gwas分析

~/biotools/emmax/emmax-intel64 -v -d 10 -t at_snp -p pheno01.txt -k at_snp.aBN.kinf -c PCA_snp.txt -o at_snp

cat at_snp.map | paste - at_snp.ps | awk '{print $2,$1,$4,$8}' | sed '1i\SNP CHR BP P' | sed 's/ /\t/g' > gwas.output

Rscript manhattan_qq.R gwas.output gwas.png 5

manhattan_qq.R 这个脚本是论文中提供的

最后5是显著性的阈值，是自己随便写的，

整个代码能够跑通，但其中有一些细节自己还不是很明白，需要再多看几遍