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

1、测序质量评估

## conda安装fastqc
fastqc *.fq.gz

 2、合并评估报告

## conda安装multiqc (可能会要求安装latex,为生成PDF文件所需)
multiqc .

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

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

4、下载参考基因组(以玉米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

5、构建索引

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

##构建基因组索引
samtools faidx genome.fasta

##构建字典序文件
picard CreateSequenceDictionary -R genome.fasta

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

## ID:S027为样品ID，SM:S027为样品名称，PL:illumina为测序平台
bwa mem -t 10 -R '@RGtID:S027tSM:S027tPL:illumina' genome.fasta S027.R1.fq.gz  S027.R2.fq.gz > test.sam

##将sam文件转换成bam文件
samtools view -bS test.sam > test.bam

7、查看及统计比对结果

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

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

8、去除重复

##对bam文件排序
samtools sort -@ 4 -o test.sorted.bam test.bam

##去除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 

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