變異資訊那些事（中）

劉小澤寫於2019.1.2
今天花花一直在全力備課，我也緊張地學了一天，第一次接觸bcftools，真的很強大，感覺陷入其中不能自拔。作為GATK前期訓練，bcftools是個很好的入門工具，可以帶領我們去了解VCF的更多內容，因此還需要繼續鞏固熟練

如何得到VCF？

VCF一般是利用"variant caller"或者"SNP caller"的工具對BAM比對檔案操作得到的

bowtie2+samtools+bcftools

之前要生成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就不輸出

GATK

# 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基本操作

得到了VCF檔案只是第一步，下面還需要一定的技巧獲得想要的資料

利用bcftools

背景

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 # 先解壓，一會做個錯誤示範

下面?重點來啦，學習一天的乾貨們！

構建索引 index =>需要用到bgzip格式

# 如果我們使用上面解壓的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

檢視、篩選、過濾 => view

## 檢視
# 同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 "#"

標準/DIY格式轉換 => view / query

## 簡單的格式轉換可以用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/BCF => merge

用於將單個樣本的VCF檔案合併成一個多樣本的（輸入檔案需要是bgzip壓縮並且有.tbi索引）

# 假設有許多vcf的檔案
$ cat sample.list
sp1.vcf.gz
sp2.vcf.gz
sp3.vcf.gz
...
$ bcftools merge -l sample.list > multi-sample.vcf

取交集 => isec

## 例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中都沒有

方便的統計=>stat

# 例如要統計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，強大的工具需要足夠的知識積澱來消化
我們想要找到可靠的結果，一般會採用兩種以上的軟體進行分析，對於找變異也是如此。

想象幾個場景：我們用不同的軟體對同一個樣本進行分析，然後想找到它們共同包含的變異資訊；對於不同的樣本，我們想知道它們有哪些相同變異，又有哪些不同變異；我們自己的樣本中有沒有比較特殊的變異...

GATK SelectVariants

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

GATK CombineVariants

上面是說提取符合特定條件的子集，許多時候還需要合併多個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

Combine聯合Select

例如，我們想找到某個病人特有變異

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瞭解地越來越深，加油⛽️?

---

歡迎關注我們的公眾號～_～　　
我們是兩個農轉生信的小碩，打造生信星球，想讓它成為一個不拽術語、通俗易懂的生信知識平臺。需要幫助或提出意見請後臺留言或傳送郵件到Bioplanet520@outlook.com