【变异检测】SNP INDEL检测
重测序SNP和INDEL变异检测主要有3套方案，samtools+bcftools、GATK、assembly-versus-assembly
1.samtools+bcftools方法
samtools和bcftools是我们做生信最常用的软件之一，功能非常的强大。这套软件组合可以单call，也可以群call来检测样本变异。单call提供一个样本，群call就提供多个样本。samtools+bcftools组合效率高，占用计算资源不大，但准确性略低于GATK，但得到变异数量往往多GATK。
BWA分析可以看：用BWA进行序列比对 - 知乎 (zhihu.com)
输入文件是bam文件（BWA比对且去重后的bam结果，基于有参比对），格式如下
如果要看bam里面的信息可以用如下命令：
samtools view -h LS10.rmdup.bam|less
由于效率原因，群call时候可以通过切割基因组或按染色体块形式分析来加快效率
# call variant
samtools mpileup -q 1 -C 50 -S -D -m 2 -F 0.002 -uf ref.fa -b bam.list -l split1.txt | bcftools call -mv -f GQ |filter varFilter filteropt - > split1.raw.vcf
awk '/^#/ || /\tINDEL;/' split1.raw.vcf| gzip -c >split1.INDEL.vcf.gz
awk '/^#/ || !/\tINDEL;/' split1.raw.vcf| gzip -c >split1.SNP.vcf.gz
其中split1.txt格式
bam.list格式
ref.fa格式
后续对切割的calling的变异vcf进行merge在一起就OK了
2.picard+GATK
目前GATK版本到了4版本，使用GATK更加方便。支持单call和群call模式，给一个样本为单call，同时给多个样本为群call
GATK优点，1.可以对于不同批次测序样本分别call得到中间文件gVCF（HaplotypeCaller），只要参考基因组相同，gvcf都可以随时合并进行后续分析，不需要再重新群call，后续可以只保存gvcf即可；2.准确性高，由于GATK call的结果，可以通过GATK参数进行低质量过滤，得到高质量SNP；3.功能全面且强大
缺失也很明显，计算资源消耗大，耗时长（java编写的软件），之前3版本总是会出现跑中断的情况，4版本已经解决。
输入文件同样是bam文件，可以提供bam文件list
变异检测前需要对基因组建立picard的索引（字典），需要安装下picard
java -jar picard.jar ref.fa O=ref.dict TMP_DIR=tmp
GATK使用需要多个模块
# HaplotypeCaller
gatk --java-options "-Xmx45G" HaplotypeCaller -R ref.fa -ERC GVCF -I MG_16B.rmdup.bam -O MG_16B.rmdup.bam.gvcf.gz --tmp-dir tmp --native-pair-hmm-threads 10
# CombineGVCFs 可以按区域进行合并，加快效率
gatk --java-options "-Xmx45G" CombineGVCFs -R ref.fa -L 1:1-50000000 -O rawVariantSplit.1:1-50000000.gvcf.gz --tmp-dir=tmp -V xx1.gvcf.gz -V xx2.gvcf.gz …
# GenotypeGVCFs
gatk --java-options "-Xmx45G" GenotypeGVCFs -R ref.fa -V rawVariantSplit.1:1-50000000.gvcf.gz -O GenoType.1:1-50000000.vcf.gz -G StandardAnnotation --tmp-dir=tmp
# MergeVcfs 提供区块genotype vcf文件路径的list
gatk --java-options "-Xmx45G" MergeVcfs -I GenotypeVCFs.list -O rawVariants.vcf.gz
# 切分snp和indel
gatk --java-options "-Xmx45G" SelectVariants -select-type SNP -V rawVariants.vcf.gz -O rawSnp.vcf.gz
gatk --java-options "-Xmx45G" SelectVariants -select-type INDEL -V rawVariants.vcf.gz -O rawIndel.vcf.gz
# VariantFiltration，GATK硬过滤，得到高质量变异位点
#snp
gatk --java-options "-Xmx45G" VariantFiltration -V rawSnp.vcf.gz --filter-expression "QD < 2.0" --filter-name QDFilter --filter-expression "MQ < 40.0" --filter-name MQFilter --filter-expression "FS > 60.0" --filter-name FSFilter --filter-expression "MQRankSum < -12.5" --filter-expression "ReadPosRankSum < -8.0" --filter-name "PosRankFilter" --filter-name "MQRankFilter" -O GATKfilter.snp.vcf.gz
#indel
gatk --java-options "-Xmx45G" VariantFiltration -V rawIndel.vcf.gz --filter-expression "QD < 2.0" --filter-name QDFilter --filter-expression "FS > 200.0" --filter-name FSFilter --filter-expression "ReadPosRankSum < -20.0" --filter-name "PosRankFilter" -O GATKfilter.Indel.vcf.gz
注：对于samtools call变异vcf和GATK call vcf文件后续需要过滤，过滤掉条件缺失率miss>0.1（标准，可以根据自己情况而定），次级等位基因maf<0.05（标准，可以根据自己情况而定），个体深度最低和最高条件过滤，基因型质量GQ<20，双等位基因过滤（如果后续分析不考虑双等位基因情况）等。
3.assembly-versus-assembly method（其实就是基因组之间比对检测snp）
基因组之间比对得到变异准确性理论是要高于重测比对（仅于单call）,但除了软件因素外，基因组的组装错误是导致比对的准确性因素之一。
lastz软件下载和使用教程：LASTZ (psu.edu)
#建索引
faSplit sequence ref.fa 50 ref
faToTwoBit ref.fa ref.fa.2bit
faSplit sequence query.fa 50 query
faToTwoBit ref.fa query.fa.2bit
lastdb -uNEAR -cR11 ref.split1.fa
faToTwoBit ref.split1.fa ref.split1.2bit
lastdb -uNEAR -cR11 query.split1.fa
faToTwoBit query.split1.fa query.split1.2bit
…(切割多个区域)
# chain（或者不需要分割区域，直接全长基因组之间比对得到genomeAlignment.axt）
lastz refsplit_out/ref.split1.2bit[multi] querysplit_out/query.split1.2bit[multi] M=254 K=4500 L=3000 Y=15000 E=150 H=2000 O=600 T=2 --format=axt > refsplit1.2bit_querysplit1.2bit.axt
axtChain -linearGap=medium refsplit1.2bit_querysplit1.2bit.axt ref.split1.2bit query.split1.2bit refsplit1.2bit_querysplit1.2bit.axt.chain
# merge+共线性
chainMergeSort *.chain >all.chain
chainPreNet all.chain ref.size query.size all.chain.filter
chainNet all.chain.filter ref.size query.size all.chain.filter.tnet all.chain.filter.qnet
netSyntenic
# 变异
axt_correction genomeAlignment.axt > all.chain.filter.tnet.synnet.axt.correct
axtSort all.chain.filter.tnet.synnet.axt.correct all.chain.filter.tnet.synnet.axt.correct.sort
axtBest all.chain.filter.tnet.synnet.axt.correct.sort all all.chain.filter.tnet.synnet.axt.correct.sort.best
lastz_postprocess all.chain.filter.tnet.synnet.axt.correct.sort.best 3 > all.chain.filter.tnet.synnet.axt.correct.sort.best.pp
trash_comments all.chain.filter.tnet.synnet.axt.correct.sort.best.pp > all.chain.filter.tnet.synnet.axt.correct.sort.best.pp.tn
fa_stat.pl zgmh_maskx.fa >genome.fa.length
trans_pos.pl genome.fa.length all.chain.filter.tnet.synnet.axt.correct.sort.best.pp.tn > all.chain.filter.tnet.synnet.axt.correct.sort.best.pp.tn.tp
msort -L4 -km5,n6,rn7 all.chain.filter.tnet.synnet.axt.correct.sort.best.pp.tn.tp > all.chain.filter.tnet.synnet.axt.correct.sort.best.pp.tn.tp.ms
best_hit all.chain.filter.tnet.synnet.axt.correct.sort.best.pp.tn.tp.ms > all.chain.filter.tnet.synnet.axt.correct.sort.best.pp.tn.tp.ms.bh
size格式，第一列为染色体编号，第二列为染色体长度
### fa_stat.pl
my $fa_file=shift;
die "perl $0 fa_file\n" unless $fa_file;
open IN,$fa_file or die $!;
$/=">";$/=<IN>;$/="\n";
while (<IN>){
chomp;
(my $id=$_)=~s/\s+.*$//;
$/=">";
my $seq=<IN>;
chomp $seq;
$seq=~s/\s+//g;
$/="\n";
my $length = length($seq);
$seq =~ s/N//g;
my $seq_len = length($seq);
print "$id\t$length\t$seq_len\n";
}
close IN;
### trans_pos.pl
use strict;
use warnings;
use Data::Dumper;
my (%info);
my ($info_file, $align_file) = @ARGV;
die "perl $0 scaffold_info_file alignment_file\n" unless $info_file && $align_file;
read_info($info_file, \%info);
open my $align_fh, $align_file or die $!;
while (<$align_fh>){
chomp;
my @tmp = split;
my $target_seq = <$align_fh>;
my $queue_seq = <$align_fh>;
my $carritage = <$align_fh>;
if ($tmp[7] eq "+"){
print join(" ", @tmp) , "\n";
print "$target_seq";
print "$queue_seq";
print "$carritage";
}else{
my $len = $info{$tmp[4]};
my ($start, $end) = ($len - $tmp[6] + 1, $len - $tmp[5] + 1);
print join(" ", @tmp[0..4]) . " $start $end $tmp[7] $tmp[8]\n";
print "$target_seq";
print "$queue_seq";
print "$carritage";
}
}
close $align_fh;
sub read_info{
my ($file, $ref) = @_;
open my $fh, "<", $file or die $!;
while (<$fh>){
chomp;
my @tmp = split;
$ref->{$tmp[0]} = $tmp[1];
}
close $fh;
}
注： 对于基因组call变异方法很多，下次详细讲基因组之间变异检测SNP，INDEL，PAV，CNV，SV等（最近流行的分析方法）
