GATK HaplotypeCaller只能处理单样本文件，但是我使用的bam文件也是单样本，为确定样本信息，通过samtools查看bam文件。bam文件标准格式如下，引用自 bam文件的理解 - (jianshu.com)
没有返回值，查看bam文件信息：
发现bam文件中没有RG信息，因此通过picard修改BAM文件的Read Group。
建议使用conda进行安装，省了很多麻烦。
用Picard的 AddOrReplaceReadGroups命令来完成：
--VALIDATION_STRINGENCY LENIENT这里是为了避免Picard因为检查Header 不符合规范而报错终止。
-Xms1g 是 Java Virtual Machine的参数，设置最小内存为1G。
-Xmx10g 也是 Java Virtual Machine的参数，设置最大内存为10G。
-XX:ParallelGCThreads=4 还是 Java Virtual Machine的参数，设置并行 线程为4。
再次查看RG信息：
参考材料：
https://wolfsonliu.github.io/archive/2020/gatk4-germline-variant-calling-fen-xi-liu-cheng.html
https://gatk.broadinstitute.org/hc/en-us/articles/4404607557787-AddOrReplaceReadGroups-Picard-
samtools是一个用于 *** 作sam和bam文件（通常是短序列比对工具如bwa,bowtie2,hisat2,tophat2等等产生的，具体格式可以在消息框输入“SAM”查看）的工具合集，包含有许多命令。以下是常用命令的介绍。
1.View
view命令的主要功能是：将sam文件与bam文件互换；然后对bam文件进行各种 *** 作，比如数据的排序(sort)和提取(这些 *** 作 是对bam文件进行的，因而当输入为sam文件的时候，不能进行该 *** 作)；最后将排序或提取得到的数据输出为bam或sam（默认的）格式。
bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。
view命令中，对sam文件头部（序列ID）的输入(-t或-T）和输出(-h)是单独的一些参数来控制的。
Usage: samtools view [options] <in.bam>|<in.sam>[region1 [...]]
默认情况下不加 region，则是输出所有的 region.
options:
-b output BAM
默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式
-h print header for the SAM output
默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息
-H print header only (no alignments)
仅仅输出文件的头
-S input is SAM
默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
-u uncompressed BAM output (force -b)
该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。
-c Instead of printing the alignments, only count them and print the
total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ , are taken into account.
过滤功统计功能
-c print only the count of matching records
-L FILE output alignments overlapping the input BED FILE [null]
-t FILE list of reference names and lengths (force -S) [null]
使用一个list文件来作为header的输入
-T FILE reference sequence file (force -S) [null]
使用序列fasta文件作为header的输入
-o FILE output file name [stdout]
-F INT filtering flag, 0 for unset [0]
Skip alignments with bits present in INT [0]
数字4代表该序列没有比对到参考序列上
数字8代表该序列的mate序列没有比对到参考序列上
过滤功能。如F12过滤只有双端map的
-q INT minimum mapping quality [0]
比对的最低质量值，一般认为20就为unique比对了，可以结合上述-bF参数使用使用提取特定的比对结果
例子：
将sam文件转换成bam文件
samtools view -bS abc.sam >abc.bam
BAM转换为SAM
samtools view -h -o out.sam out.bam
提取比对到参考序列上的比对结果
samtools view -bF 4 abc.bam >abc.F.bam
提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可
samtools view -bF 12 abc.bam >abc.F12.bam
提取没有比对到参考序列上的比对结果
samtools view -bf 4 abc.bam >abc.f.bam
提取bam文件中比对到caffold1上的比对结果，并保存到sam文件格式
samtools view abc.bam scaffold1 >scaffold1.sam
提取scaffold1上能比对到30k到100k区域的比对结果
samtools view abc.bam scaffold1:30000-100000 $gtscaffold1_30k-100k.sam
根据fasta文件，将 header 加入到 sam 或 bam 文件中
samtools view -T genome.fasta -h scaffold1.sam >scaffold1.h.sam
2.Sort
sort对bam文件进行排序。一些软件需要sort的bam或者sam文件，如stringtie，所以必须要sort使用；求depth时，也必须要sort；
Usage: samtools sort [-n] [-m <maxMem>] <in.bam><out.prefix>
-m 内存参数默认下是 500,000,000 即500M（不支持K，M，G等缩写）。对于处理大数据时，如果内存够用，则设置大点的值，以节约时间。
-n 设定排序方式按short reads的ID排序。默认下是按序列在fasta文件中的顺序（即header）和序列从左往右的位点排序。
例子：
samtools sort accept.bam accept.sort最终产生accept.sort.bam
3.merge
将2个或2个以上的已经sort了的bam文件融合成一个bam文件。融合后的文件已经sort过了的。
Usage: samtools merge [-nr] [-h inh.sam] <out.bam><in1.bam><in2.bam>[...]
Options: -n sort by read names
-r attach RG tag (inferred from file names)
-u uncompressed BAM output
-f overwrite the output BAM if exist
-1 compress level 1
-R STR merge file in the specified region STR [all]
-h FILE copy the header in FILE to <out.bam>[in1.bam]
例子：
4.index
必须对bam文件进行默认情况下的排序后，才能进行index。否则会报错。
建立索引后将产生后缀为.bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。比如samtool的tview命令就需要；gbrowse2显示reads的比对图形的时候也需要。IGV显示比对情况也需要。
Usage: samtools index <in.bam>[out.index]
例子：
以下两种命令结果一样
$ samtools index abc.sort.bam
$ samtools index abc.sort.bam abc.sort.bam.bai
5.faidx
对fasta文件建立索引,生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条（子）序列
Usage: samtools faidx <in.bam>[ [...]]
对基因组文件建立索引，方便提取序列。
例子：$ samtools faidx genome.fasta
由于有索引文件，可以使用以下命令很快从基因组中提取到fasta格式的子序列
$ samtools faidx genome.fasta scffold_10 >scaffold_10.fasta
6.tview
tview能直观的显示出reads比对基因组的情况，和基因组浏览器有点类似。
需要事先利用利用上面讲的sort和建index命令执行完后，用下述命令。
Usage: samtools tview <aln.bam>[ref.fasta]
出参考基因组的时候，会在第一排显示参考基因组的序列，否则，第一排全用N表示。
按下 g ，则提示输入要到达基因组的某一个位点。例子“scaffold_10:1000"表示到达第
10号scaffold的第1000个碱基位点处。
使用H(左）J（上）K（下）L（右）移动显示界面。大写字母移动快，小写字母移动慢。
使用空格建向左快速移动（和 L 类似），使用Backspace键向左快速移动（和 H 类似）。
Ctrl+H 向左移动1kb碱基距离； Ctrl+L 向右移动1kb碱基距离
可以用颜色标注比对质量，碱基质量，核苷酸等。30～40的碱基质量或比对质量使用白色表示；
20～30黄色；10～20绿色；0～10蓝色。
使用点号'.'切换显示碱基和点号；使用r切换显示read name等
还有很多其它的使用说明，具体按 ？ 键来查看。
7.flagstat
给出BAM文件的比对结果
Usage: samtools flagstat <in.bam>
$ samtools flagstat example.bam
11945742 + 0 in total (QC-passed reads + QC-failed reads)
#总共的reads数
0 + 0 duplicates
7536364 + 0 mapped (63.09%:-nan%)
#总体上reads的匹配率
11945742 + 0 paired in sequencing
#有多少reads是属于paired reads
5972871 + 0 read1
#reads1中的reads数
5972871 + 0 read2
#reads2中的reads数
6412042 + 0 properly paired (53.68%:-nan%)
#完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
6899708 + 0 with itself and mate mapped
#paired reads中两条都比对到参考序列上的reads数
636656 + 0 singletons (5.33%:-nan%)
#单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。
469868 + 0 with mate mapped to a different chr
#paired reads中两条分别比对到两条不同的参考序列的reads数
243047 + 0 with mate mapped to a different chr (mapQ>=5)
#同上一个，只是其中比对质量>=5的reads的数量
8.depth
得到每个碱基位点的测序深度,并输出到标准输出，所以要用大于号追加到一个文件。
Usage: bam2depth [-r reg] [-q baseQthres] [-Q mapQthres] [-b in.bed] <in1.bam>[...]
-r 后面跟染色体号（region）
-q ：计算深度时要求测序碱基质量最低质量值
-Q ：计算深度时要求比对的最低质量值
注意：做depth之前必须做samtools index；
例子
samtools depth accept.bam >depth
9.其他命令
reheader：替换bam文件的头
$ samtools reheader <in.header.sam><in.bam>
idxstats ：统计一个表格，4列，分别为”序列名，序列长度，比对上的reads数，unmapped reads number。第4列应该是paired reads中有一端能匹配到该scaffold上，而另外一端不匹配到任何scaffolds上的reads数。
$ samtools idxstats <aln.bam>
rmdup：将由PCR duplicates获得的reads去掉，并只保留最高比对质量的read。
Usage: samtools rmdup [-sS]
-s 对single-end reads。默认情况下，只对paired-end reads
-S 将Paired-end reads作为single-end reads处理。
10. 将bam文件转换为fastq文件
有时候，我们需要提取出比对到一段参考序列的reads，进行小范围的分析，以利于debug等。这时需要将bam或sam文件转换为fastq格式。
该网站提供了一个bam转换为fastq的程序：http://www.hudsonalpha.org/gsl/information/software/bam2fastq
$ wget http://www.hudsonalpha.org/gsl/static/software/bam2fastq-1.1.0.tgz
$ tar zxf bam2fastq-1.1.0.tgz
$ cd bam2fastq-1.1.0
$ make
$ ./bam2fastq <in.bam>
11. mpileup
samtools还有个非常重要的命令mpileup，以前为pileup。该命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件，在samtools的安装文件夹中可以找到。
最常用的参数有2：
-f 来输入有索引文件的fasta参考序列； -g 输出到bcf格式。用法和最简单的例子如下
Usage: samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]$ samtools mpileup -f genome.fasta abc.bam >abc.txt
$ samtools mpileup -gSDf genome.fasta abc.bam >abc.bcf
$ samtools mpileup -guSDf genome.fasta abc.bam | \
bcftools view -cvNg - >abc.vcf
mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。比如：
scaffold_1 2841 A 11 ,,,...,.... BHIGDGIJ?FF
scaffold_1 2842 C 12 ,$,,...,....^I. CFGEGEGGCFF+
scaffold_1 2843 G 11 ,,...,..... FDDDDCD?DD+
scaffold_1 2844 G 11 ,,...,..... FA?AAAA<AA+
scaffold_1 2845 G 11 ,,...,..... F656666166*
scaffold_1 2846 A 11 ,,...,..... (1.1111)11*
scaffold_1 2847 A 11 ,,+9acggtgaag.+9ACGGTGAAT.+9ACGGTGAAG.+9ACGGTGAAG,+9acggtgaag.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG.+9ACGGTGAAG %.+....-..)
scaffold_1 2848 N 11 agGGGgGGGGG !!$!!!!!!!!
scaffold_1 2849 A 11 c$,...,..... !0000000000
scaffold_1 2850 A 10 ,...,..... 353333333
mpileup生成的结果包含6行：参考序列名；位置；参考碱基；比对上的reads数；比对情况；比对上的碱基的质量。其中第5列比较复杂,解释如下：
1 ‘.’代表与参考序列正链匹配。
2 ‘,’代表与参考序列负链匹配。
3 ‘ATCGN’代表在正链上的不匹配。
4 ‘atcgn’代表在负链上的不匹配。
5 ‘*’代表模糊碱基
6 ‘^’代表匹配的碱基是一个read的开始；’^'后面紧跟的ascii码减去33代表比对质量；这两个符号修饰的是后面的碱基，其后紧跟的碱基(.,ATCGatcgNn)代表该read的第一个碱基。
7 ‘$’代表一个read的结束，该符号修饰的是其前面的碱基。
8 正则式’\+[0-9]+[ACGTNacgtn]+’代表在该位点后插入的碱基；比如上例中在scaffold_1的2847后插入了9个长度的碱基acggtgaag。表明此处极可能是indel。
9 正则式’-[0-9]+[ACGTNacgtn]+’代表在该位点后缺失的碱基；
12. 使用bcftools
bcftools和samtools类似，用于处理vcf(variant call format)文件和bcf(binary call format)文件。前者为文本文件，后者为其二进制文件。
bcftools使用简单，最主要的命令是view命令，其次还有index和cat等命令。index和cat命令和samtools中类似。此处主讲使用view命令来进行SNP和Indel calling。该命令的使用方法和例子为：
$ bcftools view [-AbFGNQSucgv] [-D seqDict] [-l listLoci] [-s listSample]
[-i gapSNPratio] [-t mutRate] [-p varThres] [-P prior]
[-1 nGroup1] [-d minFrac] [-U nPerm] [-X permThres]
[-T trioType] in.bcf [region]
$ bcftools view -cvNg abc.bcf >snp_indel.vcf
生成的结果文件为vcf格式，有10列，分别是：1 参考序列名；2 varianti所在的left-most位置；3 variant的ID（默认未设置，用’.'表示)；4 参考序列的allele；5 variant的allele(有多个alleles，则用’,'分隔)6 variant/reference QUALity7 FILTers applied8 variant的信息，使用分号隔开；9 FORMAT of the genotype fields, separated by colon (optional)； 10 SAMPLE genotypes and per-sample information (optional)。
例如：
scaffold_1 2847 . A AACGGTGAAG 194 . INDELDP=11VDB=0.0401AF1=1AC1=2DP4=0,0,8,3MQ=35FQ=-67.5 GT:PL:GQ 1/1:235,33,0:63
scaffold_1 3908 . G A 111 . DP=13VDB=0.0085AF1=1AC1=2DP4=0,0,5,7MQ=42FQ=-63 GT:PL:GQ 1/1:144,36,0:69
scaffold_1 4500 . A G 31.5 . DP=8VDB=0.0034AF1=1AC1=2DP4=0,0,1,3MQ=42FQ=-39 GT:PL:GQ 1/1:64,12,0:21
scaffold_1 4581 . TGGNGG TGG 145 . INDELDP=8VDB=0.0308AF1=1AC1=2DP4=0,0,0,8MQ=42FQ=-58.5 GT:PL:GQ 1/1:186,24,0:45
scaffold_1 4644 . G A 195 . DP=21VDB=0.0198AF1=1AC1=2DP4=0,0,10,10MQ=42FQ=-87 GT:PL:GQ 1/1:228,60,0:99
scaffold_1 4827 . NACAAAGA NA 4.42 . INDELDP=1AF1=1AC1=2DP4=0,0,1,0MQ=40FQ=-37.5 GT:PL:GQ 0/1:40,3,0:3
scaffold_1 4854 . A G 48 . DP=6VDB=0.0085AF1=1AC1=2DP4=0,0,2,1MQ=41FQ=-36 GT:PL:GQ 1/1:80,9,0:16
scaffold_1 5120 . A G 85 . DP=8VDB=0.0355AF1=1AC1=2DP4=0,0,5,3MQ=42FQ=-51 GT:PL:GQ 1/1:118,24,0:45
第8列中显示了对variants的信息描述，比较重要，其中的 Tag 的描述如下：
Tag Format Description
AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
DP int Raw read depth (without quality filtering)
DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
FQ int Consensus quality. Positive: sample genotypes differentnegative: otherwise
MQ int Root-Mean-Square mapping quality of covering reads
PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
QCHI2 int Phred-scaled PCHI2
RP int # permutations yielding a smaller PCHI2
CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
UGT string Most probable genotype configuration without the trio constraint
CGT string Most probable configuration with the trio constraint
使用bcftools得到variant calling结果后。需要对结果再次进行过滤。主要依据比对结果中第8列信息。其中的 DP4 一行尤为重要，提供了4个数据：1 比对结果和正链一致的reads数、2 比对结果和负链一致的reads数、3 比对结果在正链的variant上的reads数、4 比对结果在负链的variant上的reads数。可以设定 （value3 + value4）大于某一阈值，才算是variant。比如：
$ perl -ne 'print $_ if /DP4=(\d+),(\d+),(\d+),(\d+)/ &&($3+$4)>=10 &&($3+$4)/($1+$2+$3+$4)>=0.8' snp_indel.vcf >snp_indel.final.vcf
目录
samtools 和 picard 都有对SAM/BAM文件进行排序的功能，一般都是基于坐标排序（还提供了 -n 选项来设定用reads名进行排序），先是对chromosome/contig进行排序，再在chromosome/contig内部基于start site从小到大排序，对start site排序很好理解，可是对chromosome/contig排序的时候是基于什么标准呢？
基于你提供的 ref.fa 文件中的chromosome/contig的顺序 。当你使用比对工具将fastq文件中的reads比对上参考基因组后会生成SAM文件，SAM文件包含头信息，其中有以 @SQ 开头的头信息记录，reference中有多少条chromosome/contig就会有多少条这样的记录，而且它们的顺序与 ref.fa 是一致的。
当使用samtools或picard对SAM/BAM文件进行排序时，这些工具就会读取头信息，按照头信息指定的顺序来排chromosome/contig。所以进行排序时需要提供包含头信息的SAM/BAM文件。
那么 普通情况下我们的chromosome/contig排序情况是什么样的?
一般情况下在进行SAM文件的排序时，染色体的排序到底是按照哪种规则进行排序的，不是一个很重要的问题，也不会对后续的分析产生影响，但是在执行GATK流程时，GATK对染色体的排序是有要求的，必须按照从chr1开始到chr22，最后是chrX和chrY这样的顺序，否则会报错
面对这样变态的要求，我们怎么解决？
在构造ref.fa文件时，让它按照从chr1开始到chr22，最后是chrX和chrY这样的顺序进行组织就可以了：
FLAG列在SAM文件的第二列，这是一个很重要的列，包含了很多mapping过程中的有用信息，但很多初学者在学习SAM文件格式的介绍时，遇到FLAG列的说明，常常会一头雾水
what?还二进制，这也太反人类的设计了吧！
不过如果你站在开发者的角度去思考这个问题，就会豁然开朗
在mapping过程中，我们想记录一条read的mapping的信息包括：
这些信息总结起来总共包括以下12项：
而每一项又只有两种情况，是或否，那么我可以用一个12位的二进制数来记录所有的信息，每一位表示某一项的情况，这就是原始FLAG信息的由来，但是二进制数适合给计算机看，不适合人看，需要转换成对应的十进制数，也就有了我们在SAM文件中看到的FLAG值
但是FLAG值所包含信息的解读还是要转换为12位的二进制数
SAM格式文件的第3和第7列，可以用来判断某条reads是否比对成功到了基因组的染色体，左右两条reads是否比对到同一条染色体
有两个方法可以提取未比对成功的测序数据：
对于PE数据，在未比对成功的测序数据可以分成3类：
看完这一部分，是不是有一个感觉： FLAG玩得溜，SAM文件可以处理得出神入化
首先，思考一个问题： 对于PE数据，一条测序片段（fragment）有read1和read2两条测序片段，它们俩的名字相同，那么对于这一条测序片段，对它进行mapping之后得到的SAM文件中会出现几条记录呢？
对于我的这个假设可以用以下的方法进行验证：
上面的测试结果与我们的假设吻合
但是在一次处理三代测序数据（三代测序数据是Single-End）中发现了不同：
在输出中出现了一些不太和谐的结果：有极少部分的QNAME对应2条以上的记录，这意味着存在一条read会有多条比对记录的情况，why？
对这个与预期不完全相符的结果，尝试去寻找里面的原因，其间进行了各种各样的推理、假设、验证，最终在 李恒的github 中找到了答案
这种情况容易在三代测序数据中出现
如果你用的是Single-End的数据，那么差异应该比较小，不会太明显，而在Pair-End上差异可能会比较大，之所以会产生这些差异，原因有两点：
从上面列出的两点差异可以看出，mpileup默认输出的是高质量的覆盖深度，这是有历史原因的：当场mpileup功能被开发出来就是为了与bcftools组合，将samtools mpileup的输出作为bcftools的输入用于下游的snp-calling，当然需要保证数据的质量
当然可以通过设置对应的参数使得它的属于结果与depth的一致，但是不推荐这么做
下面是对同一个样本的paired-end Fastaq文件比对结果（比对使用hisat2），hisat2和samtools分别给出的mapping rate的统计
hisat2：
samtools flagstat：
从上面可以看出，hisat2给出的mapping rate为85.40%，而samtools给出的为86.32%，两个的统计结果不一样，而且samtools的统计会大一些，what？
介四什么鬼？
我们来简单地分析一下：
hisat2中，
samtools中，
计算没问题，那问题出在哪呢？
有没有注意到上面的两个式子中的分子和分母，计算它们的差值：
分子：
分母：
发现了没有，它们的差值正好都等于samtools flagstat的输出结果的第二行：
所以，hisat2和samtools计算mapping rate的公式实际上分别为：
一般来说，我们想得到的是hisat2计算公式所得到的统计结果，hisat2统计结果在比对结束后会以标准错误形式给出，我们可以将标准错误重定向到一个log文件中，但是如果我们忘了保持这个统计结果，怎么办？
最简单的办法就是重新跑一遍hisat2，但是这样太耗费时间和计算资源了，这时我们可以利用samtools flagstat对SAM文件的统计结果，以及它的部分统计值与hisat2计算公式的关系，快速地算出准确的mapping rate：
参考资料：
(1) 【】从零开始完整学习全基因组测序数据分析：第5节 理解并 *** 作BAM文件
(2) 【生信技能树】【直播】我的基因组（十五）:提取未比对的测序数据
(3) BWA's README in github
(4) 【】黄树嘉《样本量重要，还是测序深度重要? 生物信息工程师可以分为多少种类型? |《解螺旋技术交流圈》精华第3期》
(5) 生信媛《HISAT2的比对率计算结果和SAMTools flagstat不同，你想过原因吗？ 》
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）