什么是基因組變異
基因組變異是一個(gè)定義比較模糊的概念. 所謂的變異是相對(duì)于一個(gè)完美的“參考基因組”而言萝映。但是其實(shí)完美的“參考基因組”并不存在，因?yàn)槲覀冎皇沁x擇某一個(gè)物種里的其中似乎比較正常的個(gè)體進(jìn)行測序組裝痢法，然后基于它進(jìn)行后續(xù)的研究萨脑。簡單的說究西，參考和變異是相對(duì)而言，變異也可能完全正常图云。
常見的基因組變異一般可以歸為如下幾類：
- SNP惯悠， 單核苷酸多態(tài)性, 一個(gè)堿基的變化
- INDEL，插入或缺失, 一個(gè)堿基的增加或移除
- SNV竣况， 單核苷酸變異克婶，一個(gè)堿基的改變，可以是SNP，也可以是INDEL
- MNP情萤， 多核苷酸多態(tài)性鸭蛙，一個(gè)區(qū)塊中有多個(gè)保守的SNP
- MNV，多核苷酸變異筋岛，一個(gè)區(qū)塊中有多個(gè)SNP或INDEL
- short variations, 小于50bp的變異
- large-scale variation, 大于50bp的大規(guī)模變異
- SV, 結(jié)構(gòu)變異娶视，通常是上千個(gè)堿基，甚至是染色體級(jí)別上的變異
研究這些變異需要用到不同的手段睁宰，其中普通的DNA二代測序在尋找20bp以下的變異比較靠譜肪获，對(duì)于大于1kb的結(jié)構(gòu)變異而言，采用光學(xué)圖譜(Bionanogenomics)可能更加靠譜一點(diǎn)柒傻。因此孝赫，對(duì)于目前最常用的二代測序而言，還是盡量就找SNP和INDEL吧红符，幾個(gè)堿基的變化找起來還是相對(duì)容易些和靠譜些青柄。
對(duì)于SNP的定義這里也要注意下，最初的SNP的定義指的是單個(gè)堿基導(dǎo)致的多態(tài)性违孝，在群體中廣泛存在(1%)刹前，可用來作為分子標(biāo)記來區(qū)分不同個(gè)體。目前的定義比較粗暴一點(diǎn)雌桑，就是那個(gè)和“參考基因組”不同的單個(gè)位點(diǎn)。值得注意的是祖今，這個(gè)概念可能不同人還有不同的定義校坑，當(dāng)你和別人就某個(gè)問題爭執(zhí)的時(shí)候，最好問問他是如何定義這個(gè)基本概念千诬。由于SNP的廣泛存在耍目，并且變異可能會(huì)導(dǎo)致疾病，也就是存在某些SNP會(huì)導(dǎo)致疾病徐绑。Online Mendelian Inheritance in Man,就是一個(gè)人類遺傳疾病數(shù)據(jù)庫邪驮，建議去看下。
最后說下genotype和haplotype傲茄。genotype毅访，基因型指的是一個(gè)個(gè)體的遺傳組成。但是對(duì)于基因組變異而言盘榨，基因型通常指的是個(gè)體在某個(gè)位點(diǎn)上的等位基因情況喻粹。haplotype, 單倍型最初指的是從單個(gè)親本中遺傳的一組基因，而在基因組變異背景下草巡，則是指一組變異守呜。
一次簡單的變異檢測實(shí)戰(zhàn)
變異檢測(variant calling)即通過比較參考序列和比對(duì)結(jié)果來找到兩者的不同并記錄，基本上可以分為如下幾步：
- 序列比對(duì)
- 比對(duì)后處理(可選)
- 從聯(lián)配中確定變異
- 根據(jù)某些標(biāo)準(zhǔn)進(jìn)行過濾
- 對(duì)過濾的變異注釋
這里面的每一個(gè)可選的工具都有很多，不同工具組合后的分析流程得到的結(jié)果可能會(huì)有很大差異查乒。在變異檢測這一部分目前就有很多軟件弥喉，但是常用并且相對(duì)比較可靠的工具有如下幾個(gè)：
- bcftools: http://www.htslib.org/doc/bcftools.html
- FreeBayes: https://github.com/ekg/freebayes
- GATK: https://software.broadinstitute.org/gatk/
- VarScan2: http://varscan.sourceforge.net/
當(dāng)然這些工具最初都是用于人類基因組。
以埃博拉基因組為例完成一次簡單的Variant Calling玛迄，所需工具為efetch, fastq-dump, emboss/seqret档桃，bwa, samtools和FreeBayes和snpEff。這些都可以通過conda快速安裝憔晒。
第一步： 獲取參考基因組序列藻肄，并建立索引
# 建立文件夾
mkdir -p refs
# 根據(jù)Accession下載
ACC=AF086833
REF=refs/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC | seqret -filter -sid $ACC > $REF
bwa index $REF
第二步： 獲取需要比對(duì)的測序數(shù)據(jù)， 以前10w條為例
# 僅要前10w條read
SRR=SRR1553500
fastq-dump -X 100000 --split-files $SRR
第三步：序列比對(duì)
BAM=$SRR.bam
R1=${SRR}_1.fastq
R2=${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM
第四步：使用freebayes或HaplotypeCaller(GATK4)檢測變異
freebayes -f $REF $BAM > ${SRR}_freebayes.vcf
gatk HaplotypeCaller -R $REF-I $BAM -stand-call-conf 30 \
-bamout bamout.bam
--genotyping-mode DISCOVERY
-O ${SRR}_haplotypecaller.vcf
用IGV可視化的效果如下:
這是最簡單的變異檢測流程拒担，對(duì)于找到的變異還可以進(jìn)一步過濾嘹屯，這一部分內(nèi)容見call variant中關(guān)于snp篩選的一些思考
第五步：變異標(biāo)準(zhǔn)化（可選）
由于VCF文件的靈活性，同一種變異可以通過不同的形式表示, 如下圖
變異標(biāo)準(zhǔn)化按照如下規(guī)則對(duì)變異位點(diǎn)表示進(jìn)行簡化
- 盡可能以少的字符表示變異
- 無等位基因可以標(biāo)識(shí)為長度為0
- 變異位點(diǎn)必須左對(duì)齊
看起來很復(fù)雜从撼，其實(shí)操作起來很簡單
bcftools norm -f $REF SRR1553500_freebayes.vcf > SRR1553500_freebayes_norm.vcf
# Lines total/split/realigned/skipped: 493/0/0/0
大部分軟件州弟，如GATK， freebayes已經(jīng)是標(biāo)準(zhǔn)化的結(jié)果低零。
變異檢測那么簡單嗎婆翔？
經(jīng)過簡單的實(shí)戰(zhàn)之后，似乎變異檢測是一件非常容易的事情掏婶，只要敲幾行命令就行了啃奴。當(dāng)然最開始我也是想的，畢竟無知者無畏雄妥，但是了解的越多最蕾，你就會(huì)發(fā)現(xiàn)事情并沒有那么簡單。**大部分基因組相關(guān)的DNA序列有一些特性是人類的直覺所不能理解的老厌，因?yàn)檫@需要考慮一些背景瘟则。
- DNA序列可以非常的長
- A/T/G/C能夠構(gòu)建任意組合的DNA序列，因此在完全隨機(jī)情況下枝秤，即使隨機(jī)分配也能產(chǎn)生各種各樣的模式醋拧。
- DNA序列只有部分會(huì)受到隨機(jī)影響，基本上這部分序列都是有功能的淀弹。
- 不同物種的不同的DNA序列受到不同的隨機(jī)性影響
- 我們按照實(shí)驗(yàn)流程將大片段DNA破碎成小的部分丹壕，并嘗試通過和參考基因組比對(duì)找到原來的位置。只有它依舊和原來的位置非晨岩常靠近雀费，才能進(jìn)一步尋找變異。
因此即便序列和基因組某個(gè)序列非常接近痊焊，從算法的角度是正確比對(duì)盏袄，但其實(shí)偏離了原來正確的位置忿峻，那么從這部分找到的變異也是錯(cuò)誤的。那我們有辦法解決這個(gè)問題嘛辕羽？基本上不可能逛尚，除非技術(shù)進(jìn)步后，我們可以一次性通讀所有序列刁愿。當(dāng)然目前比較常用的方法是找到最優(yōu)變異绰寞，并且那個(gè)變異能更好的解釋問題，且和每條read中的變異都是一致的铣口。這就是目前變異檢測軟件常用策略：realignment或probabilistic alignment滤钱。
變異注釋
變異注釋意味著猜測遺傳變異(SNP, INDEL, CNY, SV)對(duì)基因功能，轉(zhuǎn)錄本和蛋白序列以及調(diào)控序列的影響脑题。為了對(duì)變異進(jìn)行預(yù)測件缸，預(yù)測軟件需要你提供基因組注釋信息，并且注釋信息的完善程度決定了預(yù)測的準(zhǔn)確性叔遂。變異預(yù)測一般會(huì)提供如下結(jié)果
- 變異位點(diǎn)所在基因組注釋的位置他炊，是轉(zhuǎn)錄本上游，還是編碼區(qū)已艰，還是非編碼RNA
- 列舉出收到影響的轉(zhuǎn)錄本和基因
- 確定變異在蛋白序列上的影響痊末， 如stop_gained(終止密碼子提前), missense(錯(cuò)義), stop_lost(終止密碼子缺失), frameshift(移碼)等
- 對(duì)于人類，還可以和已知的位點(diǎn)進(jìn)行匹配
這些效應(yīng)的定義可以在序列本體論查詢.
變異注釋常用軟件有：VEP, snpEFF, AnnoVar, VAAST2. 其中VEP是網(wǎng)頁工具http://asia.ensembl.org/Tools/VEP, 使用很方便哩掺，可惜支持的物種有限凿叠。snpEFF可以說是支持物種最多的工具，這里使用它疮丛。
snpEff databases > listing.txt
# 確認(rèn)物種名
grep -i ebola databases
# 下載
snpEff download ebola_zaire
下載之后還需要檢查一下snpEFF提供的注釋是否和我們所使用參考基因組一致.
snpEff dump ebola_zaire | less
很不幸的是幔嫂，snpEFF提供的是基于KJ660346構(gòu)建的數(shù)據(jù)庫，而我們使用的是AF086833誊薄。因此需要重新下載對(duì)應(yīng)的參考基因組重新比對(duì)，進(jìn)行注釋锰茉。對(duì)飲參考基因組的地址為https://www.ncbi.nlm.nih.gov/nuccore/KJ660346.1?report=fasta
snpEff ebola_zaire ${SRR}_freebayes.vcf > ${SRR}_annotated.vcf
最終會(huì)生成注釋后的VCF文件以及變異位點(diǎn)的描述性報(bào)告呢蔫。
