重測序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等（最近流行的分析方法）