GWAS分析docker镜像使用的方法

这篇文章主要介绍"GWAS分析docker镜像使用的方法"，在日常操作中，相信很多人在GWAS分析docker镜像使用的方法问题上存在疑惑，小编查阅了各式资料，整理出简单好用的操作方法，希望对大家解答"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 /pop-evol-gwas:v1.1    #启动docker容器并交互式进入#docker desktop方法#docker run --rm -it -m 4G --cpus 1  -v D:\gwas:/work /pop-evol-gwas:v1.2#docker toolbox方法#docker run --rm -it -m 4G --cpus 1  -v /d/gwas:/work /pop-evol-gwas:v1.2##################################################### 创建目录，以及准备一些路径变量，方便后续调用#####################################################docker run --rm -it -v /share/work/huangls/piplines//pop-evol-gwas/scripts/:/work/scripts -v /share/nas1/huangls/test/gwas_rice:/work /pop-evol-gwas:v1.2workdir=/work/  #设置工作路径refdir=$workdir/refdatadir=$workdir/datascriptdir=$workdir/scriptstmpdir=$workdir/tmpexport TMPDIR=$tmpdir      # 设置临时目录 防止容器中默认临时目录写满报错export PATH=$scriptdir:$PATH    #组学大讲堂提供的脚本 添加到环境变量中方便使用#一些关键文件所在的路径变量GFF=$refdir/all.gff3REF=$refdir/all.con.faFAI=$refdir/all.con.fa.faiGROUP=$datadir/pop_group.txt#水稻GWAS数据来源：https://www.nature.com/articles/ng.3596#分组：根据自己的分组名称进行拆分cat $GROUP |grep GROUP1 |cut -f 1 >$datadir/pop1.txtcat $GROUP |grep GROUP2 |cut -f 1 >$datadir/pop2.txt#####################################################################数据过滤#####################################################################cd $workdir  #回到工作目录mkdir 00.filtercd 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..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.traitcd 04.traitRscript $scriptdir/trait_plot.r  -i $datadir/traits_gapit.txt #如果有多年或者多点的数据可以计算BLUP值再做GWAS分析：https://www..com/article/1380Rscript $scriptdir/blup_year_or_site.r  -i heading_days.rm_outers.txt#查看 性状数据分布，去掉一些极端值样本或者不去标记为NA####################################################################GWAS分析  数据来源：doi:10.1038/ng.3596####################################################################基因型填充cd $workdir  #回到工作目录mkdir 05.imputecd 05.imputejava  -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_tasselcd 06.gwas_tasseltassel_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.txtstructure=$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_GLMcut -f 3,4,6 gwas_hd_GLM1.txt>glm_pvalue.txtRscript $scriptdir/gwas_manhattan_plot.r -i glm_pvalue.txt -F $FAI -T heading_days -n heading_days_glm_manhattan -c 1.5e-5Rscript $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 -runfork4cut -f 3,4,7 gwas_hd_MLM2.txt|grep -v "NaN" >mlm_pvalue.txtRscript $scriptdir/gwas_manhattan_plot.r -i mlm_pvalue.txt -F $FAI -T heading_days -n heading_days_mlm_manhattan -c 1.5e-5Rscript $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 -runfork4cut -f 3,4,7 gwas_hd_CMLM2.txt|grep -v "NaN" >cmlm_pvalue.txtRscript $scriptdir/gwas_manhattan_plot.r -i cmlm_pvalue.txt -F $FAI -T heading_days -n heading_days_cmlm_manhattan -c 1.5e-5Rscript $scriptdir/qq_plot.r  -i cmlm_pvalue.txt -n heading_days_cmlm_qq####################################################################利用emmax进行GWAS分析###################################################################cd $workdir  #回到工作目录mkdir 08.gwas_emmaxcd 08.gwas_emmax## prepare vcfplink --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 phenotypeperl $scriptdir/sort_trait.pl snp_filter.tfam ../06.gwas_tassel/hd.txt > hd.sort.txt#kinshipemmax-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.pvalueRscript $scriptdir/gwas_manhattan_plot.r -i emmax.pvalue -F $FAI -T heading_days -n heading_days_emmax_manhattan -c 1.5e-5Rscript $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.Qmatrixemmax-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.errpaste  snp_filter.map  emmax.out.Q.ps | awk  'BEGIN{print "CHR\tPOS\tP"}{if($2==$5){print $1"\t"$4"\t"$NF}}'  > emmax.q.pvalueRscript $scriptdir/gwas_manhattan_plot.r -i emmax.q.pvalue -F $FAI -T heading_days -n heading_days_emmax.q_manhattan -c 1.5e-5Rscript $scriptdir/qq_plot.r  -i emmax.q.pvalue -n heading_days_emmax.q_qq####################################################################利用gapit进行GWAS分析###################################################################cd $workdir  #回到工作目录mkdir 07.gwas_gapitcd 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 $mdone#可视化绘图cd GLMcut -d "," -f 2-4 GAPIT.GLM.heading_days.GWAS.Results.csv|sed 's#,#\t#g'>pvalue.txtRscript $scriptdir/gwas_manhattan_plot.r -i pvalue.txt -F $FAI -T heading_days -n heading_days_manhattan -c 1.5e-5Rscript $scriptdir/qq_plot.r  -i pvalue.txt -n heading_days_qq####################################################################利用GEMMA进行GWAS分析###################################################################cd $workdir  #回到工作目录mkdir 09.gwas_GEMMAcd 09.gwas_GEMMAvcftools   --gzvcf $workdir/00.filter/clean.sorted.vcf.gz --plink --out outplink --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.fam1mv -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.txtRscript $scriptdir/gwas_manhattan_plot.r -i pvalue.txt -F $FAI -T heading_days -n heading_days_manhattan -c 1.5e-5Rscript $scriptdir/qq_plot.r  -i pvalue.txt -n heading_days_qq

到此，关于"GWAS分析docker镜像使用的方法"的学习就结束了，希望能够解决大家的疑惑。理论与实践的搭配能更好的帮助大家学习，快去试试吧！若想继续学习更多相关知识，请继续关注网站，小编会继续努力为大家带来更多实用的文章！