劉小澤寫(xiě)于18.8.10着憨，補(bǔ)充于18.8.14-15

這之間經(jīng)歷了第一期授課結(jié)課墩衙，（回中山辦購(gòu)房手續(xù)）遙墻機(jī)場(chǎng)讀了100多頁(yè)生信書(shū)籍务嫡，飛往珠海機(jī)場(chǎng)被臺(tái)風(fēng)滯留（差點(diǎn)回不了中山的家）甲抖，從頭淋到尾打了順風(fēng)車(chē)，前往海陵島的家游泳度假（撿了水晶石）心铃，珠海又要來(lái)16號(hào)臺(tái)風(fēng)“貝碧嘉“（準(zhǔn)備趕快回去）

也許你會(huì)問(wèn)准谚，call是什么鬼？不能通俗易懂點(diǎn)嗎去扣？好的～我的理解是召喚柱衔，BCFtools樊破、gatk等就是你手中的魔法棒，揮一揮唆铐，召喚變異位點(diǎn)

學(xué)一項(xiàng)工具哲戚，掌握一種技術(shù)的目的是為了能夠更好地發(fā)現(xiàn)并解決問(wèn)題，而不是為了跑流程而跑流程艾岂，所有的數(shù)據(jù)分析都要建立在個(gè)人背景知識(shí)儲(chǔ)備和對(duì)實(shí)際項(xiàng)目的認(rèn)知上顺少。因此今天，我也只是用這一個(gè)工具來(lái)做一個(gè)引子王浴，從變異開(kāi)始學(xué)起脆炎。

要學(xué)習(xí)任何事物，都要從“是什么氓辣？為什么秒裕？怎么做？“入手钞啸，那么～

一般來(lái)講几蜻，做一個(gè)全基因組、外顯子組或者簡(jiǎn)化基因組變異分析注釋流程是這樣的：

1. 前期數(shù)據(jù)處理：
   • 預(yù)處理：原始數(shù)據(jù)準(zhǔn)備下載体斩、質(zhì)控過(guò)濾
   • 比對(duì)參考基因組（需要先構(gòu)建索引）
   • 處理得到的比對(duì)bam文件（包括排序入蛆、索引、合并等等）
2. 尋找變異：
   • 利用GATK, bcftools 或 freebayes找到初步的raw variants
   • 對(duì)raw Variants 進(jìn)行篩選
3. 變異文件深度挖掘硕勿、注釋【這里才是分析的精華哨毁，可以說(shuō)是剛剛開(kāi)始】
   • 統(tǒng)計(jì) variant 的各種分布情況和基因型信息
   • 變異注釋
   • 結(jié)合相關(guān)數(shù)據(jù)庫(kù)，尋找特定位點(diǎn)

變異是什么

你我之間源武，只差1%的不同扼褪，而這點(diǎn)不同，卻造就了幾十億人口的千差萬(wàn)別

變異是相對(duì)的粱栖，是在彼此的比較中產(chǎn)生的话浇。我們檢測(cè)基因組變異，都需要參考基因組闹究，對(duì)于人類(lèi)而言幔崖，一般就是用hg19、hg38渣淤、b37版本赏寇。但是參考基因組也是選志愿者測(cè)的，并不完美价认，因此也不能百分之百就說(shuō)參考基因組變異就是0嗅定，所有的都要參考這個(gè)基因組來(lái)確定。

如果只有這一個(gè)條件用踩，不足以讓我們信服渠退。因?yàn)槿绻麢z測(cè)出來(lái)一下假陽(yáng)性的變異位點(diǎn)忙迁，后果可能病急亂投醫(yī)。因此碎乃，在做全基因組檢測(cè)過(guò)程中姊扔，還需要選用其他幾個(gè)大樣本群體的變異位點(diǎn)數(shù)據(jù)庫(kù)進(jìn)行對(duì)照，盡可能排除那些假陽(yáng)性變異梅誓。

基因組變異一般包括：

