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.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鏡像使用的方法”的學習就結束了，希望能夠解決大家的疑惑。理論與實踐的搭配能更好的幫助大家學習，快去試試吧！若想繼續學習更多相關知識，請繼續關注億速云網站，小編會繼續努力為大家帶來更多實用的文章！