GWAS全基因组关联分型脚本记录

GWAS分析docker镜像使用





############################################################################################
#北京组学生物科技有限公司
#author huangls
#date 2020.05.06
#version 1.0
#linux 基础:
# https://study.163.com/course/introduction/1006346005.htm?share=1&shareId=1030291076
#更多docker使用及Linux服务器搭建:
# https://study.163.com/course/introduction/1209757831.htm?share=1&shareId=1030291076
###############################################################################################

######################################################################################
#下载基因组GWAS分析docker镜像
#docker pull omicsclass/pop-evol-gwas:v1.1    
#启动docker容器并交互式进入
#docker desktop方法
#docker run --rm -it -m 4G --cpus 1  -v D:\gwas:/work omicsclass/pop-evol-gwas:v1.2
#docker toolbox方法
#docker run --rm -it -m 4G --cpus 1  -v /d/gwas:/work omicsclass/pop-evol-gwas:v1.2
####################################################
# 创建目录,以及准备一些路径变量,方便后续调用
####################################################
#docker run --rm -it -v /share/work/huangls/piplines/omicsclass/pop-evol-gwas/scripts/:/work/scripts -v /share/nas1/huangls/test/gwas_rice:/work omicsclass/pop-evol-gwas:v1.2
workdir=/work/  #设置工作路径
refdir=$workdir/ref
datadir=$workdir/data
scriptdir=$workdir/scripts
tmpdir=$workdir/tmp
export TMPDIR=$tmpdir      # 设置临时目录 防止容器中默认临时目录写满报错
export PATH=$scriptdir:$PATH    #组学大讲堂提供的脚本 添加到环境变量中方便使用
#一些关键文件所在的路径变量
GFF=$refdir/all.gff3
REF=$refdir/all.con.fa
FAI=$refdir/all.con.fa.fai
GROUP=$datadir/pop_group.txt
#水稻GWAS数据来源:https://www.nature.com/articles/ng.3596
#分组:根据自己的分组名称进行拆分
cat $GROUP |grep GROUP1 |cut -f 1 >$datadir/pop1.txt
cat $GROUP |grep GROUP2 |cut -f 1 >$datadir/pop2.txt
####################################################################
#数据过滤
#####################################################################
cd $workdir  #回到工作目录
mkdir 00.filter
cd 00.filter
#对数据进行过滤
#
#过滤1:vcfutils.pl  过滤掉indel附近的snp
# -w INT    SNP within INT bp around a gap to be filtered [3]
# -W INT    window size for filtering adjacent gaps [10]
#
#vcfutils.pl varFilter -w 5 -W 10 "gzip -dc $datadir/rice.raw.vcf.gz|" |gzip - >all.varFilter.vcf.gz
#
#过滤2:vcftools
#--max-missing Exclude sites on the basis of the proportion of missing data 
#(defined to be between 0 and 1, where 0 allows sites that are completely missing 
#and 1 indicates no missing data allowed).
#
#vcftools --gzvcf all.varFilter.vcf.gz --recode --recode-INFO-all --stdout     --maf 0.05  --max-missing 0.7  --minDP 2  --maxDP 1000      \
#    --minQ 30 --minGQ 0 --min-alleles 2  --max-alleles 2 --remove-indels |gzip - > clean.vcf.gz  
#
#排序  , 如果不排序有可能会导致后面tassel的命令出错
#run_pipeline.pl -Xmx30G  -SortGenotypeFilePlugin -inputFile clean.vcf.gz \
#    -outputFile clean.sorted.vcf.gz -fileType VCF
##########################################################
#进化树分析
###########################################################
#cd $workdir  #回到工作目录
#mkdir 01.phylo_tree
#cd 01.phylo_tree
#文件格式转换
#run_pipeline.pl  -Xmx5G -importGuess  $workdir/00.filter/clean.sorted.vcf.gz  \
#    -ExportPlugin -saveAs supergene.phy -format Phylip_Inter
#
#最大似然法构建进化树
#方法1:fasttree 构建进化树
#fasttree -nt -gtr  supergene.phy   >  fasttree.nwk
#
#进化树绘制
#ggtree.r -t fasttree.nwk -f $GROUP -g Group -n fasttree
#进化树手动美化:https://www.omicsclass.com/article/671
#
########################################################
#PCA分析
#########################################################
#cd $workdir  #回到工作目录
#mkdir 02.PCA
#cd 02.PCA
## plink分析PCA
#plink --vcf  $workdir/00.filter/clean.sorted.vcf.gz --pca 10 --out  plink_pca   \
#    --allow-extra-chr --set-missing-var-ids @:#    --vcf-half-call missing
#绘图
#pca_plink_plot.r -i plink_pca.eigenvec -f $GROUP -g Group --name plink_pca
#
#
#########################################################
#群体结构STRUCTURE分析
######################################################
#cd $workdir  #回到工作目录
#mkdir 03.STRUCTURE
#cd 03.STRUCTURE
## filter LD 
#50 10 0.2   50个SNP 窗口  step 10个  r2 0.2  ; 50 5 0.4
#plink --vcf  $workdir/00.filter/clean.sorted.vcf.gz  --indep-pairwise 50 10 0.2 --out ld   \
#    --allow-extra-chr --set-missing-var-ids @:# 
#plink --vcf  $workdir/00.filter/clean.sorted.vcf.gz  --make-bed --extract ld.prune.in  \
#    --out LDfiltered --recode vcf-iid  --keep-allele-order  --allow-extra-chr --set-missing-var-ids @:#  
#
#转换成plink格式
#vcftools --vcf LDfiltered.vcf --plink \
#    --out plink
#转换成admixture要求的bed格式
#plink --noweb --file plink  --recode12 --out admixture \
#     --allow-extra-chr  --keep-allele-order
#
#admixture 群体结构分析
#for k in {2..10};do
#    echo "admixture -j8 -C 0.01 --cv admixture.ped $k >admixture.log$k.out"
#    admixture -j8 -C 0.01 --cv admixture.ped $k >admixture.log$k.out
#done
#绘图展示 
#structure_plot.r  -d ./ -s admixture.nosex 
#确定最佳K,CV值最小时对应的K值为最佳K
#grep "CV error" *out
#
##########################################################
#连锁不平衡分析 LDdecay
########################################################
#
#cd $workdir  #回到工作目录
#mkdir 04.LDdecay
#cd 04.LDdecay
#
#分组:根据自己的分组名称进行拆分
#cat $GROUP |grep Line |cut -f 1 >popid.txt
#
#PopLDdecay -InVCF  $workdir/00.filter/clean.sorted.vcf.gz \
#        -SubPop  popid.txt -MaxDist 500 -OutStat ld.stat
#Plot_OnePop.pl -inFile ld.stat.gz -output ld


