我梳理了GWAS全基因组关联分析的整个流程，并提供了基本的命令，用到的软件包括BWA、samtools、gatk、Plink、Admixture、Tassel等，在此分享出来给大家提供参考。
一、BWA比对
1.构建索引
bwa index -a is example.fasta
#构建索引 -a is算法 （BWT构造算法：bwtsw、is或rb2）
2.进行比对
bwa mem -t 6 -R '@RG\tID:foo\tPL:Illumina\tSM:example'
example.fasta example_1.fq.gz example_2.fq.gz > example.sam
# 进行对比 mem算法 -t 运行的核数目
# -R添加头部
ID：这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；
PL：指的是所用的测序平台，这个信息不要随便写，在GATK中，PL只允许被设置为：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN这几个信息。如果不知道，那么必须设置为UNKNOWN。
SM：样本ID。
LB：测序文库的名字，如果上面的lane ID足够用于区分的话，也可以不用设置LB；
（用GATK检测变异 其中ID,PL和SM信息是必须的）
二、samtools格式转换
1.sam格式转换为bam格式
samtools view -bS example.sam -o example.bam
# -b 输出bam格式文件 -S 输入sam格式文件
2.质控
samtools view -h -b -q30 example.bam > example.q30.bam
# -q 比对的最低质量值 -h 输出的文件包含头部信息 -b 输出bam格式文件
3.构建索引
samtools faidx base/example.fasta
# 该命令会在example.fasta所在目录下创建一个example.fai索引文件
gatk CreateSequenceDictionary -R example.fasta -O example.dict
# 创建gatk索引 生产dict文件
三、gatk变异检测
1.排序
gatk SortSam -I example.q30.bam -O example.q30.sort.bam
-R base/example.fasta -SO coordinate --CREATE_INDEX true
# -I 输入文件 -O 输出文件 -R参考基因组 --CREATE_INDEX 是否建立索引
将sam文件中同一染色体对应的条目按照坐标顺序从小到大进行排序
2.标记重复序列
gatk MarkDuplicates -I example.q30.sort.bam -O example.q30.sort.markdup.bam -M example.q30.sort.markdup_metrics.txt
# -I 输入排序后的文件 -O 输出文件 -M 输出重复矩阵
注意：
samtools index example.q30.sort.markdup.bam
# 为gatk建立索引这里非常重要
3.检测变异
gatk HaplotypeCaller --emitRefConfidence GVCF -R base/example.fasta -I example.q30.sort.markdup.bam -L chrX -O chrX.g.vcf.gz
# HaplotypeCaller同时检测snp和indel -R 参考基因组 -I 输入文件 -L 仅检测该染色体的变异（分染色体检测变异，加快速度）-O 输出文件
合并文件(g.vcf)
gatk CombineGVCFs -R base/example.fasta --variant example1.g.vcf.gz --variant example2.g.vcf.gz -O con.vcf.gz
# -R 参考基因组 --variant 输入变异文件 可以输入多个文件 -O 输出文件
检测变异
gatk GenotypeGVCFs -R ref.fa -V test.g.vcf -O test.vcf
4.提取SNP变异
gatk SelectVariants -R base/example.fasta -V test.vcf -O test.snp.vcf --select-type-to-include SNP
# -R 参考基因组 -O 输出vcf文件 -V 输入vcf文件 --select-type-to-include 选取提取的变异类型(#SNP,MNP,INDEL,SYMBOLIC,MIXED)
5.对SNP进行过滤
gatk VariantFiltration -O chr5.Filt.vcf -V chr5.vcf
--cluster-window-size 10
--filter-expression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" --filter-name "HARD_TO_VALIDATE"
--filter-expression "DP < 5 " --filter-name "LowCoverage"
--filter-expression "QUAL < 30.0 " --filter-name "VeryLowQual"
--filter-expression "QUAL > 30.0 && QUAL < 50.0 "
--filter-name "LowQual" --filter-expression "QD < 1.5 "
--filter-name "LowQD"
# -O 输出文件 -V输入变异文件 --filter-expression 过滤条件 --filter-name 被过滤掉的SNP不会删除，而是给一个标签，例如 "LowCoverage" 。
这里可以将过滤条件合并，仅给出一个标签。--cluster-window-size 以10个碱基为一个窗口
这里通过设定相应的参数值进行了硬过滤，实际应用时还要根据数据类型及自己的需求设定相应的参数。
6.合并文件(vcf)
删除掉被过滤的SNP
grep -v "LowCoverage" Filt.vcf > Filt1.vcf
# -v显示不包含匹配文本的所有行 "LowCoverage"上一步给出的标签
合并成一个文件
gatk MergeVcfs -I chr1.Filt.vcf -I chr2.Filt.vcf -I chr3.Filt.vcf -O Filt.vcf
# -I 输入文件 可以多次指定 -O 输出文件
至此为止就得到了整个群体的VCF变异文件，后续都是基于此文件来进行相应的分析。
四、Plink格式转换及主成分分析
1.VCF格式转换为 ped/map格式
vcftools --vcf snp.vcf --plink --out snp
2.ped/map格式转换为bed/bim/fam格式
plink --file snp --make-bed --out snp_test
3.主成分分析
Plink --threads 8 --bfile snp --pca 10 --out pca
# --threads 线程数 --bfile 输入.bed文件 --pca 主成分的成分数 --out输出的文件名
五、Admixture 群体结构
1.群体结构分析
for K in 2 3 4 5 6 7 8 9 10; \
do admixture --cv hapmap3.bed $K | tee log${K}.out; done
#2 3 4 5 6 7 8 9 10分成的群体结构数 hapmap3.bed 输入文件
注意:
如果你的数据格式是plink的bed文件, 比如a.bed, 那么你应该包含a.bim, a.fam
如果你的数据格式是plink的ped文件, 比如b.ped, 那么你应该包括b.map
K值根据实际情况进行设置，通过比较得到最佳K值。
grep -h CV log*.out
#查看最佳K值 输出最佳K值文件：hapmap3.3.Q
2.R语言作图
tbl=read.table("hapmap3.3.Q")
barplot(t(as.matrix(tbl)), col=rainbow(K),
xlab="Individual", ylab="Ancestry", border=NA)
hapmap.3.Q文件为群体结构的结果，作为协变量进行关联分析
六、Tassel关联分析
Tassel的管道命令不允许有回车符号，使用以下命令时需要将#注释及换行删除。
1.VCF格式转换为 hmp格式
run_pipeline.pl -SortGenotypeFilePlugin -inputFile example.vcf -outputFile example -fileType VCF
#给vcf文件排序，排成tassel认可的序列
#-inputFile 输入的文件名 -outputFile 输出的文件名 -fileType 输出的文件格式
run_pipeline.pl -fork1 -vcf example.vcf -export example -exportType Hapmap -runfork1
#vcf文件转换为hapmap格式
#-vcf 输出的文件 -export 输出的文件 -exportType 输出的文件格式
2.亲缘关系
run_pipeline.pl -importGuess genotype.hmp.txt #打开数据文件
-KinshipPlugin -method Centered_IBS -endPlugin #计算亲缘关系
-export genotype_kinship.txt #输出文件名
-exportType SqrMatrix #输出格式
3.关联分析
混合线性模型
run_pipeline.pl -fork1 -h genotype.hmp.txt -filterAlign -filterAlignMinCount 150 -filterAlignMinFreq 0.05 -filterAlignMaxFreq 1
#-h读取hapmap格式的基因型数据 -filterAlign 打开过滤选项
#-filterAlignMinCount -filterAlignMinFreq -filterAlignMaxFreq 过滤的条件 需要按照实际情况进行修改
-fork2 -r traits.txt
#-r 读取表型数据
-fork3 -q population_structure.txt -excludeLastTrait
#-q 读取群体结构数据 -excludeLastTrait 去掉最后一个群体结构 （当分成的群体结构＞2时）
-fork4 -k kinship.txt
#读取亲缘关系数据
-combine5 -input1 -input2 -input3 -intersect -combine6 -input5 -input4 -mlm -export example
#-combine合并数据 -mlm混合线性模型 -export输出文件名
一般线性模型
run_pipeline.pl
-fork1 -h genotype.hmp.txt -filterAlign
-filterAlignMinCount 150 -filterAlignMinFreq 0.05 -filterAlignMaxFreq 1
-fork2 -r traits.txt
-fork3 -q population_structure.txt -excludeLastTrait
-combine4 -input1 -input2 -input3 -intersect -glm -export example
#-glm 一般线性模型
4.R语言作图
Library(qqman)
#加载qqman包
曼哈顿图
manhattan(example,ylim=c(0,10),col = color_set,annotatePval = 0.01)
# ylim Y轴范围 col 颜色 annotatePval 标记最高位点 CHR==1 绘制每个染色体的曼哈顿图
Q-Q plot
qq(example$P)
七、其他
1.基因组统计工具
可以统计fasta和fastq中的信息。
seqkit fx2tab example.fasta -l -n
-l 统计序列长度 -n 统计染色体
2.提取文本文档中某列
用于Tassel关联分析后的结果文件，提取相应的列进行R语言绘图。
cat MLM.txt | awk '{print $1" "$3" "$4" "$7}' > manhattan.txt
# $提取的列数
3.删除文本文档中不包含匹配文本的行
grep -v "LowCoverage" Filt.vcf > Filt1.vcf
参考：
https://www.jianshu.com/p/0b0c4ab4c38a GATK4全基因组数据分析最佳实践...
http://starsyi.github.io/2016/05/25/%E5%8F%98%E5%BC%82%E6%A3%80%E6%B5%8B%EF%BC%88BWA-SAMtools-picard-GATK%EF%BC%89/ 变异检测（BWA+SAMtools+picard+GATK）
https://www.jianshu.com/p/c41e8f3266b4 GATK4.1 call SNP
网友评论