注意：版本不同，命令会不一致。一定要用对应的版本。

需要登录google账户 GATK3的下载地址 GATK4 各个版本之间，有些命令会有变化。 1.GATK的安装、使用别人的教程腾讯云的教程 GATK4数据预处理变异检测(BWA+SAMtools+picard+GATK) bwa+samtools+picardtools+GATK call SNP 流程 因我们服务器渣渣的网络问题，内容都是下载到本地的win10之后，再上传到服务器上。 服务器配置：cpu 4个6核2线程=48 内存：396G gatk是Java程序，下载到本地后解压缩即可使用。在win10使用IDM下载gatk4.0.10.1地址存放目录

/home/chaim/disk/gatk/

unzip gatk4.0.10.1
cd gatk4.0.10.1/
chmod 777 gatk
./gatk --list          //显示gatk的所有子命令

2.GATK4.0.10.1简介常用的pipeline有5种 1. Germline SNPs + Indels
种系SNP+Indel

Somatic SNVs + Indels 体细胞单碱基突变
RNAseq SNPs + Indels
Germline CNVs 种系拷贝数变异 (Copy numbervariations, CNVs)主要指大于1kb 以上的DNA片段的缺失、插入、重复等。一般是结构性变异
Somatic CNVs 体细胞拷贝数变异 1、2、4、5适合DNA测序分析，3适合RNA测序分析。’ 官方文档
开始分析

GATK4.0全基因组和外显子组分析实战

软件： - fastqc检测质量 - fastq/trimmomatic质控 - bwa比对 - samtool格式转换

数据存放位置

所有数据环境前提是在/home/chaim/disk/BSA/目录该目录文件有

119-8-1 //119-8测序原始数据1
-A23-16551278-1279119-8_combined_R1.fastq.gz
-A23-16551278-1279119-8_combined_R2.fastq.gz
119-8-2 //119-8测序原始数据2
-A23-16551278-1279-119-8_combined_R1.fastq.gz
-A23-16551278-1279-119-8_combined_R2.fastq.gz
origin 
- B17SF2447-20_L1_358358.R1.clean.fastq_2.gz  
- B17SF2447-20_L2_358358.R1.clean.fastq.gz
- B17SF2447-20_L1_358358.R2.clean.fastq.gz
- B17SF2447-20_L2_358358.R2.clean.fastq.gz
//四个原始数据

1. 质控检测

fastqc *.fastq.gz -t 8 -o fastqc_out/

安装fastp

wget http://opengene.org/fastp/fastp
chmod 755 fastp

使用fastp质控数据

~~据传，fastp比trimmomatic速度快，效果好。姑且信之。

./fastp -i in.R1.fq.gz -o out.R1.fq.gz -I in.R2.fq.gz -O out.R2.fq.gz

#运行目录于/BSA/119-8   质控119-8的数据
../origin/fastp -i ../119-8-1/A23-16551278-1279119-8_combined_R1.fastq.gz -o ./fastp_out/119-8_1.R1.clean.fastq.gz -I ../119-8-1/A23-16551278-1279119-8_combined_R2.fastq.gz -O ./fastp_out/119-8_1.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &
 ../origin/fastp -i ../119-8-2/A23-16551278-1279-119-8_combined_R1.fastq.gz -o ./fastp_out/119-8_2.R1.clean.fastq.gz -I ../119-8-2/A23-16551278-1279-119-8_combined_R2.fastq.gz -O ./fastp_out/119-8_2.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &

#运行于/BSA/origin/
./fastp -i ./B17SF2447-20_L1_358358.R1.clean.fastq.gz -o ./fastp_out/2447-20_L1$.R1.clean.fastq.gz -I ./B17SF2447-20_L1_358358.R2.clean.fastq.gz -O ./fastp_out/2447-20_L1.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &

