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
- 阅读 ( 2410 )
- 分类：GWAS
你可能感兴趣的文章
- 暑期线下培训班开始报名啦！提前报名享大优惠！ 170 浏览
- GWAS，TWAS，eQTL,EWAS等的关系定义 2700 浏览
- GWAS模型描述方法 943 浏览
- 基因（加性、显性、上位性）效应 5205 浏览
- 群体遗传进化与GWAS分析动植物文献 2342 浏览
相关问题
- 你好，我想请问一下qq图做出来的这种结果是不是不太好呢，横坐标超过4之后突然就变得很陡峭，这大概是什么原因呢？ 1 回答
- LD为什么突然终止呢？自动就killed 1 回答
- 请问我想用GWAS的自然群体进行KAKS分析，请问有什么具体的操作流程吗？ 0 回答
- BSR-seq 结果中的曼哈顿图怎么看？ 0 回答
- 想问一下，RNA-seq测序之后，进行snpcalling得到的一些snp位点是杂合的，关联分析时按缺失小于0.2,第二等位基因频率大于0.05过滤后，得到的显著位点正好是杂合位点。想问一下，进行关联分析时，怎么处理的分子数据中的snp杂合位点，是直接当成缺失，还是过滤掉，还是基因型填充，还是直接导入tassel进行关联分析？谢谢各位老师 1 回答
- GWAS分析模型的选择 1 回答
0 条评论
请先 登录 后评论