• 單堿基變異旱眯，學(xué)名單核苷酸多態(tài)性（SNP），原來(lái)的定義是單個(gè)堿基導(dǎo)致的群體中廣泛存在的（約1%）的多態(tài)性证九，后來(lái)指與參考基因組不同的位點(diǎn)删豺，它最常見(jiàn)也最簡(jiǎn)單。正常人的全基因組測(cè)序結(jié)果中大概有幾百萬(wàn)個(gè)SNP愧怜，外顯子組中也存在數(shù)萬(wàn)個(gè)SNP呀页。對(duì)于人來(lái)說(shuō)，SNP之間的平均距離在1.2M左右拥坛，但部分SNP位點(diǎn)僅間隔數(shù)個(gè)堿基甚至相鄰蓬蝶。SNP一般可以分為：
  • 發(fā)生在編碼區(qū)的SNP，由于密碼子具有簡(jiǎn)并性不一定會(huì)引起氨基酸的改變：引起氨基酸變化的叫做 Non-Synonymous SNP猜惋，不引起改變的叫做Synonymous SNP丸氛。
  • 如果氨基酸發(fā)生了改變，又有兩種情況：氨基酸的密碼子變成另一種著摔，因而導(dǎo)致多肽鏈的氨基酸種類(lèi)和順序發(fā)生改變缓窜，這就是錯(cuò)義突變。如果突變導(dǎo)致編碼氨基酸的密碼子變成了終止子谍咆，蛋白質(zhì)合成進(jìn)行到該突變位點(diǎn)時(shí)會(huì)提前終止禾锤，就導(dǎo)致了無(wú)義突變
• 小片段的插入與缺失（合稱InDel），一般發(fā)生在基因組上短的有序的基因片段摹察，長(zhǎng)度小于50bp
• 更大范圍的結(jié)構(gòu)性變異（SV）恩掷，研究的熱門(mén)，長(zhǎng)度大于50bp的片段的插入供嚎、缺失（Big Indel）黄娘，染色體倒位（Inversion）、染色體之間或者內(nèi)部發(fā)生易位（Translocation）克滴、拷貝數(shù)變異（CNV）逼争、串聯(lián)重復(fù)（Tandem repeat）、嵌合體（chimera）

為什么要找變異

• 人類(lèi)基因組的變異與人類(lèi)起源發(fā)展偿曙、疾病防控等密切相關(guān)氮凝，二代測(cè)序的低成本帶來(lái)了很大的便利羔巢，讓我們有能力去做這件事望忆；

• 變異不等于突變罩阵，有的變異缺失能夠引起疾病，比如75%的癌癥患者都存在TP53基因突變;TP53變異與近一半以上的癌癥發(fā)生有關(guān)启摄，包括肺癌稿壁、胃癌、肝癌歉备、膀胱癌傅是、食管癌、乳腺癌蕾羊、宮頸癌等多種癌癥喧笔。但是還有一些基因，比如球場(chǎng)彎道超車(chē)的球員就是11號(hào)染色體上的ACTN3基因上MA000028位點(diǎn)發(fā)生變異龟再。找變異书闸，可以解釋更多生物現(xiàn)象

  img

• 結(jié)構(gòu)變異的尋找更為重要，因?yàn)樗袝r(shí)會(huì)與環(huán)境互作增加出生缺陷利凑、癌癥的風(fēng)險(xiǎn)

如何尋找變異

將高通量數(shù)據(jù)在某個(gè)位點(diǎn)上的堿基于與概率統(tǒng)計(jì)相結(jié)合進(jìn)行檢驗(yàn)

通過(guò)全基因組重測(cè)序以及簡(jiǎn)化基因組測(cè)序(GBS浆劲，RAD-seq)得到變異數(shù)據(jù)

尋找SNP、InDel變異一般流程

1. 將reads比對(duì)到參考基因組
2. 矯正比對(duì)（可選）
3. 從比對(duì)結(jié)果中確定變異位點(diǎn)
4. 根據(jù)某些需求過(guò)濾
5. 變異位點(diǎn)注釋

