67.《Bioinformatics Data Skills》之
pileup文件是指通过BAM文件每个位置重叠的read对比对结果进行的总结，可用于判断各个位点突变的可能性。这里通过命令行工具演示使用pileup文件进行突变识别（仅供参考，突变识别有更专业的工具与流程）。
数据下载地址：
NA12891_CEU_sample.bam
human_g1k_v37.fasta(ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz)
我们可以通过samtools的子命令mpileup来展示例子数据NA12891_CEU_sample.bam文件的pileup结果：
$ samtools mpileup --no-BAQ --fasta-ref human_g1k_v37.fasta --region 1:215906528-215906567 NA12891_CEU_sample.bam #1
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
1 215906528 G 21 ,,,,,,,,.,,,.,,..,.,, ;=?./:?>>;=7?>>@A?==: #2
...
1 215906545 G 21 ,.,,,.,,..,.,,,.,...^S. 88?<;5>=9?>=>>6:=>B2?
1 215906546 G 23 ,.,$,,.,,..,.,,,.,.C...^F, :7=727>=94;<==3:<;5C?<9
1 215906547 C 15 gGg$,GggGG,,.... <;80;><9=86=C>= #3
1 215906548 G 19 c$,ccC.,.,,,.,....,^]. ;58610=7=>75=7<463;
1 215906549 G 22 .,.$,$,..,.,,,.,.....,.^]. 13;==9;><=>76=79C?;8<6
1 215906550 G 21 .,,..,.,,,.,.....,..^:. 7:>6>><==;9=69C?=.;?;
1 215906551 G 18 .,,$.,.,,,.,....... 75;:>;==79=7C?E;:;
1 215906552 G 19 .,.$.,.,,,.,........ 9597><==69=<9D8E;99
1 215906553 G 16 .,$,.,,,.,....... 54:=<<69<<C?=;9;
1 215906554 C 15 G,,,,A,....,... 2@<<4/>4E854>9<
1 215906555 G 16 .$aaaaaA.AAAaAAA^:A 2@>?8?;<:335?:A> #4
...
解释：
#1.
mpileup参数--no-BAQ代表不考虑read比对质量（最后介绍），--fasta-ref接参考基因组，--region代表区域（不使用这个参数则展示整个基因组）。#2. pileup结果包括以下列：第一列为染色体编号；第二列为位置；第三列为参考碱基；第四列为深度；第五列为比对结果，其中”.“代表匹配而“,”代表反向匹配，大小写的字母分别代表正向链与负向链上的错配，“+”与"-"号代表插入与缺失，紧跟着插入与缺失的长度与序列，“^”与“$”代表read的开始与结束，"^"后的ASCII字符代表此read比对质量；最后一列为碱基测序质量。
#3. 此行同时存在大小写字母，代表正向或负向的错配都成立。
#4. 此行大部分的碱基都是错配的，“:”号代表其比对质量只有25。
mpileup命令加上-v参数可以获取每个位置变异的基因型与可能性，返回压缩的VCF文件（VCF文件是常见的突变结果存储格式）：
$ samtools mpileup -v --no-BAQ --fasta-ref human_g1k_v37.fasta --region 1:215906528-215906567 NA12891_CEU_sample.bam >
NA12891_CEU_sample.vcf.gz
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
此VCF文件的主要内容为：
$ zgrep -v "^##" NA12891_CEU_sample.vcf.gz | \
awk 'BEGIN {OFS="\t"}{split($8, a, ";"); print $1,$2,$3,$4,$5,$6,a[1],$9,$10}'
#CHROM POS ID REF ALT QUAL INFO FORMAT NA12891
1 215906528 . G <*> 0 DP=21 PL 0,63,236
...
1 215906547 . C G,<*> 0 DP=22 PL 123,0,103,144,127,233
1 215906548 . G C,<*> 0 DP=22 PL 23,0,163,68,175,207
1 215906549 . G <*> 0 DP=22 PL 0,66,248
1 215906550 . G <*> 0 DP=22 PL 0,63,247
1 215906551 . G <*> 0 DP=22 PL 0,54,242
1 215906552 . G <*> 0 DP=21 PL 0,57,239
1 215906553 . G <*> 0 DP=20 PL 0,48,222
1 215906554 . C G,A,<*> 0 DP=18 PL 0,25,192,27,192,195,39,195,198,198
1 215906555 . G A,<*> 0 DP=19 PL 184,7,0,190,42,204
1 215906556 . G <*> 0 DP=18 PL 0,51,252
...
ALT列值为“<*>”代表该位点没有可靠的突变，当位点除了“<*>”外还有其它碱基的出现（例如215906547，215906548）代表可能存在碱基突变。
使用bcftools call工具我们可以推断并保留可信的突变位点：
$ bcftools call -v -m NA12891_CEU_sample.vcf.gz | \
grep -v "^##" | \
awk 'BEGIN{OFS="\t"}{split($8,a,";");print $1,$2,$3,$4,$5,$6,a[1],$9,$10}'
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
#CHROM POS ID REF ALT QUAL INFO FORMAT NA12891
1 215906547 . C G 90 DP=22 GT:PL 0/1:123,0,103
1 215906555 . G A 157 DP=19 GT:PL 1/1:184,7,0
其中参数-v代表只保留变异位点，-m代表使用multiallelic caller。这里只有两个变异被最终认为是可信的，QUAL列给出变异的可靠性得分，其计算方式采用Phred公式（见碱基质量得分）。如果ALT列为“.”代表此碱基没有变异，此时QUAL列则为此列无突变的可靠性。
去掉-v参数可以保留所有结果：
$ bcftools call -m NA12891_CEU_sample.vcf.gz | grep -v "^##" | \
awk 'BEGIN{OFS="\t"}{split($8,a,";");print $1,$2,$3,$4,$5,$6,a[1],$9,$10}' > NA12891_CEU_sample_calls.vcf
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
$ cat NA12891_CEU_sample_calls.vcf
#CHROM POS ID REF ALT QUAL INFO FORMAT NA12891
1 215906528 . G . 266 DP=21 GT 0/0
...
1 215906547 . C G 90 DP=22 GT:PL 0/1:123,0,103
1 215906548 . G . 11.796 DP=22 GT 0/0
1 215906549 . G . 278 DP=22 GT 0/0
1 215906550 . G . 277 DP=22 GT 0/0
1 215906551 . G . 272 DP=22 GT 0/0
1 215906552 . G . 269 DP=21 GT 0/0
1 215906553 . G . 252 DP=20 GT 0/0
1 215906554 . C . 26.9824 DP=18 GT 0/0
1 215906555 . G A 157 DP=19 GT:PL 1/1:184,7,0
...
215906548位点的可靠性得分只有11.796，意味着这个位置变异可信的概率只有。
FORMAT列包含大量的信息，这里我们只保留了由“;”分隔的第一组信息，即经常出现GT:PL信息，它是什么意思呢？通过检索VCF的头部注释得到相应的解释：
$ bcftools call -m NA12891_CEU_sample.vcf.gz > NA12891_CEU_sample_calls.vcf
$ grep "FORMAT=<ID=GT" NA12891_CEU_sample_calls.vcf
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
$ grep "FORMAT=<ID=PL" NA12891_CEU_sample_calls.vcf
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
根据解释GT代表基因型而PL代表此基因型的可能性得分（同样使用Phred标准）。作为二倍体，基因型存在如下可能性：ref/ref(0/0)，ref/alt(0/1)，alt/ref(1/0)，alt/alt(1/1)。
最后补充一下，由于前面samtools mpileup命令使用了参数"--no-BAQ"表明不关注比对质量的问题，如果我们去掉这个参数的话结果如下：
$ samtools mpileup -v -u --fasta-ref human_g1k_v37.fasta --region 1:215906528-215906567 NA12891_CEU_sample.bam |\
grep -v "^##" | \
awk 'BEGIN{OFS="\t"}{split($8, a, ";");print $1,$2,$3,$4,$5,$6,a[1],$9,$10}'
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
#CHROM POS ID REF ALT QUAL INFO FORMAT NA12891
1 215906528 . G <*> 0 DP=21 PL 0,63,236
...
1 215906547 . C <*> 0 DP=22 PL 0,21,141
1 215906548 . G <*> 0 DP=22 PL 0,42,200
1 215906549 . G <*> 0 DP=22 PL 0,45,221
1 215906550 . G <*> 0 DP=22 PL 0,48,228
1 215906551 . G <*> 0 DP=22 PL 0,45,229
1 215906552 . G <*> 0 DP=21 PL 0,45,231
1 215906553 . G <*> 0 DP=20 PL 0,45,220
1 215906554 . C A,<*> 0 DP=18 PL 0,27,199,39,202,203
1 215906555 . G A,<*> 0 DP=19 PL 194,36,0,194,36,194
...
可以看到215906547，215906548位置之前可能存在的变异更可能是由错配引起从而被直接去掉了。