./fastp -i ./B17SF2447-20_L2_358358.R1.clean.fastq.gz -o ./fastp_out/2447-20_L2.R1.clean.fastq.gz -I ./B17SF2447-20_L2_358358.R2.clean.fastq.gz -O ./fastp_out/2447-20_L2.R2.clean.fastq.gz -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6 &

fastp参数参考地址 -i R1输入双端测序数据的R1端 -o outputR1质控后输出的R1端 -I R2输入R2原始测序数据 -O outputR2质控后输出的R2端 -Q禁用质量过滤 --thread=5设置线程数为5 --length_required=50设置过滤的最短的序列长度50bp --n_base_limit=6一个reads中N的次数大于6，则舍弃该reads --compression=6输出的gzip文件压缩程度为6，1-9，压缩程度加大。

bwa流程参数 bwa参考1 bwa参考2

2.1. bwa建立索引文件

#/bwa的命令一定不要使用nohup。nohup 的输出信息会被bwa输出到目标文件，会影响后续步骤/ B73序列地址位置

/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa /BSA/bwa/zm437软连接到上述文件

 //工作目录/BSA/bwa/
bwa index -a bwtsw -p zm437 /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa
 ```
`index -a  bwtsw`设置模式，适合大基因组
`-p zm437`设置输出文件名
################分割线#####################################
>/**#注意：此处的2.2.1和2.2.2这两步的bwa 一定要用同一版本的bwa，不然后面会报错**/
#*2.2和2.3*二选一即可，建议使用2.3. bwa的mem的效率更高，且更加准确。
#2.2 bwa寻找输入reads文件的SA坐标

//工作目录/BSA/bwa/ bwa aln zm437 read1.fq.gz -l 30 -k 2 -t 8 -I > read1.fq.gz.sai bwa aln zm437 read2.fq.gz -l 30 -k 2 -t 8 -I > read2.fq.gz.sai

//前4个是本次2447-20样品 bwa aln zm437 ../origin/fastp_out/2447-20_L1.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L1.R1.fq.gz.sai & bwa aln zm437 ../origin/fastp_out/2447-20_L1.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L1.R2.fq.gz.sai & bwa aln zm437 ../origin/fastp_out/2447-20_L2.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L2.R1.fq.gz.sai & bwa aln zm437 ../origin/fastp_out/2447-20_L2.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >2447-20_L2.R2.fq.gz.sai &

//后4个是119-8的数据 bwa aln zm437 ../119-8/fastp_out/119-8_1.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_1.R1.fq.gz.sai & bwa aln zm437 ../119-8/fastp_out/119-8_1.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_1.R2.fq.gz.sai & bwa aln zm437 ../119-8/fastp_out/119-8_2.R1.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_2.R1.fq.gz.sai & bwa aln zm437 ../119-8/fastp_out/119-8_2.R2.clean.fastq.gz -l 30 -k 2 -t 8 -I >119-8_2.R2.fq.gz.sai &

#2.2.2 sai转sam

bwa sampe -r "@RG\tID:\tLB:\tSM:\tPL:ILLUMINA" read1.fq.gz.sai read2.fq.gz.sai read1.fq.gz read2.fq.gz > read.sam

*注释*：`SAMPLE_NAME`应替换为对应样品名称，否则会被当做一个样品处理。

//2447-20数据 bwa sampe zm437 -r "@RG\tID:2447-20\tLB:B73\tSM:2447-20_L1\tPL:ILLUMINA" 2447-20_L1.R1.fq.gz.sai 2447-20_L1.R2.fq.gz.sai ../origin/fastp_out/2447-20_L1.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L1.R2.clean.fastq.gz >2447-20_L1.sam & bwa sampe zm437 -r "@RG\tID:2447-20\tLB:B73\tSM:2447-20_L2\tPL:ILLUMINA" 2447-20_L2.R1.fq.gz.sai 2447-20_L2.R2.fq.gz.sai ../origin/fastp_out/2447-20_L2.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L2.R2.clean.fastq.gz >2447-20_L2.sam &

//119-8数据 bwa sampe zm437 -r "@RG\tID:119-8\tLB:B73\tSM:119-8_1\tPL:ILLUMINA" 119-8_1.R1.fq.gz.sai 119-8_1.R2.fq.gz.sai ../119-8/fastp_out/119-8_1.R1.clean.fastq.gz ../119-8/fastp_out/119-8_1.R2.clean.fastq.gz >119-8_1.sam & bwa sampe zm437 -r "@RG\tID:119-8\tLB:B73\tSM:119-8_2\tPL:ILLUMINA" 119-8_2.R1.fq.gz.sai 119-8_2.R2.fq.gz.sai ../119-8/fastp_out/119-8_2.R1.clean.fastq.gz ../119-8/fastp_out/119-8_2.R2.clean.fastq.gz >119-8_2.sam &

#2.3 BWA的mem的使用，好用快速一步到位。~~[参考地址](https://blog.csdn.net/u013553061/article/details/53120973)~~(注意2.2和2.3，二选一即可，建议使用2.3)

bwa mem的使用

/工作目录在/home/chaim/disk/BSA/bwa// /zm437是B73基因组序列/ /比对的参数-R一定不能省略或写错/ bwa mem -t 24 -M -P -R '@RG\tID:2447-20\tSM:2447-20\tLB:WES\tPL:Illumina' zm437 ../origin/fastp_out/2447-20_L1.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L1.R2.clean.fastq.gz >2447-20_L1.sam & bwa mem -t 24 -M -P -R '@RG\tID:2447-20\tSM:2447-20\tLB:WES\tPL:Illumina' zm437 ../origin/fastp_out/2447-20_L2.R1.clean.fastq.gz ../origin/fastp_out/2447-20_L2.R2.clean.fastq.gz >2447-20_L2.sam & bwa mem -t 12 -M -P -R '@RG\tID:119-8\tSM:119-8\tLB:WES\tPL:Illumina' zm437 ../119-8/fastp_out/119-8_1.R1.clean.fastq.gz ../119-8/fastp_out/119-8_1.R2.clean.fastq.gz >119-8_1.sam & bwa mem -t 8 -M -P -R '@RG\tID:119-8\tSM:119-8\tLB:WES\tPL:Illumina' zm437 ../119-8/fastp_out/119-8_2.R1.clean.fastq.gz ../119-8/fastp_out/119-8_2.R2.clean.fastq.gz >119-8_2.sam &

#3. 对Sam文件进行重排序(recorder)
下载安装最新版[picard](https://github.com/broadinstitute/picard/releases/download/2.18.15/picard.jar)
保存到路径`/home/chaim/disk/BSA/bwa/`
在该路径运行`java -jar  picard.jar -h`,会列出picard包含的所有工具。
3.1 构建索引序列
`nohup samtools faidx zm437 &`
3.2对Sam文件进行重排序

java -jar picard.jar CreateSequenceDictionary R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa O=zm437.dict java -jar picard.jar ReorderSam I=2447-20_L1.sam O=2447-20_L1.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa & java -jar picard.jar ReorderSam I=2447-20_L2.sam O=2447-20_L2.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa & java -jar picard.jar ReorderSam I=119-8_1.sam O=119-8_1.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa java -jar picard.jar ReorderSam I=119-8_2.sam O=119-8_2.reordered.sam R=/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa &

#4.将sam文件转换成bam文件。

samtools view --threads 24 -bS 2447-20_L1.reordered.sam -o 2447-20_L1.bam & samtools view --threads 24 -bS 2447-20_L2.reordered.sam -o 2447-20_L2.bam & samtools view --threads 8 -bS 119-8_1.reordered.sam -o 119-8_1.bam samtools view --threads 8 -bS 119-8_2.reordered.sam -o 119-8_2.bam &

#5. 对bam文件进行sort排序

java -jar picard.jar SortSam INPUT=2447-20_L1.bam OUTPUT=2447-20_L1.sort.bam SORT_ORDER=coordinate & java -Xmx48G -jar picard.jar SortSam INPUT=2447-20_L2.bam OUTPUT=2447-20_L2.sort.bam SORT_ORDER=coordinate & java -Xmx96G -jar picard.jar SortSam INPUT=119-8_1.bam OUTPUT=119-8_1.sort.bam SORT_ORDER=coordinate java -jar picard.jar SortSam INPUT=119-8_2.bam OUTPUT=119-8_2.sort.bam SORT_ORDER=coordinate &

#6. Merge

\合并一个样本的多个lane的bam文件。 java -jar picard.jar MergeSamFiles I=2447-20_L1.sort.bam I=2447-20_L2.sort.bam O=2447-20.bam & java -jar picard.jar MergeSamFiles I=119-8_1.sort.bam I=119-8_2.sort.bam O=119-8.bam

#7. Duplicates Marking
一般情况是RNA-seq不需要这一步，chip-seq,dap-seq,call snp需要这一步骤。
[具体原理](https://www.jianshu.com/p/5781e7d74c40)
测序原理是随机打断，那么理论上出现两条完全相同的read的概率是非常低的，而且建库时PCR扩增存在偏向性，因此标出完全相同的read。

java -jar picard.jar MarkDuplicates REMOVE_DUPLICATES= false MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 INPUT=2447-20.bam OUTPUT=2447-20.repeatmark.bam METRICS_FILE=2447-20.bam.metrics java -jar picard.jar MarkDuplicates REMOVE_DUPLICATES= false MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 INPUT=119-8.bam OUTPUT=119-8.repeatmark.bam METRICS_FILE=119-8.bam.metrics

#8. 生成上一步的结果的索引文件

samtools index 2447-20.repeatmark.bam samtools index 119-8.repeatmark.bam



/*#因前面的bwa的mem的R参数，我第一次运行时未设置完整，导致此处需要二次更改头文件**/
使用picard更改头文件
ID str：输入reads集ID号；LB：read集文库名；PL：测序平台（illunima或solid）；PU：测序平台下级单位名称（run的名称）；SM：样本名称。

java -Xmx96g -jar picard.jar AddOrReplaceReadGroups I=2447-20.repeatmark.bam O=2447-20.repeat.bam LB=lib2447-20 PL=illumina PU=2447-20 SM=2447-20 & java -Xmx96g -jar picard.jar AddOrReplaceReadGroups I=119-8.repeatmark.bam O=119-8.repeat.bam LB=lib119-8 PL=illumina PU=119-8 SM=119-8 &

/*一定不要手动加头文件，手动的后续无法识别。*/
#9.Base (Quality Score) Recalibration
Tools involved: BaseRecalibrator, Apply Recalibration, AnalyzeCovariates (optional)
[参考地址](https://www.jianshu.com/p/0e6162104294)
[流程参考地址](http://kuaibao.qq.com/s/20180414G0F47900?refer=cp_1026)
碱基质量分数重校准（Base quality score recalibration，BQSR)，就是利用机器学习的方式调整原始碱基的质量分数。它分为两个步骤:

- 利用已有的snp数据库，建立相关性模型，产生重校准表( recalibration table)
- 根据这个模型对原始碱基进行调整，只会调整非已知SNP区域。
参数列表
-R ： 参考基因组
-I ： 输入的BAM文件
--known-sites 已知SNP的vcf文件
-O ： 输出的重校准表

java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar BaseRecalibrator -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 2447-20.repeat.bam --known-sites /home/guo/maize/zm437/zea_mays_vcfsort.vcf -O 2447-20_recal_data.table & java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyBQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 2447-20.repeat.bam -bqsr 2447-20_recal_data.table -O 2447-20_recal_reads.bam &

java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar BaseRecalibrator -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 119-8.repeat.bam --known-sites /home/guo/maize/zm437/zea_mays_vcfsort.vcf -O 119-8_recal_data.table & java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyBQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 119-8.repeat.bam -bqsr 119-8_recal_data.table -O 119-8_recal_reads.bam &

*#检测上述生成的bam文件是否可用。*

java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ValidateSamFile -I 2447-20_recal_reads.bam java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ValidateSamFile -I 119-8_recal_reads.bam

如果显示`no errors found`,则可以用HaplotypeCaller call SNP/Indel.






#*二、GATK变异检测*
[参考教程地址](https://www.jianshu.com/p/6f3198b7a070)
**说明：后续会有部分命令有`-nt 24`这个参数，代表使用24个进程。并不是每一个命令都可以开多进程的，需要到gatk官网查询文档，搜索命令后，在命令的API文档里搜索`thread`即可快速查找是否能使用多线程**
（这一步在V4.1.4.0中，java的内存分配最好不要超过16g,建议8g足够。hmm进程也少设置点4-8个。如果Java内存设置太高，某些服务器**即使在服务器内存范围也会出错**会因为内存分配不足而中途挂掉程序，如果中途挂掉，同样有修复办法。[picard修复中断的gvcf文件](https://www.jianshu.com/p/4e9cf85db722)）
1.生成raw vcf 文件（3-5天能跑完的，都是大服务器）
参数说明
先用`HaplotypeCaller`生成gvcf文件，然后再运行`CombineGVCFs`。
java -Xmx96G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar \ #Xmx96G  使用的最大内存
    HaplotypeCaller \   #使用HaplotypeCaller模式，比较吃配置
    -R /home/chaim/disk/BSA/bwa/zm437 \ #参考B73基因组
    -I 2447-20.repeatmark.bam \ #若多样品，则-I sample1.bam -I sample2.bam
    --dbsnp zm437vcf \  #参考B73的snp
    -stand_emit_conf 10 \
    -stand_call_conf 30 \
    -O 2447-20.vcf
1.1生成Gvcf文件（此步非常浪费时间，需要多（3-5）个昼夜）

java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar HaplotypeCaller -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 119-8_recal_reads.bam --dbsnp /home/guo/maize/zm437/zea_mays_vcfsort.vcf --native-pair-hmm-threads 96 -stand-call-conf 30 -O 119-8.g.vcf.gz -ERC GVCF

java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar HaplotypeCaller -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -I 2447-20_recal_reads.bam --dbsnp /home/guo/maize/zm437/zea_mays_vcfsort.vcf --native-pair-hmm-threads 96 -stand-call-conf 30 -O 2447-20.g.vcf.gz -ERC GVCF

1.2合并之前生成的GVCF文件到一个文件。

java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar CombineGVCFs -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa --dbsnp /home/guo/maize/zm437/zea_mays_vcfsort.vcf --variant 2447-20.g.vcf.gz --variant 119-8.g.vcf.gz -O cohort.g.vcf.gz

1.3 GenotypeGVCFs

java -Xmx128G -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar GenotypeGVCFs -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V cohort.g.vcf.gz -O common.vcf.gz

2.select SNP

java -Xmx96g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar SelectVariants -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -O common_SNP.vcf --variant common.vcf.gz --select-type-to-include SNP 2>select_snp.err

3.select indel

java -Xmx96g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar SelectVariants -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -O common_INDEL.vcf --variant common.vcf.gz --select-type-to-include INDEL 2>select_indel.err

/*#4.1和4.2变异过滤，是不同算法的过滤。4.1是机械参数过滤，4.2是机器学习过滤。二者选一即可*/
4.1 filter SNP（变异过滤，硬过滤。）[参数讲解](https://www.jianshu.com/p/938d362fc48d)

java -Xmx4g -jar $GATK -R $REF -T VariantFiltration --variant $Slect_SNP --clusterSize 4 --clusterWindowSize 10 --maskName aroundIndel --mask $Slest_INdel -maskExtend 3 --filterName "lowMQRankSum" --filterExpression "MQRankSum < -12.5" --filterName "highFS" --filterExpression "FS > 60.0" --filterName "lowReadPosRankSum" --filterExpression "ReadPosRankSum < -8.0" --filterName "lowMQ" --filterExpression "MQ < 40.0" --filterName "lowQD" --filterExpression "QD < 2.0" --out $Filter_SNP --genotypeFilterName "lowDP" --genotypeFilterExpression "DP < 8.0" >filter_snp.err

java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar VariantFiltration -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa --variant common_SNP.vcf --cluster-size 4 --cluster-window-size 10 --mask-name aroundIndel --mask common_INDEL.vcf -mask-extension 3 --filter-name "lowMQRankSum" --filter-expression "QUAL < 30" --filter-name "qua130" --filter-expression "MQRankSum < -12.5" --filter-name "highFS" --filter-expression "FS > 60.0" --filter-name "lowReadPosRankSum" --filter-expression "ReadPosRankSum < -8.0" --filter-name "lowMQ" --filter-expression "MQ < 40.0" --filter-name "lowQD" --filter-expression "QD < 2.0" -O common_filtration.vcf --genotype-filter-name "lowDP" --genotype-filter-expression "DP < 8.0" >filter_snp.err

参数详解（[参考](https://mp.weixin.qq.com/s?__biz=MzAxOTUxOTM0Nw==&mid=2649798455&idx=1&sn=67a7407980a57ce138948eb46992b603&chksm=83c1d52bb4b65c3dde31df94e9686654bf616166c7311b531213ebf0010f67a32ce827e677b1&scene=21#wechat_redirect)）
######以下的过滤参数可以根据自己的实际vcf的值分布，自己修改。
QD,描述单位深度的变异值，越大可信度越高。一般过滤掉<2的值。
FS，描述正负链特异性，差异性较大，说明测序或组装的过程中不够随机。FS越小越好。一般过掉掉>40（严格）或60（普通）
MQ 使用bwa-mem的话，正常值应该是60，描述某个位点测序reads的质量值的离散程度。MQ< 40.0
MQRankSum < -12.5
SOR，也是表示正负链特异性，正常值在0-3，过滤掉>3的值。


4.2 变异质控VQSR,共分为两步（##此步本实验不适用，未运行。）
#/*本此实验不能使用该模型过滤，该模型适应于多样本的vcf质控*/
- 1. VariantRecalibrator
构建重新校准模型以评估变体质量以进行过滤
（[VariantRecalibrator](https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_vqsr_VariantRecalibrator.php)）
  2. 变异质量得分重新校准[ApplyVQSR](https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_vqsr_ApplyVQSR.php)
---
VariantRecalibrator

gatk VariantRecalibrator \ -R Homo_sapiens_assembly38.fasta \ -V input.vcf.gz \ --resource hapmap，known = false，training = true，truth = true，prior = 15.0：hapmap_3.3.hg38.sites.vcf.gz \ --resource omni，known = false，training = true，truth = false，prior = 12.0：1000G_omni2.5.hg38.sites.vcf.gz \ --resource 1000G，known = false，training = true，truth = false，prior = 10.0：1000G_phase1.snps.high_confidence.hg38.vcf.gz \ --resource dbsnp，known = true，training = false，truth = false，prior = 2.0：Homo_sapiens_assembly38.dbsnp138.vcf.gz \ -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \ -mode SNP -O output.recal \ --tranches-file output.tranches \ --rscript-file output.plots.R

VQSR

gatk ApplyVQSR \ -R Homo_sapiens_assembly38.fasta \ -V input.vcf.gz \ -O output.vcf.gz \ --truth-sensitivity-filter-level 99.0 \ --tranches-file output.tranches \ --recal-file output.recal \ -mode SNP

java -Xmx128g jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyVQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V 2447-20.vcf -O 2447-20.vqsr.vcf & java -Xmx128g jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar ApplyVQSR -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V 119-8.vcf -O 119-8.vqsr.vcf &

5. 转换vcf为table格式

java -Xmx128g -jar /home/chaim/disk/gatk/gatk4/gatk-package-4.0.10.1-local.jar VariantsToTable -R /home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa -V common_filtration.vcf -F CHROM -F POS -F REF -F ALT -GF AD -GF DP -GF GQ -GF PL -O common.table

三、[使用qtl-seqr包](https://www.jianshu.com/p/6ff6c3a44fb6)（可以和第四步同时进行）
[QTLseqr在github](https://github.com/bmansfeld/QTLseqr)
#*以下代码为R语言代码*

安装使用QTL-seqr

安装

installed.packages('devtools') library('devtools') install_github("bmansfeld/QTLseqr")

load the package

library("QTLseqr")

set sample and file name

HighBulk <- "样品名称1" LowBulk <- "样品名称2" file <- "common.table"

choose which chromosomes will be included in the analysis

chroms <- paste0 Chroms <- paste0(rep("Chr", 10), 1:10)

Import SNP data from file

df <- importFromGATK( file = file, highBulk = HighBulk, lowBulk = LowBulk, chromList = Chroms )

Filter SNPs based on some criteria

df_filt <- filterSNPs( SNPset = df, refAlleleFreq = 0.20, minTotalDepth = 100, maxTotalDepth = 400, minSampleDepth = 40, minGQ = 99 )

Run G' analysis

df_filt <- runGprimeAnalysis( SNPset = df_filt, windowSize = 1e6, outlierFilter = "deltaSNP")

Run QTLseq analysis

df_filt <- runQTLseqAnalysis( SNPset = df_filt, windowSize = 1e6, popStruc = "F2",#根据群体结构调整F2或者RIL bulkSize = c(25, 25),#根据池中的样本单株数量 replications = 10000, intervals = c(95, 99) )

Plot

plotQTLStats(SNPset = df_filt, var = "Gprime", plotThreshold = TRUE, q = 0.01) plotQTLStats(SNPset = df_filt, var = "deltaSNP", plotIntervals = TRUE)

export summary CSV

getQTLTable(SNPset = df_filt, alpha = 0.01, export = TRUE, fileName = "my_BSA_QTL.csv")

四. 使用snpEff处理gatk4输出的vcf。（此步骤和第三步骤的输出结果可以互相对比，两者具有相同的功能）
[snpEff的教程](https://www.jianshu.com/p/967006e0089c)
1.  [下载snpEff地址](https://sourceforge.net/projects/snpeff/files/snpEff_latest_core.zip/download)
    解压`unzip snpEff_latest_core.zip`
    我的路径是/home/chaim/bsa/snpEff/
2.  配置玉米zm437版本的数据库
    在路径`/home/chaim/bsa/snpEff/snpEff/`目录下创建文件夹`data`,

cd /home/chaim/bsa/snpEff/snpEff/ mkdir data cd data mkdir genomes mkdir zm437

在zm437目录存放基因组注释文件genes.gff, 蛋白库，protein.fa

在genomes目录放置基因组参考序列 zm437.fa


注意上述的基因组注释文件是gff3格式。
修改snpEff.config的参数
添加如下内容

maize genome,version zm437

zm437.genome:maize


回到snpEFF目录，运行命令
`java -jar snpEff.jar build -gff3 -v zm437`

3.  对vcf格式文件进行注释：
    `bwa`目录存放着GATK4处理之后的文件`common_filtration.vcf`
    在`/home/chaim/bsa/bwa`目录执行下面命令

java -Xmx128g -jar /home/chaim/bsa/snpEff/snpEff/snpEff.jar zm437 common_filtration.vcf > common.eff.vcf

```

会输出三个文件， snpEff_genes.txt snpEff_summary.html common.eff.vcf

如果想更改使用其他注释文件，删除/home/chaim/bsa/snpEff/snpEff//data/zm437/snpEffectPredictor.bin该文件删除即可。重新从步骤1开始即可。