GWAS数据分析流程—SNP、Indel提取

1、测序质量评估

## conda安装fastqc
fastqc *.fastq.gz -t 8 -o fastqc_out/

可选（合并质控报告）
## conda安装multiqc (可能会要求安装latex,为生成PDF文件所需)
multiqc .

2、数据过滤(重复、接头等)

## conda安装fastp（fastp也会生成过滤前后的质量报告，如果需要的话，需要对报告文件重命名）
fastp -i in.R1.fq.gz -I in.R2.fq.gz -o out.R1.fq.gz -O out.R2.fq.gz

3、下载参考基因组(以玉米v5参考基因组为例)

## 下载参考基因组数据
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/902/167/145/GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0/GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0_genomic.fna.gz

##解压基因组文件
gunzip GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0_genomic.fna.gz

##更改名字
mv GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0_genomic.fna genome.fasta

4、构建索引

## 构建BWA索引
bwa index genome.fasta

bwa index的用法和参数：
bwa index [-p prefix] [-a algoType] <in.db.fasta>
-p  #输出数据库的前缀，默认与输入文件名一致
-a  #设置参考基因组构建index的算法（is和bwtsw两种）

bwa构建参考基因组index的两种算法：
-is 用于构建后缀数组的线性时间算法：缺点:比较耗内存（需要参考基因组5.37倍的内存）；优点:速度快。
 is是构建index的默认算法，但是不适合 “参考基因组” 大于2GB的文件。
-bwtsw 使用BWT-SW算法，多用于人类基因组。

##构建基因组索引，该命令会在运行目录下创建一个genome.fai索引文件
samtools faidx genome.fasta

##构建字典序文件，生成dict文件
gatk CreateSequenceDictionary -R example.fasta -O example.dict  或
picard CreateSequenceDictionary -R genome.fasta

GATK (全称The Genome Analysis Toolkit)是Broad Institute开发的二代重测序数据分析软件，是基因分析的工具集。在4.0以后，GATK包含有Picard工具集，所有Picard工具都能够使用GATK完成。

5、将过滤后的reads比对到参考基因组上，并将生成的sam文件转换成bam文件

## ID:S027为样品ID，SM:S027为样品名称，PL:illumina为测序平台，输入文件为参考基因组和质控过滤后的文件
bwa mem -t 10 -R '@RG\tID:S027\tSM:S027\tPL:illumina' genome.fasta S027.R1.fq.gz S027.R2.fq.gz | samtools sort -@ 4 -m 1G -o s027.sort.bam -
# bwa比对命令最好不要使用nohup，nohup的输出信息会被bwa输出到目标文件，可能会影响后续步骤。
# mem算法 -t 运行线程数 -R添加头部信息，-R一定不能省略或写错。
# bwa的mem的效率更高，且更加准确。 

6、查看及统计比对结果

##查看bam文件
samtools view -h test.bam | less -S

##统计bam文件信息
samtools flagstat test.bam > test.bam.flagstat 

7、去除重复

##去除bam文件中的重复
picard -Xmx4g MarkDuplicates I=test.sorted.bam O=test.sorted.rmdup.bam CREATE_INDEX=true REMOVE_DUPLICATES=true M=test.marked_dup_metrics.txt 

# I 输入排序后的文件; O 输出文件; M 输出重复矩阵。 

9、检测变异位点

##检测变异
nohup gatk --java-options "-Xmx10g" HaplotypeCaller -R genome.fasta -I test.sorted.rmdup.bam -ERC GVCF -O test.g.vcf &

10、合并变异位点

##生成合并列表
ls *.g.vcf > gvcf.list 

##合并GVCF文件
nohup gatk --java-options "-Xmx10g"  CombineGVCFs -R genome.fasta -V gvcf.list -O all.g.vcf &

11、提取所有样本变异位点

nohup gatk --java-options "-Xmx4g" GenotypeGVCFs -R genome.fasta -V all.g.vcf -O all.var.vcf &

12.1、分选SNP

## 从变异里提取所有SNP
gatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.var.vcf --select-type SNP -O all.raw.snp.vcf

## 过滤不符合条件的SNP
gatk  --java-options "-Xmx4g" VariantFiltration -R genome.fasta -V all.raw.snp.vcf  --filter-expression "QD < 2.0 || MQ<40.0 || FS > 60.0 ||  SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"  --filter-name 'LowQualFilter' -O all.snp.fil.vcf  

## 提取过滤后的SNP
gatk  --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.snp.fil.vcf  --exclude-filtered -O all.snp.filed.vcf

12.2、分选Indel

## 从变异里提取所有Indel
gatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.var.vcf  --select-type INDEL -O all.raw.indel.vcf

## 过滤不符合条件的Indel
gatk  --java-options "-Xmx4g" VariantFiltration -R genome.fasta -V all.raw.indel.vcf  --filter-expression "QD < 2.0 ||  FS > 200.0 ||  SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0"  --filter-name 'LowQualFilter' -O all.indel.fil.vcf  

##提取过滤后的Indel
gatk  --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.indel.fil.vcf   --exclude-filtered -O all.indel.filed.vcf