###################################################################
#性状数据分析
###################################################################
#基因型填充
cd $workdir  #回到工作目录
mkdir 04.trait
cd 04.trait
Rscript $scriptdir/trait_plot.r  -i $datadir/traits_gapit.txt 
#如果有多年或者多点的数据可以计算BLUP值再做GWAS分析:https://www.omicsclass.com/article/1380
Rscript $scriptdir/blup_year_or_site.r  -i heading_days.rm_outers.txt
#查看 性状数据分布,去掉一些极端值样本或者不去标记为NA
###################################################################
#GWAS分析  数据来源:doi:10.1038/ng.3596
###################################################################
#基因型填充
cd $workdir  #回到工作目录
mkdir 05.impute
cd 05.impute
java  -Xmx4g -Djava.io.tmpdir=$TMPDIR -jar /biosoft/beagle.18May20.d20.jar \
       gt=$workdir/00.filter/clean.sorted.vcf.gz   out=all.impute impute=true window=10 nthreads=2
###################################################################
#利用tassel进行GWAS分析
###################################################################
cd $workdir  #回到工作目录
mkdir 06.gwas_tassel
cd 06.gwas_tassel
tassel_trait=$workdir/04.trait/heading_days_BLUP_value.tsv
#计算PCA和kinship, 也可以使用前面的群体遗传进化分析的结果
run_pipeline.pl -Xmx30G -fork1 -vcf $workdir/00.filter/clean.sorted.vcf.gz \
    -PrincipalComponentsPlugin  -ncomponents 3  -covariance true -endPlugin \
    -export pca -runfork1 
run_pipeline.pl -Xmx30G -vcf $workdir/00.filter/clean.sorted.vcf.gz \
    -KinshipPlugin -method Centered_IBS -endPlugin -export kinship.txt -exportType SqrMatrix 
kinship=$workdir/06.gwas_tassel/kinship.txt
structure=$workdir/06.gwas_tassel/pca1.txt
#也可以用structure计算的Q矩阵,但是需要去掉最后一列
#GLM方法,不考虑亲缘关系的影响
run_pipeline.pl -Xmx30G -fork1 -vcf $workdir/00.filter/clean.sorted.vcf.gz \
    -fork2 -t $tassel_trait -fork3 -q $structure -excludeLastTrait \
   -combine5 -input1 -input2 -input3 -intersect -FixedEffectLMPlugin \
   -endPlugin -export gwas_hd_GLM
