1. 前言
随着基因组测序价格的降低,近几点发表了越来越多的动植物基因组,因此基于基因组重测序进行群体遗传学的研究也逐渐火热。再加上后续的研究计划也会涉及群体遗传学,所以来提前学习一下全部的分析流程。今天就从Call SNP开始学习。
2. 软件下载
console.log( 'Code is Poetry' );
3. 开始分析
首先第一步就是要准备分析需要的数据,也就是(1)参考基因组(2)重测序数据(要是Clean data,raw data去跑一下fastp就可以)
3.1 为参考基因组构建索引
samtools faidx genome.fasta
bwa index genome.fasta
picard.jar CreateSequenceDictionary R=genome.fasta
3.2 将reads比对到参考基因组
bwa mem -t 2 -R '@RG\tID:S1\tSM:S1\tPL:illumina' Genome.fasta S1_1.fq.gz S1_2.fq.gz | samtools sort -@ 2 -o S1.sort.bam -
java -Xmx4g -XX:ParallelGCThreads=2 -jar picard.jar MarkDuplicates I=S1.sort.bam O=S1.sort.markdup.bam CREATE_INDEX=true REMOVE_DUPLICATES=true M=S1.marked_dup_metrics.txt
3.3 GATK 检测变异
GATK检测变异的时一般是针对于单个染色体完成,所以要分染色体进行分析。
gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" HaplotypeCaller -R genome.fasta -I S1.sort.markdup.bam -L chr1 -ERC GVCF -O S1.chr1.g.vcf.gz
接下来有两种路线可以选,(1)合并成1个gvcf文件然后进行变异检测(2)分样本进行变异检测,本研究中只对多样本进行研究,并且样本数量不超过1000,所以选择(1),但是(2)也会写下教程,但是不会验证脚本是否有错误
3.3.1 (1)
gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" CombineGVCFs -R genome.fasta -V gvcf.chr1.list -O chr1.g.vcf.gz
# -R 参考基因组
# -V 样品list
# -O 输出结果
cat gvcf.chr1.list
S1/S1.chr1.g.vcf.gz
S2/S2.chr1.g.vcf.gz
检测变异
gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" GenotypeGVCFs -R genome.fasta -V chr1.g.vcf.gz -O chr1.raw.vcf.gz
# -R 参考基因组
# -V 输入vcf文件
#-O 输出结果文件
合并不同染色体的VCF文件,形成全基因组的VCF文件
gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" MergeVcfs -I raw_vcf.list -O all.merge_raw.vcf.gz
cat raw_vcf.list
chr1.raw.vcf.gz
chr2.raw.vcf.gz
4. 提取SNP和INDEL并且进行过滤
4.1 SNP
gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" SelectVariants -R genome.fasta -V all.merge_raw.vcf --select-type SNP -O all.raw.snp.vcf
gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" 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 'SNP_filter' -O all.filter.snp.vcf
gatk --java-options "-Xmx4g -Djava.io.tmpdir=./tmp" SelectVariants -R genome.fasta -V all.filter.snp.vcf --exclude-filtered -O all.filtered.snp.vcf