尋找SV變異的幾種方法：

• 從頭組裝好哀澈，不管質(zhì)量牌借，再與參考基因組比對(duì)
• 根據(jù)reads覆蓋深度（一般用于CNV尋找）：先假設(shè)reads比對(duì)到基因組得到的覆蓋深度符合泊松分布，再檢測(cè)與該分布不一致的部分
• 配對(duì)法：先將雙端reads比對(duì)到參考基因組割按，然后對(duì)兩端reads的方向和距離信息是否與建庫(kù)結(jié)果一致做判斷
• 分割reads法：先確保雙端測(cè)序的一個(gè)read 可以完整的比對(duì)膨报，然后把另一條read進(jìn)行拆分進(jìn)行比對(duì)

最常用的“魔法棒”

• bcftoolshttp://www.htslib.org/doc/bcftools.html

• FreeBayes: https://github.com/ekg/freebayes

  freebayes 的用法比較簡(jiǎn)單，但是不支持多線程适荣，相對(duì)較慢丙躏。輸入文件為和samtoosl 格式一樣，提前去重復(fù)束凑、分組信息做好

• GATK: https://software.broadinstitute.org/gatk/

  尋找變異使用HaplotypeCaller晒旅，它利用實(shí)時(shí)de novo對(duì)有可能是變異的區(qū)域進(jìn)行局部組裝，尋找SNP和InDel

• VarScan2: http://varscan.sourceforge.net/

• 當(dāng)然還有根據(jù)體細(xì)胞汪诉、生殖細(xì)胞废恋、家系數(shù)據(jù)優(yōu)化的軟件，比如GATK就針對(duì)體細(xì)胞扒寄、生殖細(xì)胞做了不同的流程鱼鼓，目前算得上是找變異的黃金搭檔，例如檢測(cè)體細(xì)胞突變就可以用Mutect2

實(shí)例分析

實(shí)例分析

根據(jù)這一篇跑一個(gè)小流程该编，加深理解迄本。唉，現(xiàn)在這學(xué)習(xí)课竣，不跑個(gè)數(shù)據(jù)嘉赎，就像沒(méi)學(xué)習(xí)一樣

第一步 下載數(shù)據(jù)

需要用到EMBOSS工具包置媳，詳細(xì)介紹https://wenku.baidu.com/view/5bdf981c55270722192ef7e3.html?pn=NaN

官網(wǎng)：https://www.ebi.ac.uk/seqdb/confluence/display/THD/EMBOSS+seqret

#軟件安裝Emboss-6.6.0 
cd ~/.biosoft
wget ftp://ftp.ebi.ac.uk/pub/software/unix/EMBOSS/EMBOSS-6.6.0.tar.gz
tar zxvf EMBOSS-6.6.0.tar.gz && cd EMBOSS-6.6.0
chmod -R 755 EMBOSS-6.6.0
./configure
make
make install
#下載參考數(shù)據(jù)
HOME=~/project/ebola
mkdir -p $HOME/ref
cd ~/project/ebola/ref
ACC=AF086833 #ebola病毒在genebank的序列號(hào)
REF=ref/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC |seqret -filter -sid $ACC> $REF
#efetch是根據(jù)genebank accessoin number下載
#Seqret可以讀入一種格式的序列，然后改寫(xiě)成另外一種格式公条。-sid是輸入序列登記號(hào)拇囊，也就是AF086833；-filter是
bwa index $REF

#下載測(cè)試數(shù)據(jù)
mkdir -p $HOME/raw
cd ~/project/ebola/raw
fastq-dump -X 100000 --split-files SRR1553500 #取前10萬(wàn)條
# 27M Aug 10 22:58 SRR1553500_1.fastq
# 27M Aug 10 22:58 SRR1553500_2.fastq

下載測(cè)試數(shù)據(jù)