cut -f 3,4,6 gwas_hd_GLM1.txt>glm_pvalue.txt
Rscript $scriptdir/gwas_manhattan_plot.r -i glm_pvalue.txt -F $FAI -T heading_days -n heading_days_glm_manhattan -c 1.5e-5
Rscript $scriptdir/qq_plot.r  -i glm_pvalue.txt -n heading_days_glm_qq
#MLM方法
#单个表型数据关联 
run_pipeline.pl -Xmx30g -fork1 -vcf $workdir/00.filter/clean.sorted.vcf.gz \
    -fork2 -r $tassel_trait -fork3 -q $structure -excludeLastTrait \
    -fork4 -k $kinship -combine5 -input1 -input2 -input3 \
    -intersect -combine6 -input5 -input4 -mlm -mlmVarCompEst P3D \
    -mlmCompressionLevel None  -export gwas_hd_MLM -runfork1 -runfork2 -runfork3 -runfork4
cut -f 3,4,7 gwas_hd_MLM2.txt|grep -v "NaN" >mlm_pvalue.txt
Rscript $scriptdir/gwas_manhattan_plot.r -i mlm_pvalue.txt -F $FAI -T heading_days -n heading_days_mlm_manhattan -c 1.5e-5
Rscript $scriptdir/qq_plot.r  -i mlm_pvalue.txt -n heading_days_mlm_qq


#CMLM 
run_pipeline.pl -Xmx30g -fork1 -vcf $workdir/00.filter/clean.sorted.vcf.gz \
    -fork2 -r $tassel_trait -fork3 -q $structure -excludeLastTrait \
    -fork4 -k $kinship -combine5 -input1 -input2 -input3 \
    -intersect -combine6 -input5 -input4 -mlm -mlmVarCompEst P3D \
    -mlmCompressionLevel Optimum  -export gwas_hd_CMLM -runfork1 -runfork2 -runfork3 -runfork4
cut -f 3,4,7 gwas_hd_CMLM2.txt|grep -v "NaN" >cmlm_pvalue.txt
Rscript $scriptdir/gwas_manhattan_plot.r -i cmlm_pvalue.txt -F $FAI -T heading_days -n heading_days_cmlm_manhattan -c 1.5e-5
Rscript $scriptdir/qq_plot.r  -i cmlm_pvalue.txt -n heading_days_cmlm_qq



###################################################################
#利用emmax进行GWAS分析
###################################################################
cd $workdir  #回到工作目录
mkdir 08.gwas_emmax
cd 08.gwas_emmax
## prepare vcf
plink --vcf $workdir/00.filter/clean.sorted.vcf.gz --maf 0.05 --geno 0.1  --recode12  --output-missing-genotype 0 --transpose --out snp_filter   --set-missing-var-ids @:#  --allow-extra-chr
## prepare phenotype
perl $scriptdir/sort_trait.pl snp_filter.tfam ../06.gwas_tassel/hd.txt > hd.sort.txt
#kinship
emmax-kin-intel64 -v -d 10  -o ./pop.kinship  snp_filter
#关联分析
emmax-intel64 -v -d 10 -t snp_filter  -p hd.sort.txt -k pop.kinship   -o emmax.out 1> emmax.log 2>emmax.err
#合并数据
paste  snp_filter.map  emmax.out.ps | awk  'BEGIN{print "CHR\tPOS\tP"}{if($2==$5){print $1"\t"$4"\t"$NF}}'  > emmax.pvalue
Rscript $scriptdir/gwas_manhattan_plot.r -i emmax.pvalue -F $FAI -T heading_days -n heading_days_emmax_manhattan -c 1.5e-5
Rscript $scriptdir/qq_plot.r  -i emmax.pvalue -n heading_days_emmax_qq

