这篇文章主要介绍“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.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.亿速云.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.亿速云.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分析docker镜像使用的方法”的学习就结束了，希望能够解决大家的疑惑。理论与实践的搭配能更好的帮助大家学习，快去试试吧！若想继续学习更多相关知识，请继续关注亿速云网站，小编会继续努力为大家带来更多实用的文章！
免责声明：本站发布的内容（图片、视频和文字）以原创、转载和分享为主，文章观点不代表本网站立场，如果涉及侵权请联系站长邮箱：is@yisu.com进行举报，并提供相关证据，一经查实，将立刻删除涉嫌侵权内容。
