最近跑通了一遍GWAS分析，全程在linux操作，虽然具体还有好多需要微调的地方，先把代码整理分享出来mark一下
前期准备
1.理论知识
强烈推荐百迈客云课堂课程GWAS生物信息培训课程
或者可以看看我的gwas相关文章
GWAS基本分析内容 （课程学习笔记）
常用GWAS统计方法和模型简介 （课程学习笔记）
临床生物信息学中的GWAS分析 （内附扩展阅读）
精细定位——降低 GWAS的复杂度 （文献研读）
2.数据下载
如果你没有自己的数据又想做gwas分析的话，可以选择3000水稻基因组的http://snp-seek.irri.org/数据库直接下载，vcf、表型数据甚至plink bed/bim/fam文件直接下载，gwas结果也做到了可视化
另外推荐华农谢为博团队开发的这个网站 http://ricevarmap.ncpgr.cn/v2/ 非常好用，提供数据下载和gwas可视化结果）
基因型数据可以根据bioproject accession编号从NCBI上下载：
表型数据直接从网站上以excel或csv格式导出：
533份样品下载和比对挺耗时和占内存的，建议保留bam文件，其他的包括fastq（由bam能转成fq）、中间文件都可以删掉，数据过滤质控很快，整个项目做下来大概耗时一个多月吧
第一步 SNP calling
需要安装的软件：BWA和GATK/samtools
一、BWA比对
1.构建index
bwa index -a is ref.fa 或bwa index -a bwtsw ref.fa (>2G)
samtools faidx ref.fa
java -jar $picard/CreatSequenceDictionary.jar R=ref.fa O=ref.dict
2.每个样分别比对到参考基因组
bwa mem -t 5 -M -R "@RG\tID:A\tSM:A" ref.fa A1_1.fq A1_2.fq > A1.sam &
bwa mem -t 5 -M -R "@RG\tID:A\tSM:A" ref.fa A2_1.fq A2_2.fq > A2.sam & 以此类推.......（-M 将shorter split hits标记为次优，可以兼容Picard．-R 每个标记号需不同，方便后面合并）
3.SortSam
java -jar $picard/SortSam.jar I=A.sam O=A1.sort.bam SO=coordinate
4.MarkDuplicates
java -jar $picard/MarkDuplicates.jar I=A1.sort.bam O=A1.Mark.bam M=A1.metrics
二、SNP检测
1.RealignerTargetCreator
java -jar $GATK -R $ref -T RealignerTargetCreator -I A1.Mark.bam -o A1.realign.interval_list
2.IndelRealigner
java -jar $GATK -R $ref -T IndelRealigner -I A1.Mark.bam -o A1.realn.bam -targetIntervals A1.realign.interval_list
3.HaplotypeCaller
java -jar $GATK -T HaplotypeCaller -R $ref -ERC GVCF -I A1.realn.bam --variant_index_type LINEAR --variant_index_parameter 128000 -o A1.gvcf
4.CombineGVCFs
java -jar $GATK -T CombineGVCFs -R $ref --disable_auto_index_creation_and_locking_when_reading_rods --variant A1.gvcf A2.gvcf A3.gvcf .... -o combine.gvcf (这一步需要把每个样品的gvcf合并)
5.GenotypeGVCFs
java -jar $GATK -T GenotypeGVCFs -nt 4 -R $ref --disable_auto_index_creation_and_locking_when_reading_rods -o test_final.vcf --variant combine.gvcf
第二步 基因型填补
需要安装的软件：Tassel/beagle等
1. tassel (-LDKNNilmputationPlugin参数有误，没有跑成功, 有朋友指出LDKNNi后面是大写的i，不是小写的L,以后再试试)
perl /home/user/soft/tassel_v5/run_pipeline.pl -Xms512m -Xmx5g -importGuess test_final.vcf -LDKNNiImputationPlugin -highLDSSites 50 -knnTaxa 10 -maxLDDistance 100000000 -endPlugin -export test.imputed.vcf -exportType VCF
也可以java -jar sTASSEL.jar 在窗口操作 记得给服务器接显示屏
2.beagle
java -jar beagle.08Jun17.d8b.jar gt=test_final.vcf out=test.imputed.vcf
第三步 数据筛选及格式转换
需要安装的软件：plink等
1.按MAF>0.05和缺失率<0.1过滤
/home/user/soft/plink --vcf test.imputed.vcf --maf 0.05 --geno 0.1 --recode vcf-iid --out test.filter --allow-extra-chr （非数字染色体号ChrUn/Sy用此参数, 建议尽量把染色体号转成数字，另外需要对vcf中的标记ID进行编号）
2.对标记进行LD筛选
/home/user/soft/plink --vcf test.filter.vcf --indep-pairwise 50 10 0.2 --out test.filter --allow-extra-chr （.in文件里是入选的标记id）
3.提取筛选结果
/home/user/soft/plink --vcf test.filter.vcf --make-bed --extract test.filter.in --out test.filter.prune.in
4.转换成structure/admixture格式
/soft/plink --bfile test.filter.prune.in --recode structure --out test.filter.prune.in #生成. recode.strct_in为structure输入格式
/soft/plink --bfile test.filter.prune.in --recode 12 --out test.filter.prune.in #生成.ped为admixture输入格式
第四步 群体结构
需要安装的软件：structure/admixture等
这里我选了比较简单admixture来做 k值范围1到13
/soft/admixture_linux-1.3.0/admixture --cv test.filter.ped 1 >>log.txt
/soft/admixture_linux-1.3.0/admixture --cv test.filter.ped 2 >>log.txt
.......
/soft/admixture_linux-1.3.0/admixture --cv test.filter.ped 13 >>log.txt
wait
grep "CV error" log.txt >k_1to13
取CV error最小时的k值=10, 其中test.filter.prune.in.10.Q结果文件作为关联分析的输入源文件（去掉最后一列 添加表头和ID）
第五步 亲缘关系/PCA(选做)
需要安装的软件: tassel等（ PCA分析可以用R包，已经做了群体结构这里就没做PCA分析）
perl /tassel_v5/run_pipeline.pl -importGuess test_impute.vcf -KinshipPlugin -method Centered_IBS -endPlugin -export test_kinship.txt -exportType SqrMatrix
第六步 关联分析
需要安装的软件: tassel/GAPIT/FaSt-LMM等（ GAPIT很强大，要装很多R包，自动做图可视化。FaSt-LMM以后打算尝试一下）
输入文件的格式需要手动修改一下，也比较简单 (如图)
test.best_k10.txt
test.trait.txt
test_kinship.txt
1.vcf转hapmap格式
perl /tassel_v5/run_pipeline.pl -fork1 -vcf test.imputed.vcf -export test -exportType Hapmap -runfork1
2.SNP位点排序
perl /tassel_v5/run_pipeline.pl -SortGenotypeFilePlugin -inputFile test.hmp.txt -outputFile test_sort -fileType Hapmap
得到test.hmp.txt
3. GLM模型
perl /tassel_v5/run_pipeline.pl -fork1 -h test_sort.hmp.txt -fork2 -r test.trait.txt -fork3 -q test.best_k10.txt -excludeLastTrait -combine4 -input1 -input2 -input3 -intersect -glm -export test_glm -runfork1 -runfork2 -runfork3
4.MLM模型
perl /tassel_v5/run_pipeline.pl -fork1 -h test_sort.hmp.txt -fork2 -r test.trait.txt -fork3 -q test.best_k10.txt -excludeLastTrait -fork4 -k test_kinship.txt -combine5 -input1 -input2 -input3 -intersect -combine6 -input5 -input4 -mlm -mlmVarCompEst P3D -mlmCompressionLevel None -export test_mlm -runfork1 -runfork2 -runfork3
最后得到关联分析的结果文件，喜大普奔！
当然GWAS分析需要根据实际项目材料的需要，灵活地选择分析方法，理解统计学、群体遗传学等原理，认识GWAS的特点和局限性还是很有必要的。这里只是简单地跑通了一遍常用的流程，有很多不懂的地方需要多学习，毕竟大牛都是自己手动写脚本来分析，不怎么用软件的。。/(ㄒoㄒ)/
网友评论