###########################emmax +Q################################## 
#注意Q矩阵要去掉最后一列
paste LDfiltered.nosex  admixture.5.Q |awk '{print $1" "$2" 1 "$3" "$4" "$5" "$6}'  > pop.Qmatrix
emmax-intel64 -v -d 10 -t snp_filter  -p trait.sort.txt -k pop.kinship -c pop.Qmatrix   -o emmax.out.Q 1> emmax.log 2>emmax.err
paste  snp_filter.map  emmax.out.Q.ps | awk  'BEGIN{print "CHR\tPOS\tP"}{if($2==$5){print $1"\t"$4"\t"$NF}}'  > emmax.q.pvalue
Rscript $scriptdir/gwas_manhattan_plot.r -i emmax.q.pvalue -F $FAI -T heading_days -n heading_days_emmax.q_manhattan -c 1.5e-5
Rscript $scriptdir/qq_plot.r  -i emmax.q.pvalue -n heading_days_emmax.q_qq
###################################################################
#利用gapit进行GWAS分析
###################################################################
cd $workdir  #回到工作目录
mkdir 07.gwas_gapit
cd 07.gwas_gapit
#vcf 转换成 hmp格式:
run_pipeline.pl -Xms10G -Xmx64G -fork1 -vcf $workdir/00.filter/clean.sorted.vcf.gz \
    -export ./hapmap  -exportType Hapmap
#排序
run_pipeline.pl -Xms10g -Xmx100g  -vcf genotype.filter.vcf -sortPositions -export genotype.filter.hapmap -exportType HapmapDiploid
#关联分析
for m in GLM MLM MLMLM CMLM ECMLM SUPER FarmCPU Blink FaST EMMA EMMAx;do
    echo "gapit_gwas.r -i hapmap.hmp.txt -t hd.txt -o $m -m $m"
    gapit_gwas.r -i hapmap.hmp.txt -t hd.txt -o $m -m $m
done
#可视化绘图
cd GLM
cut -d "," -f 2-4 GAPIT.GLM.heading_days.GWAS.Results.csv|sed 's#,#\t#g'>pvalue.txt
Rscript $scriptdir/gwas_manhattan_plot.r -i pvalue.txt -F $FAI -T heading_days -n heading_days_manhattan -c 1.5e-5
Rscript $scriptdir/qq_plot.r  -i pvalue.txt -n heading_days_qq
###################################################################
#利用GEMMA进行GWAS分析
###################################################################
cd $workdir  #回到工作目录
mkdir 09.gwas_GEMMA
cd 09.gwas_GEMMA
vcftools   --gzvcf $workdir/00.filter/clean.sorted.vcf.gz --plink --out out
plink --file out   --make-bed --out test --noweb
#因为GEMMA是通过识别plink格式的fam文件中的表型来进行关联分析的,所以我们需要预处理一下,
#把表型数据添加进去。添加到fam文件的6,7,8...等等列。
#性状排序
perl $scriptdir/sort_trait.pl test.fam ../06.gwas_tassel/hd.txt |cut -f 3 > hd.sort.txt
#合并文件
cat test.fam | cut -d " " -f 1-5 |paste -d " " - hd.sort.txt >test.fam1
mv -f test.fam1 test.fam
#获得kinship矩阵,使用混合线性模型进行分析。
#-gk 2 标准化的方法计算G矩阵
gemma -bfile test -gk 2 -o kin
#关联分析
gemma -bfile test -k output/kin.sXX.txt -lmm 1 -o hd
#其中:
#chr:SNP所在染色体号
#rs: SNP名称
#ps: SNP物理位置
#n_miss: SNP缺失个体数
#allele1: 次等位基因
#allele0: 主等位基因
#af:SNP频率
#beta: SNP效应值
#se: beta估计标准误
#l_remle: 计算该SNP效应时对应的lamda的remle估计值。
#p_wald :wald检验P值
#其中,我们最关心的三个结果是chr, ps, p_wald,我们可以借助这三个结果画曼哈顿图和QQ图。l_remle比较难理解,需要懂模型才知道它的含义,但对分析来说,不是很重要。

cat hd.assoc.txt|sed 's#S##' |awk -F "[\t_]" 'BEGIN{OFS="\t"}{print $2,$4,$NF}'>pvalue.txt
Rscript $scriptdir/gwas_manhattan_plot.r -i pvalue.txt -F $FAI -T heading_days -n heading_days_manhattan -c 1.5e-5
Rscript $scriptdir/qq_plot.r  -i pvalue.txt -n heading_days_qq




GWAS关联分析课程推荐:https://study.163.com/course/introduction/1212364803.htm?share=1&shareId=1030291076

  • 发表于 2021-12-01 14:10
  • 阅读 ( 1806 )
  • 分类:GWAS

0 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

579 篇文章

作家榜 »

  1. omicsgene 579 文章
  2. 安生水 263 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. 红橙子 78 文章
  6. CORNERSTONE 72 文章
  7. 生信老顽童 52 文章
  8. landy 37 文章