服務(wù)器100M獨(dú)立光纖還不如自己電腦下載速度快【兩個(gè)同時(shí)下載的】

EMBOSS

或者用一個(gè)小腳本就搞定下載準(zhǔn)備數(shù)據(jù)

mkdir -p ~/project/ebola
ebola=~/project/ebola
mkdir -p $/ebola/ref && cd $ebola/ref
ACC=AF086833
REF=$ebola/ref/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC |seqret -filter -sid $ACC> $REF
bwa index $REF
mkdir -p $ebola/raw && cd $ebola/raw 
fastq-dump -X 100000 --split-files SRR1553500
#-X最大觀測(cè)值（max spot id）意思是截取前多少條

準(zhǔn)備工作

第二步 比對(duì)&召喚變異

長(zhǎng)流程

#一開(kāi)始寫(xiě)好引用靶橱，方便以后
ACC=AF086833
ebola=/vol2/agis/xiaoyutao_group/liuyunze/project/ebola
REF=$ebola/ref/$ACC.fa
SRR=SRR1553500
BAM=$ebola/align/$SRR.bam
R1=$ebola/raw/${SRR}_1.fastq
R2=$ebola/raw/${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
VARI=$ebola/variant

##bwa比對(duì)寥袭，samtools排序并構(gòu)建索引，為了之后更快調(diào)用比對(duì)文件
mkdir -p $ebola/align && cd $ebola/align
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM

mkdir -p $VARI
samtools faidx $REF
##第一種方法：bcftools召喚變異
samtools mpileup -uvf $REF $BAM | bcftools call -vm -Oz > bcftools.vcf.gz
##第二種方法：freebayes
freebayes -f $REF $BAM > $ebola/align/freebayes.vcf
##第三種方法：GATK（版本：4.0.7.0）
#注意：在使用GATK之前关霸，需要先建立兩個(gè)參考基因組的索引文件.dict和.fai【具體參見(jiàn)https://gatkforums.broadinstitute.org/gatk/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference】
#.dict中包含了基因組中contigs的名字传黄，也就是一個(gè)字典；
#構(gòu)建.dict文件（原來(lái)要使用picard的CreateSequenceDictionary模塊队寇，但是現(xiàn)在gatk整合了此模塊尝江，可以直接使用）
gatk CreateSequenceDictionary -R $REF -O $ebola/ref/$ACC.dict
#.fai也就是fasta index file，索引文件英上，可以快速找出參考基因組的堿基
#構(gòu)建
samtools faidx $REF
#gatk開(kāi)始：
#必選 -I -O -R炭序，代表輸入、輸出苍日、參考
#接下來(lái)可以按照字母順序依次寫(xiě)出來(lái)惭聂，這樣比較清晰
#-bamout：將一整套經(jīng)過(guò)gatk程序重新組裝的單倍體基因型（haplotypes）輸出到文件
#-stand-call-conf :低于這個(gè)數(shù)字的變異位點(diǎn)被忽略，可以設(shè)成標(biāo)準(zhǔn)30（默認(rèn)是10）
gatk HaplotypeCaller -R $REF -I $BAM -O $ebola/align/HaplotypeCaller.vcf \
-bamout $ebola/align/$SRR.gatk.bam \
-stand-call-conf 30 
# gatk用時(shí)3.95 minutes.
#<gatk補(bǔ)充>GATK進(jìn)行變異檢測(cè)的時(shí)候相恃，是按照染色體排序順序進(jìn)行的辜纲，并非多條染色體并行檢測(cè)的。因此拦耐，如果樣本數(shù)據(jù)量比較大的話耕腾，一般多個(gè)染色體并行。

比對(duì)結(jié)果

得到的結(jié)果

文件對(duì)比

#gatk找snp
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SRR1553500
AF086833        44      .       G       C       36.77   .       AC=1;AF=0.500;AN=2;BaseQRankSum=0.318;DP=14;ExcessHet=3.0103;FS=2.963;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=2.63;ReadPosRankSum=-1.868;SOR=1.981      GT:AD:DP:GQ:PL  0/1:11,3:14:65:65,0,339

#freebayes找snp
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SRR1553500
AF086833        44      .       GCG     CCG     45.2953 .       AB=0.277778;ABP=10.7311;AC=1;AF=0.5;AN=2;AO=5;CIGAR=1X2M;DP=18;DPB=19.3333;DPRA=0;EPP=13.8677;EPPR=3.87889;GTI=0;LEN=1;MEANALT=3;MQM=60;MQMR=60;NS=1;NUMALT=1;ODDS=10.4296;PAIRED=1;PAIREDR=0.9;PAO=0;PQA=0;PQR=0;PRO=0;QA=179;QR=365;RO=10;RPL=0;RPP=13.8677;RPPR=10.8276;RPR=5;RUN=1;SAF=5;SAP=13.8677;SAR=0;SRF=6;SRP=3.87889;SRR=4;TYPE=snp GT:DP:AD:RO:QR:AO:QA:GL 0/1:18:10,5:10:365:5:179:-11.6057,0,-28.344
#其中SAF 和 SAR 展示了alles的strand信息杀糯，RPL和RPR 展示了reads的方向信息

#bcftools找snp
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SRR1553500
AF086833        44      .       G       GTC     153     .       INDEL;IDV=2;IMF=0.117647;DP=17;VDB=0.709575;SGB=-0.636426;MQSB=1;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=6,1,5,2;MQ=60       GT:PL   0/1:186,0,203

圖形對(duì)比

將得到的三種工具vcf文件扫俺、比對(duì)文件、參考基因組文件都添加進(jìn)IGV固翰，可以比較三種工具的異同點(diǎn)狼纬，多種工具同時(shí)找出的變異位點(diǎn)更具有說(shuō)服力

工具對(duì)比

放大來(lái)看，發(fā)現(xiàn)三種工具在這兩個(gè)位點(diǎn)都發(fā)現(xiàn)了變異骂际，第一個(gè)：參考基因組堿基為T(mén)疗琉，而變異成了C（顯示藍(lán)色）；第二個(gè)：參考基因組為A歉铝，變異成了G（顯示橙色）盈简；另外變異的T是紅色，A是綠色。

不知你發(fā)現(xiàn)沒(méi)有柠贤，大多數(shù)都是T->C香浩，A->G等，都是嘧啶變嘧啶种吸，嘌呤變嘌呤弃衍。這就是轉(zhuǎn)換（transition）呀非；還有一種變異方式坚俗，顛換（transversion），即嘧啶變嘌呤岸裙。相比之下猖败，轉(zhuǎn)換突變比顛換突變更常見(jiàn)，并且與顛換相比在氨基酸序列上產(chǎn)生差異的可能性更低降允。一般Transition 是 Transversion 數(shù)量的2倍恩闻。

轉(zhuǎn)換與顛換

從三種工具得到的變異結(jié)果來(lái)看，會(huì)出現(xiàn)許許多多的變異位點(diǎn)剧董，這些位點(diǎn)是真實(shí)的還是假陽(yáng)性幢尚？對(duì)實(shí)際生物學(xué)有幫助的又是那些位點(diǎn)呢？因此翅楼，有時(shí)需要去繁求簡(jiǎn)尉剩，進(jìn)行過(guò)濾

第三步 原始變異檢測(cè)結(jié)果進(jìn)行過(guò)濾

一般采用gatk或者samtools下的bcftools進(jìn)行過(guò)濾，但他們過(guò)濾的標(biāo)準(zhǔn)不同毅臊，得到的結(jié)果也不同理茎。如果使用GATK，可以采用VQSR（variant recalibration）【注意：非人物種的研究基本不用VQSR】或者直接硬過(guò)濾hard-filtering方式【看名字就知道管嬉，硬過(guò)濾簡(jiǎn)單粗暴皂林，因此能不用就不用它。但是背后的參數(shù)選擇還是要了解】

介紹一下三個(gè)工具的硬過(guò)濾方式：

• gatk的硬過(guò)濾包括VariantFiltration和 SelectVariants兩種方式【VariantFiltration是基于vcf的INFO和FORMAT兩列進(jìn)行蚯撩，而SelectVariants 是從大變異集合中取小子集】

  1. Select Variants【操作要簡(jiǎn)單一些】

     #必備參數(shù)
     -O：輸出文件
     -V：輸入的VCF文件
     #常用可選參數(shù)
     -select-type：選擇需要過(guò)濾的變異類(lèi)型础倍，包括（NO_VARIATION, SNP, MNP, INDEL,SYMBOLIC, MIXED）
     -R：參考基因組
     
     mkdir -p $FLTR && cd $FLTR
     gatk SelectVariants -O gatk_sv_snp.vcf -V $VARI/HaplotypeCaller.vcf \
     -select-type SNP -R $REF
     

  2. Variant Filtration【需要自己指定條件，前提對(duì)VCF比較了解】

     #對(duì)于SNP過(guò)濾【常見(jiàn)標(biāo)準(zhǔn)】
     # QUAL<30.0胎挎，QD<2.0,MQ<40.0, FS>60.0, SOR>3.0, MQRankSum<-12.5, ReadPosRankSum<-8.0著隆，DP>10
     #VariantFiltration只是把不符合的內(nèi)容標(biāo)記出來(lái)，在新文件中呀癣，符合條件的變異位點(diǎn)在FILTER一列會(huì)顯示pass
     

     #對(duì)于InDel過(guò)濾
     #QD<2.0美浦，ReadPosRankSum<-20.0，InbreedingCoeff<-0.8项栏，F(xiàn)S>200.0浦辨，SOR>10.0
     

     QUAL：測(cè)序質(zhì)量值【該位點(diǎn)存在variant的可能性；該值越高，則variant的可能性越大流酬，Q=30時(shí)币厕，錯(cuò)誤率就控制在了0.001】
     QD:Variant Confidence/Quality by Depth 質(zhì)量值/深度【變異質(zhì)量值（Quality）除以覆蓋深度（Depth）得到的比值。這里的變異質(zhì)量值就是VCF中QUAL的值芽腾，它是用來(lái)衡量變異的可靠程度旦装；覆蓋深度是這個(gè)位點(diǎn)上所有含變異位點(diǎn)的樣本的覆蓋深度之和，也就是說(shuō)摊滔，這個(gè)Depth值是通過(guò)將每個(gè)含變異的樣本（除了GT 0/0）的覆蓋深度（VCF中Dp值）累加的結(jié)果】
     【為何采用比值的手段阴绢？原因就是：一般我們得到的數(shù)據(jù)，深度越高的位點(diǎn)質(zhì)量值也越高艰躺，但是肯定會(huì)出現(xiàn)位點(diǎn)測(cè)序深度不一的情況呻袭，有的多測(cè)一些，有的少測(cè)一些腺兴，不管直接使用質(zhì)量值還是深度值左电，都會(huì)有偏差】
     MQ:（RMS Mapping Quality）所有比對(duì)到這個(gè)位點(diǎn)上的read的質(zhì)量值均方根，用于描述比對(duì)質(zhì)量值的離散程度
     FS：（先記住值越小檢測(cè)越嚴(yán)格页响，優(yōu)中取優(yōu)）通過(guò)Fisher檢驗(yàn)對(duì)鏈進(jìn)行偏離度bias檢測(cè)篓足，得到的概率值，越小表明變異準(zhǔn)確度越高【消除正負(fù)鏈的偏好性】闰蚕，對(duì)應(yīng)samtools的DP4【 DP4:高質(zhì)量測(cè)序堿基栈拖，位于REF或ALT前后】如果FS接近于零，就說(shuō)明不會(huì)出現(xiàn)鏈特異性結(jié)果
     SOR：（Strand Odds Ratio）用來(lái)評(píng)估是否存在鏈偏向性陪腌，值越高辱魁，鏈偏向性越大。如果完全沒(méi)有偏差诗鸭，理論等于1
     MQRankSum：（Mapping Quality Rank Sum Test）只針對(duì)雜合子
     ReadPosRankSum：計(jì)算變異的等位基因距離reads尾部的距離染簇，越接近尾部可信度越低
     DP：（Depth）這個(gè)位點(diǎn)過(guò)濾掉一些reads后的覆蓋度

• bcftools:兩種方式 bcftools filter和bcftools view

  #第一步，對(duì)結(jié)果進(jìn)行index
  tabix -p vcf bcftools.vcf.gz #輸入文件一定是使用bgzip壓縮的vcf文件
  #利用bcftools view强岸，提取snp和indel結(jié)果
  bcftools view -i "type='snps'" bcftools.vcf.gz -Oz\
  -o bcftools.snp.vcf.gz && tabix -p vcf bcftools.snp.vcf.gz
  # 提取indel只需要替換那個(gè)type='indels'
  
  #利用bcftools filter
  bcftools filter -e "MQ < 40 || QUAL < 30" -s LOWQUAL \
  bcftools.vcf.gz | bcftools view -f PASS bcftools.snp.vcf.gz
  

• freebayes

  可以參考bcftools锻弓，另外還有其他一些過(guò)濾標(biāo)準(zhǔn)PRR>1 、RPL>1蝌箍、SAF >0青灼、SAR>0、QUAL/AO>10

第四步 變異結(jié)果標(biāo)準(zhǔn)化

標(biāo)準(zhǔn)化就是對(duì)變異位點(diǎn)進(jìn)行簡(jiǎn)化妓盲，用盡量少的字符表示變異杂拨；沒(méi)有等位基因的位點(diǎn)可以標(biāo)記長(zhǎng)度0；變異位點(diǎn)左對(duì)齊悯衬。標(biāo)準(zhǔn)化前后的對(duì)比弹沽，在一篇文獻(xiàn)中有提及變異標(biāo)準(zhǔn)化

變異結(jié)果標(biāo)準(zhǔn)化

使用的GATK、freebayes的結(jié)果都已經(jīng)進(jìn)行過(guò)標(biāo)準(zhǔn)化，對(duì)于bcftools

bcftools norm -f $REF bcftools.vcf  > bcftools_norm.vcf

第五步 變異注釋

對(duì)于人來(lái)講策橘，一般根據(jù)疾病遺傳方式炸渡、變異類(lèi)型（錯(cuò)義missense、synonymous/non-synonymous丽已、終止密碼子提前/缺失stop gain/loss蚌堵、splice、移碼frameshift等）沛婴、種群頻率吼畏、疾病危害性預(yù)測(cè)等條件進(jìn)行過(guò)濾。拿到過(guò)濾后的變異位點(diǎn)信息后瘸味，進(jìn)行遺傳病檢測(cè)宫仗，查看包括SNP够挂、InDel旁仿、CNV、SV等對(duì)基因孽糖、轉(zhuǎn)錄本以及蛋白和調(diào)控元件的影響。因此變異位點(diǎn)注釋信息必須準(zhǔn)確，否則會(huì)遺漏潛或者診斷錯(cuò)誤惜犀。

最常用的四種注釋工具是VEP牙甫、Annovar、snpeff病蛉、 VAAST2炫加。其中snpeff支持物種最多，VEP直接網(wǎng)頁(yè)版使用

歡迎關(guān)注我們的公眾號(hào)～_～　　
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩铺然，打造生信星球俗孝，想讓它成為一個(gè)不拽術(shù)語(yǔ)、通俗易懂的生信知識(shí)平臺(tái)魄健。需要幫助或提出意見(jiàn)請(qǐng)后臺(tái)留言或發(fā)送郵件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!