今天是生信星球陪你的第237天
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
豆豆写于2019.1.2
今天花花一直在全力备课，豆豆也紧张地学了一天，第一次接触bcftools，真的很强大，感觉陷入其中不能自拔。作为GATK前期训练，bcftools是个很好的入门工具，可以带领我们去了解VCF的更多内容，因此还需要继续巩固熟练
VCF一般是利用'variant caller'或者'SNP caller'的工具对BAM比对文件操作得到的
之前要生成VCF或者它的二进制BCF，需要用samtools mpileup，然后再利用bcftools去进行SNP calling。但是有一个问题就是：samtools与bcftools更新速度都很快，使用mpileup+bcftools call pipeline时会出现版本冲突导致报错的问题，于是后来直接一步到位将mpileup加入了bcftools的功能中
【关于bcftools的介绍，继续向下看】image.png
## 下载bowtie2
cd ~/test
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip
unzip bowtie2-2.3.4.3-linux-x86_64.zip
cd bowtie2-2.3.4.3-linux-x86_64/example/reads
## bowtie2 比对 + samtools排序
wkd=~/test/bowtie2-2.3.4.3-linux-x86_64
$wkd/bowtie2 -x $wkd/example/index/lambda_virus -1 reads_1.fq -2 reads_2.fq | samtools sort -@ 5 -o lamda.bam -
bcftools mpileup -f $wkd/example/reference/lambda_virus.fa lamda.bam |bcftools call -mv -o lamda.call.vcf
# 其中bcftools call使用时需要选择算法，有两个选项：-c(consensus-caller)；-m(multiallelic-caller)，其中-m比较适合多个allel与罕见变异的calling情况
# 另外-v的意思是说：只输出变异位点信息，如果位点不是SNP/InDel就不输出
# From https://www.biostars.org/p/307422/
1. trim reads
2. bwa mem align to genome
3. mark duplicates (e.g. picardtools =》MarkDuplicates)
4. use HaplotypeCaller to generate gvcf
5. CombineGVCFs
6. GenotypeGVCFs on the combined gvcf
7. filter your vcf however you want
8. You can do base recalibration iteratively now if you want with the filtered vcf.
参考：http://www.bio-info-trainee.com/3577.html
bcftools 说明书 https://samtools.github.io/bcftools/bcftools.html
拓展阅读：GATK的深度学习Convolutional Neural Networks（CNN）vs. Google DeepVariant vs. GATK VQSR https://www.biostars.org/p/301751/ [其中 Andrew认为GATK是寻找变异的金标准，但是程序有点繁琐，好用不好入门。有推荐使用samtools+bcftools的，并且它们在鉴定SNVs有优势，但InDel方面就欠缺一些]
https://gatkforums.broadinstitute.org/gatk/discussion/10996/deep-learning-in-gatk4
得到了VCF文件只是第一步，下面还需要一定的技巧获得想要的数据
bcftools与samtools同出李恒之手，是用来操作VCF文件的，之前有一个软件叫vcftools，但是就更新到2015年，bcftools是利用C语言开发，因此处理速度很快，可以接替vcftools。用过samtools的都知道，主命令下还有许多的子命令，例如：samtools index 、samtools sort、samtools flagstat等等。bcftools也是如此，主要有三大功能：
对VCF/BCF构建索引：bcftools index
对VCF/BCF进行操作（查看、排序、过滤、交集、格式转换等）annotate、concat、convert、isec、merge、norm、plugin 、query、reheader、sort、view
找变异 call、consensus、cnv、csq、filter、gtcheck、mpileup、roh、stats
# 测试文件是19号染色体：400-500kb
wget http://data.biostarhandbook.com/variant/subset_hg19.vcf.gz
gunzip subset_hg19.vcf.gz # 先解压，一会做个错误示范
下面👇重点来啦，学习一天得到的干货们！
# 如果我们使用上面解压的subset_hg19.vcf，结果会如何呢？
$ bcftools index subset_hg19.vcf
# 报错！
index: the file is not BGZF compressed, cannot index: subset_hg19.vcf
# 因此我们需要按它的要求来压缩
$ bgzip subset_hg19.vcf
$ bcftools index subset_hg19.vcf.gz #默认是构建.csi索引
# 另外还有一种.tbi索引
$ bcftools index -t subset_hg19.vcf.gz
## 查看
# 同samtools一样，要想查看二进制文件(这里的.bcf)，需要用view
$ bcftools view subset_hg19.bcf
# .vcf可以直接less查看
less -S subset_hg19.vcf.gz
##筛选
# 想要筛选某些样本的VCF信息(多个样本逗号分隔)
$ bcftools view subset_hg19.vcf.gz -s HG00115,HG00116 -o subset.vcf
# 如果不想要某些样本，只需要在样本名前加^
$ bcftools view subset_hg19.vcf.gz -s ^HG00115 -o subset_rm_115.vcf
##过滤
# -k参数得到已知的突变位点(ID列中不是.的那部分)
$ bcftools view subset_hg19.vcf.gz -k -o knowm.vcf
# -n参数得到位置的突变位点(可能是novel新的位点，也就是ID列中为.的部分)
$ bcftools view subset_hg19.vcf.gz -n -o unknowm.vcf
# -c参数设置最小的allele count(AC)值进行过滤(就是说某个变异位点出现了多少次)
$ bcftools view -c 5 subset_hg19.vcf.gz | grep -v '#'
## 简单的格式转换可以用view
# 将vcf转换成bcf：-O指定输出文件类型(u表示未压缩的bcf文件；b表示压缩的bcf文件；v表示未压缩的vcf文件；z表示压缩的vcf文件)
$ bcftools view subset_hg19.vcf.gz -O u -o subset_hg19.bcf
## 想定制转换后的格式，可以用query[就是从VCF/BCF中抽取出某些部分]
# 【-i:得到指定表达式的位点；-e：排除指定表达式的位点；-f：指定输出的format；-H输出表头（也就是format信息的汇总）】
##############################################################################
# 例1：想得到染色体号、变异位置、参考碱基、第一个变异碱基
$ bcftools query -f '%CHROM %POS %REF %ALT{0}\n' subset_hg19.vcf.gz
#【结果】
19 499924 T C
##############################################################################
# 例2：在例1的基础上将空格改成tab分隔，再增加样本名和基因型信息
$ bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' subset_hg19.vcf.gz
#【结果】
19 499924 T C HG00115=0|0 HG00116=0|0 HG00117=0|0 HG00118=0|0 HG00119=0|0 HG00120=0|0
##############################################################################
# 例3：做一个BED文件：chr, pos (0-based), end pos (1-based), id
$ bcftools query -f'%CHROM\t%POS0\t%END\t%ID\n' subset_hg19.vcf
#【结果】
19 499923 499924 rs117528364
##############################################################################
# 例4：将所有变异的样本以及变异情况输出
$ bcftools query -f'[%CHROM:%POS %SAMPLE %GT\n]' -i'GT='alt'' subset_hg19.vcf
#【结果】
19:499978 HG00115 0|1
19:498212 HG00118 1|1
##############################################################################
# 例5：练习-i和-H
$ bcftools query -i'GT='het'' -f'[%CHROM:%POS %SAMPLE %GT \n]' -H subset_hg19.vcf
#【结果有省略】
[1]HG00115:CHROM:[2]HG00115:POS [3]SAMPLE [4]HG00115:GT
[5]HG00116:CHROM:[6]HG00116:POS [7]SAMPLE [8]HG00116:GT
19:400666 HG00115 1|0
19:400666 HG00116 0|1
##############################################################################
# 例6：列出VCF中所有样本
$ bcftools query -l subset_hg19.vcf.gz
##############################################################################
# 例7：提取指定区域中所有的变异信息【一个region用-r；对于多个region（如BED文件）要找变异用-R】
$ bcftools query -r '19:400300-400800' -f '%CHROM\t%POS\t%REF\t%ALT\n' subset_hg19.vcf.gz | head
#【结果有省略】
19 400410 CA C
19 400666 G C
19 400742 C T
##############################################################################
# 例8：相反，如果想去掉某个区域的变异，要用-t 【-t与-r相似，都是指定一个区域，不同的是：-r跳转到那个区域，而-t是扔掉这个区域；并且-t还需要加上^】
$ bcftools query -t ^'19:400300-400800' -f '%CHROM\t%POS\t%REF\t%ALT\n' subset_hg19.vcf.gz | head
#【结果有省略】
19 400819 C G
19 400908 G T
19 400926 C T
# 思考：【如果要去掉多个区域，可以把它们放在BED文件中，然后用-T和^；有没有和例7中的-R很像？】
$ cat >exclude.bed
19 400300 400700
19 400900 401000
$ bcftools query -T ^exclude.bed -f '%CHROM\t%POS\t%REF\t%ALT\n' subset_hg19.vcf.gz | head
#【结果有省略】
19 400742 C T
19 401076 C G
19 401077 G A,T
##############################################################################
# 例9：要找到所有样本中的变异位点(也就是排除-e GT为.或者为0|0)
$ bcftools view -e 'GT='.' | GT='0|0'' subset_hg19.vcf.gz | bcftools query -f '%POS[\t%GT\t]\n' | head
# 这里大括号[]的意思是：对其中的每个sample进行循环处理
#【结果有省略】
402556 0|1 0|1 1|1 1|0 0|1 1|1
402707 0|1 0|1 1|1 1|0 0|1 1|1
# TIP：如果想看是否有样本是纯合/杂合，可以用-g快速看一下【-g可以选择hom(homozygous), het(heterozygous ),miss(missing),还可以连用^表示反选】
$ bcftools view -g het subset_hg19.vcf.gz | bcftools query -f '%POS[\t%GT\t]\n' | head
# 因此，要选择所有样本都是纯合的基因型，就要排除掉het和miss
$ bcftools view -g ^het subset_hg19.vcf.gz | bcftools view -g ^miss | bcftools query -f '%POS[\t%GT\t]\n' | head
##############################################################################
# 例10：只选择InDel (顺便统计下有多少个)
$ bcftools view -v indels subset_hg19.vcf.gz | bcftools query -f '%POS\t%TYPE\n' | wc -l
# -v参数可选snps|indels|mnps|other，多选用逗号分隔；排除某个类型的变异，用-V
##############################################################################
# 例11：基于site quality(QUAL)和read depth(DP)选变异位点
$ bcftools query -i 'QUAL>50 && DP>5000' -f '%POS\t%QUAL\t%DP\n' subset_hg19.vcf.gz | head
#【结果有省略】
400410 100 7773
400666 100 8445
400742 100 15699
用于将单个样本的VCF文件合并成一个多样本的（输入文件需要是bgzip压缩并且有.tbi索引）
# 假设有许多vcf的文件
$ cat sample.list
sp1.vcf.gz
sp2.vcf.gz
sp3.vcf.gz
...
$ bcftools merge -l sample.list > multi-sample.vcf
## 例1：假设现在有三个样本，每个样本有自己的VCF，我们现在想找它们三者之间的共同点(-p 指定输出文件的存放目录前缀；-n指定多少样本)
$ bcftools isec -p tmp -n=3 sp1.vcf.gz sp2.vcf.gz sp3.vcf.gz
## 例2：取至少两个样本中一致的变异位点
$ bcftools isec -p tmp -n+2 sp1.vcf.gz sp2.vcf.gz sp3.vcf.gz
## 例3：只选出一个样本中有，其他样本中没有的变异位点(找到距离-C最近的那个样本中特异的位点)
$ bcftools isec -p tmp -C sp1.vcf.gz sp2.vcf.gz sp3.vcf.gz
# 这样得出来的结果就是：只在sp1中存在，sp2和sp3中都没有
# 例如要统计SNP信息(包括)
$ bcftools view -v snps subset_hg19.vcf.gz > snps.vcf
$ bcftools stats snps.vcf > snps.stats
可视化需要安装latex、matplotlib，直接上conda，然后使用plot-vcfstats
$ plot-vcfstats snps.stats -p tmp/
# -p 指定文件夹名字，结果生成的pdf文件就存放其中。然后主要看summary.pdf
参考：什么资料都比不上manual好！
bcftools的manual：https://samtools.github.io/bcftools/bcftools.html
samtools的manual：http://www.htslib.org/doc/samtools.html
bcftools examples: https://samtools.github.io/bcftools/howtos/query.html
浅显地认识一下GATK，强大的工具需要足够的知识积淀来消化
我们想要找到可靠的结果，一般会采用两种以上的软件进行分析，对于找变异也是如此。想象几个场景：我们用不同的软件对同一个样本进行分析，然后想找到它们共同包含的变异信息；对于不同的样本，我们想知道它们有哪些相同变异，又有哪些不同变异；我们自己的样本中有没有比较特殊的变异…
GATK的这个功能就是用来选择VCF文件中符合指定条件的变异，它的基本用法：
# 输入参数
-V : 原始的VCF文件 [Required参数]
-select ： 指定特定条件(e.g. AF<0.25; dp>1000)
-sn/sf/se : 指定筛选的样本名称/文件/正则条件
# 输出参数
-o : 筛选后的VCF文件
# 筛选设置
-ef : 去掉不符合筛选条件的变异
-env : 去掉没有变异的行
-selectType : 选择特定种类的变异(SNP, InDel ...)
-xl_{sn/sf/se}: 去掉特定的样本
输入参数中，sn（sampleName）表示单个样本名，sf(sampleFile)表示包含多个样本名的文件，se(sampleExpressions)表示用正则表示出样本名。
需要注意的是：SelectVariants功能只是标出了符合条件的变异，并不会默认去除不符合条件的变异。
如果只想保存符合条件的变异，使用ef(excludeFilteredVariants) ；
如果从多个样本中选取变异信息，但其实它们都存在共同的没有变异的位点，就可以用env(excludeNonVariants)去除；
利用ef与env，可以去除不是这个人种的变异位点
例如：我们想看看自己的样本中有多少变异存在于1000Genome中
java -jar GATK.jar \
-R reference.fasta \
-T SelectVariants \
-V 1000G_ALL_Genotypes.vcf.gz \
-sf samples.list \ # 将文件名都存在这个list中
-ef \ #--excludeFiltered 保留PASS与.的位点
-env \ #去除样本中没有在1000Genome中出现的位点
-o 1000G_.vcf \
# 最后生成的文件中只包含过滤没有fail并且至少一个样本中有变异的位点
更详细的参数说明参考官方：https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php
看到这里累了吧，送你个小福利：2018.11最新的GATK Worksheet 链接：https://share.weiyun.com/5vefLMG 密码：yf5322
上面是说提取符合特定条件的子集，许多时候还需要合并多个VCF文件。当合并多个文件时，如果包含同一个位点，但位点信息不一样，就需要考虑设置合并参数。基本用法：
# 输入参数
-V : 多个vcf文件(多个文件可以添加tag描述，用于区分)
# 可选参数
-genotypeMergeOptions ：#关于合并同一个样本的不同基因型
-filteredAreUncalled ： #忽略filtered变异位点
-filteredRecordsMergeType ： #如何对待filter status不一样的位点
-priority #对于多个vcf文件，设置它们的基因型优先级(需要配合之前添加-V时的tag描述，逗号分隔，来表示不同的文件)
# 输出参数
-o :
genotypeMergeOptions的设置：如果目的是Merge two separate callsets，就可以用UNIQUIFY；如果目的是Get the union of calls made on the same samples，就可以用PRIORITIZE
filteredRecordsMergeType的设置：KEEP_IF_ANY_UNFILTERED是至少有一个filter status不一致就不合并这个位点；KEEP_IF_ALL_UNFILTERED是当所有的status不一致才不合并这个位点；KEEP_UNCONDITIONAL不管filter status，只要是位点信息，就合并
例如，我们有两个vcf文件需要合并
java -jar GATK.jar \
-R reference.fa \
-T CombineVariants \
-V:hua file1.vcf \ # 如果文件比较多，可以写到list中，然后只用-V file.list
-V:dou file2.vcf \
-filteredAreUncalled \
-o output.vcf
# 结果文件中，包含了两个输入文件中的变异，每个变异位点都有一个'set=='栏，标识位点来源，可以来自file1可以来自file2，或者来自二者
官方说明：https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_CombineVariants.php
例如，我们想找到某个病人特有变异
java -jar GATK.jar \
-R reference.fa \
-T CombineVariants \
-V:1000G 1000G_ALL_Sites.vcf \ #ESP6500 SNP文件
-V:ESP ESP6500.snps.vcf \ #1000G数据库文件
-V:D2D D2D.vcf \ #假设使我们得到的vcf文件
-o patient.D2D.combined.vcf # 合并的文件
java -jar GATK.jar \
-R reference.fa \
-T SelectVariants \
-V patient.D2D.combined.vcf \ # 加载合并好的vcf文件
-select 'set==D2D' \ # 选择标记的D2D位点
-o patient.D2D.private.vcf #仅在该病人中有，在ESP6500和1000G中没有的变异位点
今天是了解VCF的第二天，今天的主要内容就是关于bcftools的了，自己可以找一些VCF练习题做一下，希望对vcf了解地越来越深，加油⛽️🤓
联系客服
