小马的生信笔记

重测序数据 Call SNP

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
				
			

4.2 INDEL

发表评论